# 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}
#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)
rect<-list(x=c(xmin,xmax,xmax,xmin),y=c(ymin,ymin,ymax,ymax))
return(in.poly(x,y,rect,bdry))
}

in.triangle<-function(x,y,ax,ay,bx,by,cx,cy,bdry=TRUE) {
stopifnot(length(x)==length(y))
tri<-list(x=c(ax,bx,cx),y=c(ay,by,cy))
return(in.poly(x,y,tri,bdry))
}

#modified from plot.ppp{spatstat}
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)
}

#TRUE: les points sur la bordure sont = inside
in.poly<-function(x,y,poly,bdry=TRUE) {
stopifnot(is.poly(poly))
xp <- poly\$x
yp <- poly\$y
npts <- length(x)
nedges <- length(xp)   # sic

score <- rep(0, npts)
on.boundary <- rep(FALSE, npts)
temp <- .Fortran(
"inpoly",
x=as.double(x),
y=as.double(y),
xp=as.double(xp),
yp=as.double(yp),
npts=as.integer(npts),
nedges=as.integer(nedges),
score=as.double(score),
onbndry=as.logical(on.boundary),
)
score <- temp\$score
on.boundary <- temp\$onbndry
score[on.boundary] <- 1
res<-rep(FALSE,npts)
res[score==(-1)]<-TRUE
if(bdry)
res[score==1]<-TRUE
return(res)
}

####################
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)
}

res<-NULL

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)
}

##############
#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))
}
```