R/util.R

In ads: Spatial Point Patterns Analysis

overlapping.polygons<-function(listpoly) {
	stopifnot(unlist(lapply(listpoly,is.poly)))
	nbpoly<-length(listpoly)
	res<-rep(FALSE,(nbpoly-1))
	for(i in 1:(nbpoly-1)) {
		for(j in (i+1):nbpoly) {
		 if(abs(overlap.poly(listpoly[[i]],listpoly[[j]]))>.Machine$double.eps^0.5)
			res[i]<-TRUE
		}
	}
	return(res)
}

#from area.xypolygon{spatstat.geom}
#return area>0 when (xp,yp) vertices are ranked anticlockwise
#or area<0 when (xp,yp) vertices are ranked clockwise
area.poly<-function(xp,yp) {
    nedges <- length(xp)
    yp <- yp - min(yp)
    nxt <- c(2:nedges, 1)
    dx <- xp[nxt] - xp
    ym <- (yp + yp[nxt])/2
    -sum(dx * ym)
}

in.circle<-function(x,y,x0,y0,r0,bdry=TRUE) {
	stopifnot(length(x)==length(y))
	l<-length(x)
	inside<-vector(mode="logical",length=l)
	for(i in 1:l) {
		if(bdry) {
			if(((x[i]-x0)^2+(y[i]-y0)^2)<=(r0^2))
				inside[i]<-TRUE
		}
		else {
			if(((x[i]-x0)^2+(y[i]-y0)^2)<(r0^2))
				inside[i]<-TRUE
		}
	}
	return(inside)
}

in.rectangle<-function(x,y,xmin,ymin,xmax,ymax,bdry=TRUE) {
	stopifnot(length(x)==length(y))
	stopifnot((xmax-xmin)>0)
	stopifnot((ymax-ymin)>0)
	l<-length(x)
	inside<-vector(mode="logical",length=l)
	for(i in 1:l) {
		if(bdry) {
			if(x[i]>=xmin&&x[i]<=xmax&&y[i]>=ymin&&y[i]<=ymax)
				inside[i]<-TRUE
		}
		else {
			if(x[i]>xmin&&x[i]<xmax&&y[i]>ymin&&y[i]<ymax)
				inside[i]<-TRUE
		}
	}
	return(inside)
}

#modified from inside.triangle(spatstat.utils)
in.triangle<-function(x,y,ax,ay,bx,by,cx,cy,bdry=TRUE) {
	stopifnot(length(x)==length(y))
	p<-length(x)
	inside<-vector(mode="logical",length=p)
	t<-length(ax)
	stopifnot(all(c(length(ay),length(bx),length(by),length(cx),length(cy))==t))
	intri<-matrix(0,p,t)
	for(j in 1:t) {
		v0x <- cx[j] - ax[j]
    	v0y <- cy[j] - ay[j]
    	v1x <- bx[j] - ax[j]
    	v1y <- by[j] - ay[j]
		dot00 <- v0x^2 + v0y^2
    	dot01 <- v0x * v1x + v0y * v1y
		dot11 <- v1x^2 + v1y^2
		for(i in 1:p) {
			if(!inside[i]) {
				v2x <- x[i] - ax[j]
    			v2y <- y[i] - ay[j]
				dot02 <- v0x * v2x + v0y * v2y
    			dot12 <- v1x * v2x + v1y * v2y
    			Denom <- dot00 * dot11 - dot01 * dot01
    			u <- dot11 * dot02 - dot01 * dot12
    			v <- dot00 * dot12 - dot01 * dot02
				if(bdry) {
          			if((u >= 0) & (v >= 0) & (u + v <= Denom))
            			inside[i]<-TRUE
        		}
        		else {
          			if((u > 0) & (v > 0) & (u + v < Denom))
              			inside[i]<-TRUE
        		}
			}
		}	
	}
	return(inside)			
}

#modified from plot.ppp{spatstat.core}
adjust.marks.size<-function(marks,window,maxsize=NULL) {
	if(is.null(maxsize)) {
		if("rectangle"%in%window$type)
			diam<-sqrt((window$xmin-window$xmax)^2+(window$ymin-window$ymax)^2)
		else if("circle"%in%window$type)
			diam<-2*window$r0
		maxsize<-min(1.4/sqrt(pi*length(marks)/area.swin(window)),diam*0.07)
	}
	mr<-range(c(0,marks))
	maxabs<-max(abs(mr))
	if(diff(mr)<4*.Machine$double.eps||maxabs<4*.Machine$double.eps) {
		ms<-rep(0.5*maxsize,length(marks))
		mp.value<-mr[1]
		mp.plotted<-0.5*maxsize
	}
	else {
		scal<-maxsize/maxabs
		ms<-marks*scal
		mp.value<-pretty(mr)
		mp.plotted<-mp.value*scal
	}
	return(ms)
}

#modified from verify.xypolygon{spatstat}
is.poly<-function (p) {
	stopifnot(is.list(p))
	stopifnot(length(p)==2,length(names(p))==2)
	stopifnot(identical(sort(names(p)),c("x","y")))
	stopifnot(!is.null(p$x),!is.null(p$y))
	stopifnot(is.numeric(p$x),is.numeric(p$y))
	stopifnot(length(p$x)==length(p$y))
	return(TRUE)
}

testInteger<-function(i) {
	if(as.integer(i)!=i) {
		warning(paste(substitute(i),"=",i," has been converted to integer : ",as.integer(i),sep=""),call.=FALSE)
		i<-as.integer(i)
	}
	return(i)
}

testIC<-function(nbSimu,lev) {
	if(lev*(nbSimu+1)<5) {
		warning(paste(
			"Low validity test: a*(n+1) < 5\n  Significance level: a = ",lev,
			"\n  Number of simulations: n = ",nbSimu,"\n",sep=""))
	}
}

#from spatstat
overlap.poly <- function(P, Q) {
  # compute area of overlap of two simple closed polygons 
 # verify.xypolygon(P)
  #verify.xypolygon(Q)
  
  xp <- P$x
  yp <- P$y
  np <- length(xp)
  nextp <- c(2:np, 1)

  xq <- Q$x
  yq <- Q$y
  nq <- length(xq)
  nextq <- c(2:nq, 1)

  # adjust y coordinates so all are nonnegative
  ylow <- min(c(yp,yq))
  yp <- yp - ylow
  yq <- yq - ylow

  area <- 0
  for(i in 1:np) {
    ii <- c(i, nextp[i])
    xpii <- xp[ii]
    ypii <- yp[ii]
    for(j in 1:nq) {
      jj <- c(j, nextq[j])
      area <- area +
        overlap.trapez(xpii, ypii, xq[jj], yq[jj])
    }
  }
  return(area)
}

#from spatstat
overlap.trapez <- function(xa, ya, xb, yb, verb=FALSE) {
  # compute area of overlap of two trapezia
  # which have same baseline y = 0
  #
  # first trapezium has vertices
  # (xa[1], 0), (xa[1], ya[1]), (xa[2], ya[2]), (xa[2], 0).
  # Similarly for second trapezium
  
  # Test for vertical edges
  dxa <- diff(xa)
  dxb <- diff(xb)
  if(dxa == 0 || dxb == 0)
    return(0)

  # Order x coordinates, x0 < x1
  if(dxa > 0) {
    signa <- 1
    lefta <- 1
    righta <- 2
    if(verb) cat("A is positive\n")
  } else {
    signa <- -1
    lefta <- 2
    righta <- 1
    if(verb) cat("A is negative\n")
  }
  if(dxb > 0) {
    signb <- 1
    leftb <- 1
    rightb <- 2
    if(verb) cat("B is positive\n")
  } else {
    signb <- -1
    leftb <- 2
    rightb <- 1
    if(verb) cat("B is negative\n")
  }
  signfactor <- signa * signb # actually (-signa) * (-signb)
  if(verb) cat(paste("sign factor =", signfactor, "\n"))

  # Intersect x ranges
  x0 <- max(xa[lefta], xb[leftb])
  x1 <- min(xa[righta], xb[rightb])
  if(x0 >= x1)
    return(0)
  if(verb) {
    cat(paste("Intersection of x ranges: [", x0, ",", x1, "]\n"))
    abline(v=x0, lty=3)
    abline(v=x1, lty=3)
  }

  # Compute associated y coordinates
  slopea <- diff(ya)/diff(xa)
  y0a <- ya[lefta] + slopea * (x0-xa[lefta])
  y1a <- ya[lefta] + slopea * (x1-xa[lefta])
  slopeb <- diff(yb)/diff(xb)
  y0b <- yb[leftb] + slopeb * (x0-xb[leftb])
  y1b <- yb[leftb] + slopeb * (x1-xb[leftb])
  
  # Determine whether upper edges intersect
  # if not, intersection is a single trapezium
  # if so, intersection is a union of two trapezia

  yd0 <- y0b - y0a
  yd1 <- y1b - y1a
  if(yd0 * yd1 >= 0) {
    # edges do not intersect
    areaT <- (x1 - x0) * (min(y1a,y1b) + min(y0a,y0b))/2
    if(verb) cat(paste("Edges do not intersect\n"))
  } else {
    # edges do intersect
    # find intersection
    xint <- x0 + (x1-x0) * abs(yd0/(yd1 - yd0))
    yint <- y0a + slopea * (xint - x0)
    if(verb) {
      cat(paste("Edges intersect at (", xint, ",", yint, ")\n"))
      points(xint, yint, cex=2, pch="O")
    }
    # evaluate left trapezium
    left <- (xint - x0) * (min(y0a, y0b) + yint)/2
    # evaluate right trapezium
    right <- (x1 - xint) * (min(y1a, y1b) + yint)/2
    areaT <- left + right
    if(verb)
      cat(paste("Left area = ", left, ", right=", right, "\n"))    
  }

  # return area of intersection multiplied by signs 
  return(signfactor * areaT)
}

#Points on boundary are considered outside. No alternative option implemented yet.
in.poly<-function(x,y,poly,bdry=FALSE) {
	if(bdry) {
		bdry<-FALSE
		warning("argument 'bdry' automatically set to FALSE. No alternative implemented yet")
	}
	stopifnot(is.poly(poly))
	npts <- length(x)
	nedg <- length(poly$x)   # sic
	xmi<-ifelse(min(x)<=min(poly$x),min(x),min(poly$x))
	ymi<-ifelse(max(x)>=max(poly$x),max(x),max(poly$x))
	polxrange<-ifelse(min(y)<=min(poly$y),min(y),min(poly$y))
	polyrange<-ifelse(max(y)>=max(poly$y),max(y),max(poly$y))
	res <- .C("pnpoly",
 				as.double(x),as.double(y),as.double(poly$x),
 				as.double(poly$y),as.integer(npts),as.integer(nedg),as.double(xmi),
		 		as.double(ymi),as.double(polxrange),as.double(polyrange),score=double(npts),
				PACKAGE="ads")
	return(as.logical(res$score))
}

####################
convert<-function(x) {
r<-alist()
x<-as.matrix(x)
for(i in 1:dim(x)[1])
	r[[i]]<-data.frame(x=c(x[i,1],x[i,5],x[i,3]),y=c(x[i,2],x[i,6],x[i,4]))
return(r)
}

convert2<-function(x){ # liste de liste vers df 6 var
x<-unlist(x)
mat<-matrix(x,ncol=6,byrow=TRUE,dimnames=list(NULL,c("ax","bx","cx","ay","by","cy")))
#mat<-cbind(tmp$ax,tmp$ay,tmp$bx,tmp$by,tmp$cx,tmp$cy)
r<-data.frame(ax=mat[,1],ay=mat[,4],bx=mat[,2],by=mat[,5],cx=mat[,3],cy=mat[,6])
return(r)
}


read.tri<-function(X) {
	res<-NULL
	tabtri<-read.table(X)	
	
	n<-length(tabtri[,1])
	
	for(i in 1:n) {
		tri<-list(x=c(tabtri[,1][i],tabtri[,3][i],tabtri[,5][i]),
				  y=c(tabtri[,2][i],tabtri[,4][i],tabtri[,6][i]))
				  
		#if(area.xypolygon(tri)>0){
		if(area.poly(tri$x,tri$y)>0){
			tri<-list(x=c(tabtri[,5][i],tabtri[,3][i],tabtri[,1][i]),
				  y=c(tabtri[,6][i],tabtri[,4][i],tabtri[,2][i]))
		}
				  
		res<-c(res,list(tri))
	}
	if(length(res)==1) {
		res<-unlist(res,recursive=FALSE)
	}
	return(res)
}

transpose<-function(x,y) {

	nbTri<-length(x)/3
	
	res<-.C("transpose",x=as.double(x),y=as.double(y),nbTri=as.integer(nbTri),
			x1=double(nbTri),y1=double(nbTri),x2=double(nbTri),y2=double(nbTri),
			x3=double(nbTri),y3=double(nbTri),PACKAGE="ads")
		
	list(x1=res$x1,y1=res$y1,x2=res$x2,y2=res$y2,x3=res$x3,y3=res$y3)
}
##############
#subsetting dist objects
#sub is a logical vector of True/False
subsetdist<-function(dis,sub) {
	mat<-as.matrix(dis)
	k<-dimnames(mat)[[1]]%in%sub
	submat<-mat[k,k]
	return(as.dist(submat))
}
	
#ordering dist objetcs on ind
sortmat<-function(dis,ind) {
	mat<-as.matrix(dis)[,ind]
	mat<-mat[ind,]
	return(as.dist(mat))
}