triangulate() does not triangulate a PSLG with two holes correctly

[davidcsterratt]

I've had a bug report from Paolo Piras:

With one hole pslg() and triangulate() properly triangulate the shape, with two holes it is not the case. Here below a fully reproducible example

library(RTriangle)
myshape<-matrix(c
(1507,1370,1470,1598,1685,1551,1470,1490,1779,2249,2517,2652,2598,2370,2209,2162,2323,2517,2611,2544,2370,2276,2081,2054,1994,1893,1745,1611,1605,1887,1794,1913,2226,2306,2140,1998,1851,1799,1856,2016,2081,1988,2333,2126,1938,1723,1515,1233,944,669,535,508,703,1025,1253,1381,1575,1703,1891,1944,2092,2179,2186,2112,2105,2280,2481,2374,2099,2206,2307,972,822,672,716,947,1009,1022,1909,1834,1681,1676,1787,1917),ncol
=2)
plot(myshape,asp=1)
text(myshape[,1],myshape[,2],c(1:nrow(myshape)))
links<-matrix(c
(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,36,30,31,32,33,34,35,42,37,38,39,40,41,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,1,30,31,32,33,34,35,36,37,38,39,40,41,42),ncol
=2)
for(i in 1:nrow(links)){
  segments(myshape[links[i,1],1],myshape[links[i,1],2],myshape[links
[i,2],1],myshape[links[i,2],2])
}
#### Let's try with one hole structure
p<-pslg(myshape,S=links,H=rbind(apply(myshape[30:36,],2,mean)))
pt<-triangulate(p)
plot(pt,asp=1)#### with one hole it works
####  try with two holes
p<-pslg(myshape,S=links,H=rbind(apply(myshape[30:36,],2,mean),apply(myshape
[37:42,],2,mean)))
pt<-triangulate(p)
plot(pt,asp=1)####    why it does not work with two holes?

[epipping]

I was puzzled by the fact that more than one hole wasn't properly triangulated at some point as well. Then I noticed that transposing H fixed the issue for me. If I continue to work with such a transposed H, I now get an error (albeit fortunately a rather helpful one). So even though the previous solution maybe shouldn't have worked, this is actually a bit of a breaking change...

The issue even goes a bit further: The case of no holes now needs to be special cased. I'm generating my list of holes through rgeos::gPointOnSurface(...)@coords inside sapply. To address the aforementioned issue, I now need to transpose the result. But if sapply is applied to an empty list, I'll get an empty list in return (previously an acceptable H), which is unchanged by t(), yet the pslg() constructor will always make sure the number of columns is equal to 2 now while the empty list has no columns.

Maybe it would be a good idea to carry out the number-of-columns check only after H is found to contain anything at all so that I don't need to pass a matrix(ncol=2,nrow=0)? I thought I'd at least suggest it...

[davidcsterratt]

Oh dear, I'd not realised this would break some things. However, the fix does now mean that pslg() is doing what it's documented to do. I'm a little anxious about adding sanitisation of empty lists, as I then feel I might need to document it. Also, perhaps then it would make sense to supply options for giving lists of hole coordinates... However, I'd be open to suggestions for patches.

[mdsumner]

Fwiw I post-process holes by testing if any triangle centroid falls in a hole, it's pretty efficient. This is in r-gris/trimesh in the tri_mesh function. It's not ready for general use yet but getting close.