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