Message Boards Message Boards

3
|
8522 Views
|
18 Replies
|
17 Total Likes
View groups...
Share
Share this post:

VoronoiMesh: user error, tetgen error, else?

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]]

enter image description here

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]]

enter image description here

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.
POSTED BY: Sander Huisman
18 Replies
Sort By:

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...
POSTED BY: Sander Huisman

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.
POSTED BY: Udo Krause

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...
POSTED BY: Sander Huisman

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

enter image description here

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.
POSTED BY: 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.

POSTED BY: Udo Krause

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!
POSTED BY: Sander Huisman

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:

enter image description here

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
POSTED BY: Henrik Schachner

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!
POSTED BY: Sander Huisman

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.

POSTED BY: Udo Krause

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!
POSTED BY: Sander Huisman

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:

enter image description here

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.
POSTED BY: Udo Krause

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!
POSTED BY: Sander Huisman

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

enter image description here

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. .
POSTED BY: Udo Krause
POSTED BY: Sander Huisman

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

enter image description here
POSTED BY: Udo Krause

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!
POSTED BY: Sander Huisman

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
POSTED BY: Henrik Schachner
POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract