VoronoiMesh: user error, tetgen error, else?

[sghuisman (Sander Huisman)]

I found some unusual error in VoronoiMesh. I have the points pts, and first I compute the VoronoiMesh for that, Note the red circle: these are two close datapoints, and it correctly draws a line between them.

dom=1.25\[Pi] (* domain size *)
domspec=Transpose[{{0,0},{dom,dom}}] 
pts={{1.880951464600067`,1.473752017061138`},{1.8913812543459405`,2.319493673670671`},{1.7674367495794008`,3.5325088938533917`},{1.9320143413972677`,3.611151717969884`},{1.8107126970778207`,3.6790935311559876`},{1.804725058474518`,3.8520729911980185`},{1.805368854133445`,3.852728507949241`},{1.563939776616803`,3.9192958794919646`}};
VoronoiMesh[pts,domspec,Epilog->{Point[pts],Red,Circle[pts[[6]],0.1]},PlotRange->domspec,ImageSize->500,PlotRangePadding->Scaled[.05]]

Now, I add points around it (to make the Voronoi periodic), leaving the old ones, and only adding more points around it. Then the voronoi is not correctly computed:

allpts=Join@@Join@@Table[p+dom {i,j},{i,-1,1},{j,-1,1},{p,pts}];  (* add groups of points around it, making it periodic *)
VoronoiMesh[allpts,Epilog->{Point[pts],Red,Circle[pts[[6]],0.1]},PlotRange->domspec,ImageSize->500,PlotRangePadding->Scaled[.05]]

I have millions and millions of points and it took me several hours to ‘minimize’ the error and go to the core of it.

Any help is welcome! Is there some option that controls when two points are considered the same? It looks like some points are deleted because they are too close?

The distance between point 6 and 7 is:

Log10[Norm[pts[[6]] - pts[[7]]]]

-3.03678

So it should not be some machine-precision issue I presume! If someone can confirm it on different versions that would be appreciated.

This is done in 10.4.1 and 10.3.1 on Mac, both show the error.

[sghuisman (Sander Huisman)]

This has as a consequence—very disturbingly—that the number of Voronoi cells does not equal the number of input-points…

[henrikschachner]

Hi Sander,

good observation! I can confirm this on my system too.

$Version
(* "10.3.1 for Linux x86 (32-bit) (December 8, 2015)"  *)

Regards – Henrik

[sghuisman (Sander Huisman)]

Thanks for checking! I think VoronoiMesh is not implemented in the kernel but through Tetgen-library. But apparently also doesn’t work on Linux.

I was driving me mad today! The main problem was that I need the mesh ordered in the original order of the points (by default it isn’t!), and I needed this periodic Voronoi. So I had a lot of data, and one of the files did not work. So I thought it was my custom PeriodicVoronoi function that had an error. I had number of point xyz-million and number of mesh-polygons was xyz-million minus one! But after long digging down I found the error…very frustrating! lost one day on it… oops!

[sukrause (Udo Krause)]

As it stands also defective with Mathematica 10.4.1.0 under Windows 10 Professional 64 Bit

[sghuisman (Sander Huisman)]

@Udo Krause thanks for checking! All 3 operating systems and various versions of Mathematica… I guess a bug…

[sukrause (Udo Krause)]

This might be fixed with a change in dom:

(* comm849935 *)
Clear[dom, domspec, pts]
dom = (* 1.25 \[Pi] *) 3.; (*domain size*)
domspec = Transpose[{{0, 0}, {dom, dom + 1}}];
pts = {{1.880951464600067`, 1.473752017061138`}, {1.8913812543459405`,
     2.319493673670671`}, {1.7674367495794008`, 
    3.5325088938533917`}, {1.9320143413972677`, 
    3.611151717969884`}, {1.8107126970778207`, 
    3.6790935311559876`}, {1.804725058474518`, 
    3.8520729911980185`}, {1.805368854133445`, 
    3.852728507949241`}, {1.563939776616803`, 3.9192958794919646`}};


VoronoiMesh[pts, domspec, Epilog -> {Point[pts], Red, Circle[pts[[6]], 0.1]}, 
 PlotRange -> domspec, ImageSize -> 500, PlotRangePadding -> Scaled[.05]]
(* snip *)

Clear[allpts]
allpts = Join @@ Join @@ Table[p + dom {i, j}, {i, -1, 1}, {j, -1, 1}, {p, pts}];

VoronoiMesh[allpts, (* domspec, *) Epilog -> {Point[allpts], Red, Circle[pts[[6]], 0.1], Circle[allpts[[3 Length[pts] + 6]], 0.1]}, 
PlotRange -> domspec, ImageSize -> 500, PlotRangePadding -> Scaled[.05]]

giving

The new dom is smaller than the original one.

Oftenly tools transform input data into a natural standard area, do their magic and transform the results back into users area, you might want to check here, whether tetgen does so.

It will not work with

dom  = N[1.25 Pi, 30]

At least VoronoiMesh or tetgen or both of them suffer from a spurious domain dependency.
.

[sghuisman (Sander Huisman)]

Hi Udo,

I indeed found the same. The results can be domain-dependent, which it shouldn’t be, unless point fall outside the domain of course. The more worrying is that the number of mesh-polygons returned is not equal to the number of points! I would not care so much if the angle of the polygon was slightly wrong because of rescaling, it is just very strange the number of points change! It should at least show a warning message!

Rescaling seem also strange to me though; the distance is still relatively large compared to the common double precision, and even single precision…

well… I send it to Wolfram bug-support, let’s see what they see, hope they fix it, not so much for me, more for other people who got frustrated!

[sukrause (Udo Krause)]

Call this

Clear[pts]
pts = Flatten[Table[{x, Re[BesselJ[\[Pi] o I, x + o I]]}, {o, 1, 8}, {x, 0, 5, 0.1}], 1];
VoronoiMesh[pts, Epilog -> Point[pts], PlotRange -> {{-1/2, 1}, {-1, 1}}]

to see some cells with more than one point inside:

These points are neither microscopically near to each other nor in some special position, it’s plainly wrong.
Someone would like to file a bug.

[sghuisman (Sander Huisman)]

O NO!!! That’s BAD! How did you find this case? just trial and error? (trial and terror!)

I’ll submit it, so they get inundated by me about the same issue, and will have a good look at it! Thanks a lot!

[sukrause (Udo Krause)]

Oh yes, this is so not fair (you spend a day to prepare an example and now it’s obvious … sorry)!

I tried to poke some fun at VoronoiMesh and was pretty sure, that it will give nice tortoise-like cells - but this result is disconcerting.

[sghuisman (Sander Huisman)]

Well… it does illustrate that it is not just 1 weird case, apparently it shows up in many more cases…
VoronoiMesh should at least spawn a error message that something went wrong!

[henrikschachner]

Hi all,

regarding the example from Udo I guess the problem arises from the wide range where the points lay in:

Clear[pts]
pts = Flatten[Table[{x, Re[BesselJ[\[Pi] o I, x + o I]]}, {o, 1, 8}, {x, 0, 5, 0.1}], 1];
MinMax /@ Transpose[pts]
(*  {{0.`,5.`},{-23917.819047547382`,33188.76634356813`}}  *)

As one can see the y-range is just huge. Playing around with a limitation of that y-range a gradual occurrence of that problem can be observed:

pts1000 = Select[pts, Abs[Last[#]] < 1000 &];
pts2000 = Select[pts, Abs[Last[#]] < 2000 &];
pts3000 = Select[pts, Abs[Last[#]] < 3000 &];
pts5000 = Select[pts, Abs[Last[#]] < 5000 &];
VoronoiMesh[#, Epilog -> Point[pts], PlotRange -> {{-1/2, 1}, {-1, 1}}, ImageSize -> 250] & /@ {pts1000, pts2000, pts3000, ts5000}

giving:

The effect also could be provoced by the sheer number of points. In any case this cannot be the reason for the original problem.

Regards – Henrik

[sghuisman (Sander Huisman)]

But even with single-precision numbers it would be more than enough to present these numbers with great accuracy! I really don’t understand what happens here, I guess it is some parsing error from and to tetgen, I assume tetgen is ‘bullet-proof’ enough…

I submitted the issue (twice) to the wolfram people, so they can fix it. Once I know something concrete, I will post it here.

Thanks for further looking in to it!

[sukrause (Udo Krause)]

Giving a domspec without selecting from the point set does not work either

Clear[pts]
pts = Flatten[Table[{x, Re[BesselJ[\[Pi] o I, x + o I]]}, {o, 1, 8}, {x, 0, 5, 0.1}], 1];
VoronoiMesh[pts, {{-1/2, 1}, {-1, 1}}, Epilog -> Point[pts], PlotRange -> {{-1/2, 1}, {-1, 1}}]

selecting from the point set works with and without domspec, but must not necessarily give the same set of cells the whole point set gives:

Clear[pts, domspec]
pts = Flatten[Table[{x, Re[BesselJ[\[Pi] o I, x + o I]]}, {o, 1, 8}, {x, 0, 5, 0.1}], 1];
domspec = {{-1/2, 1}, {-1, 1}};
VoronoiMesh[#, domspec, Epilog -> Point[#], PlotRange -> domspec]&[Select[pts, Norm[#] <= Norm[domspec] &]]

To select the right way one had roughly spoken to take into account all pairs of points whose perpendicular bisector has a non-empty intersection with the domspec area, a thing VoronoiMesh should do on it’s own if it is given a domspec.

[sukrause (Udo Krause)]

Just for the logs and the bug tracker: in this case \lambda=4, \mu=5 are minimal to invoke the error

Clear[pts]
With[{\[Lambda] = 4, \[Mu] = 5},
 pts = Flatten[Table[{x, Im[BesselJ[\[Pi] o I, x + o I]]}, {o, 0, \[Lambda]}, {x, 0, \[Mu], 0.1}], 1];
 VoronoiMesh[pts, Epilog -> Point[pts], PlotRange -> {{-1/10, 5/6}, {-1/2, 1/2}}]
]

which imposes now on VoronoiMesh[] that way

Note

In[122]:= MinMax /@ Transpose[pts]
Out[122]= {{0., 5.}, {-853.437, 3364.88}}

This is \frac{33188.819 - (-23917.866)}{3364.88 - (-853.437)}=13.5378 more than one order of magnitude smaller as in the previous example and still a failure.

[sghuisman (Sander Huisman)]

Indeed very interesting, thanks for further exploring. I also did a bit myself, and went to the simplest case of 3 points (on a diagonal):

lower=1.0;
upper=10^9;
n=40;
offset=1000;
m=5;
While[n-->0,
    test=(upper+lower)/2;
    (*Echo[{upper,lower,test},n];*)
    pts={offset+{0.0,0.0},offset+m{1.0,1.0},{test,test}};
    vm=VoronoiMesh[pts];
    cells=MeshCellCount[vm,2];
    If[cells==2,
            upper=test
        ,
            lower=test
    ]
];

(test-offset)/(m Sqrt[2])

if you set the offset to 0 or to 1000 and you set m to 1 or to 10, always the ratio between the maximum diagonal and the distance between the nearest two points is 9860.55 which is:

\frac{10^6}{100+\sqrt{2}}

probably 10^4 with some padding…

[sukrause (Udo Krause)]

Interesting, separation on a diagonal can also be a problem! Modifying your code into

lower = 1.0; upper = 10^9; n = 40; offset =(* 1000; *) 0;
m = 1; 
While[n-- > 0, test = (upper + lower)/2; 
  pts = {offset + {0.0, 0.0}, offset + m {1.0, 1.0}, {test, test}}; 
  vm = VoronoiMesh[pts]; 
  cells = MeshCellCount[vm, 2];
  (* Print["ult: ", {upper,lower,test},"| n: ",n, "| cells = ",cells]; *)
  If[cells == 2,
   Print["pts: ", pts];
   upper = test,(* else *)
   lower = test]
 ];

then

m = 0.1:   pts: {{0., 0.}, {0.1,0.1}, {1394.5, 1394.5}}
m = 1:     pts: {{0., 0.}, {1., 1.},  {13944.9, 13944.9}}
m = 5:     pts: {{0., 0.}, {5., 5.},  {69724.6, 69724.6}}

all are defective and the defect scales linear

In[22]:= 50 1394.5
Out[22]= 69725.

Well, with a trustworthy MeshCellCount the error is straightforward to detect, given it is suspected - which it is now, after your original observation.

[sghuisman (Sander Huisman)]

I forwarded the error a few days ago, I told them they should do this check nonetheless, if the number of points != MeshCellCount then it should spawn an error message. This could happen due to wrong rounding, but also because one gives 2 identical points, for which the Voronoi diagram is ill-defined…