geometry
geometry copied to clipboard
Strange behaviour of delaunayn with options
I'm seeing some hither to unseen behavior with this function. Im using version ‘0.4.7’ of geometry on Windows with R version 4.2.2. The following is a reproducible example.
source("Evaldf.txt") library(geometry) delaunayn(evaldf[,c("XCoord","YCoord")],options="") delaunayn(evaldf[,c("XCoord","YCoord")])
The first call to delaunayn returns a matrix of 0 rows and 3 columns- presumably incorrect. The second call to delaunayn with options returns a matrix of 514 rows and 3 columns- presumably correct.
I inherited this piece of code and I have to admit that I don’t fully understand why options=""
is specified. I suspect it may be legacy related to previous version of the package. 99.9% of the time the first call with options="" works absolutely fine.
It's just with this specific data set it's exhibiting this bug. Can you shed any light on this?
Evaldf.txt
From the output I only need the matrix of Delaunay triangulation. Can you recommend what is the safest option choice to specify when calling the delaunayn function which will also work reliably with older package versions?
I could just remove the options=""
in the code but I'm nervous it will cause issues with older versions of the package.
I tried removing the options="" but I saw a change in behaviour. So I am reverting back to include options="". https://github.com/WayneGitShell/GWSDAT/blob/GWWellReport/R/interpBary.R#L70
See weird behavior of resulting barycentric interpolation following removal of options="".
Having removed options=""
With options=""
I think this issue is related to the problem to a set of points that span a small distance relative to their distance to the origin (Issue #11). This issue is fixed by the update of the Qhull library to 8.0.2 (Issue #37, currently in the qhull-8.0.2 branch). I've just tested this dataset with the code on that branch:
source("Evaldf.txt")
library(geometry)
p <-evaldf[,c("XCoord","YCoord")]
t1 <- delaunayn(p,options="")
t2 <- delaunayn(p)
## Count points missing from triangulation. Should be 0 in both cases - but the first is not with geometry-0.4.7
print(length(setdiff(seq(1,nrow(p)), unique(c(t1[,1], t1[,2], t1[,3])))))
print(length(setdiff(seq(1,nrow(p)), unique(c(t2[,1], t2[,2], t2[,3])))))
In summary, I would recommend requiring the latest version of geometry (as in the master branch here, soon to be 0.5.0) for GWSDAT and not specifying options="". Additionally, you may want to centre the coordinates by subtracting the mean x and y values from the x and y coordinates.
In Issue #12 I added a patch to throw a warning in case that points are missing from a triangulation.
Many thanks for solving this issue and advising on best way forward.