1triang.list <- function (object) 2{ 3 stopifnot(inherits(object,"deldir")) 4 a <- object$delsgs[, 5] 5 b <- object$delsgs[, 6] 6 tlist <- matrix(integer(0), 0, 3) 7 for(i in seq(nrow(object$summary))) { 8 # find all Delaunay neighbours of i 9 jj <- c(b[a==i], a[b==i]) 10 jj <- sort(unique(jj)) 11 # select those with a higher index than i 12 jj <- jj[jj > i] 13 # find pairs of neighbours which are Delaunay neighbours 14 # (thus, triangles where the first numbered vertex is i) 15 if(length(jj) > 0) 16 for(j in jj) { 17 kk <- c(b[a == j], a[b == j]) 18 kk <- kk[(kk %in% jj) & (kk > j)] 19 if(length(kk) > 0) 20 for(k in kk) 21 # add (i,j,k) to list of triangles (i < j < k) 22 tlist <- rbind(tlist, c(i, j, k)) 23 } 24 } 25 x <- object$summary[,"x"] 26 y <- object$summary[,"y"] 27 xtri <- matrix(x[tlist], nrow(tlist), 3) 28 ytri <- matrix(y[tlist], nrow(tlist), 3) 29 ztri <- ytri - min(y) 30 dx <- cbind(xtri[, 2] - xtri[, 1], xtri[, 3] - xtri[, 2], 31 xtri[, 1] - xtri[, 3]) 32 zm <- cbind(ztri[, 1] + ztri[, 2], ztri[, 2] + ztri[, 3], 33 ztri[, 3] + ztri[, 1]) 34 negareas <- apply(dx * zm, 1, sum) 35 clockwise <- (negareas > 0) 36 if (any(clockwise)) { 37 xc <- xtri[clockwise,,drop=FALSE] 38 yc <- ytri[clockwise,,drop=FALSE] 39 xtri[clockwise, ] <- xc[, c(1, 3, 2)] 40 ytri[clockwise, ] <- yc[, c(1, 3, 2)] 41 } 42 rslt <- list() 43 K <- 0 44 for(i in 1:nrow(xtri)) { 45 tmp <- .Fortran( 46 "intri", 47 x=as.double(xtri[i,]), 48 y=as.double(ytri[i,]), 49 u=as.double(x), 50 v=as.double(y), 51 n=as.integer(length(x)), 52 okay=logical(3*length(x)), 53 xxx=double(3*length(x)), 54 PACKAGE="deldir") 55# 56okay <- matrix(tmp$okay,ncol=3) 57xxx <- matrix(tmp$xxx,ncol=3) 58ok <- apply(okay,1,any) 59chk <- apply(xxx,1,function(x){any(x<=0)}) 60if(!isTRUE(all.equal(ok,chk))) browser() 61# 62 if(all(ok)) { 63 K <- K+1 64 rslt[[K]] <- data.frame(x=xtri[i,],y=ytri[i,]) 65 } 66 } 67 attr(rslt,"rw") <- object$rw 68 class(rslt) <- "triang.list" 69 rslt 70} 71