Problem to build diagram on dense point set

[Mr-Question]

Hello Pierre,

I continue to experiment with your algorithm. Now I'm trying to run it on some real data.

I have attached file with points. When I read it, I have got segmentation fault during intersection stage. Seems it is caused by duplicating points. The matter is that halfEdge in do-while loop in intersect(Box box) method becomes a null pointer. If I change condition to "while (halfEdge && halfEdge != site.face->outerComponent);" it goes to infinite loop until it hangs entire system.

If duplicates are removed, it works, but produces redundant (to my view) link crossing several zones in the middle:

Best Regards! case.txt

[pvigier]

Hi,

Indeed duplicated points are not supported by the algorithm. The reason is that it makes little sense: the two cells would be exactly the same and the Voronoi diagram would not be a partitioning of the plane anymore which is not desirable.

I will add in the documentation that the points must all be unique. Maybe I will add a check at the beginning of the algorithm to ensure that.

After removing the duplicated points, I obtain the same diagram as you. I am almost sure the problem does not come from the density of points. But rather that some points are aligned along the y-axis.

For instance, if I only take your first five points (without the duplicates):

1 0
0 0
1 1
0 1
0.795082 0.385057

I obtain the same problem (a vertical line and the intersection algorithm makes a seg fault ...):

However, if I slightly change the y-coordinates so that no two points are aligned along the y-axis anymore:

1 0.0000001
0 0
1 1.0000001
0 1
0.795082 0.385057

Then, the diagram is correct:

Slightly changing the y-coordinates of (1, 0) and (1, 1) also works on your full example:

This problem is a known issue and I will work on it.

Thanks for reporting!

[Mr-Question]

Hello,

just added one more use case with points with same Y coordinates. Hope this will help to check the solution in further.

circle.txt

[Mr-Question]

Hello Pierre!

Just wanted to ask, do you plan to continue development around your project? It would be great to see it as a common and robust solution.

Best wishes in your continued success!

[pvigier]

Hello!

I was working on another project for a while but I will get back to this one shortly. I hope in the next month.

Thanks for your interest!

[pvigier]

Hi!

I have updated the project.

I fixed lots of edge cases. The algorithms are way more robust now. In particular, I have added some tests with some special configurations (line: all points aligned, grid: points are placed on a 2D grid, circle: points are on a circle) where many points share the same x or y coordinates and many cells meet on a same vertex. And it works fine.

It is still possible to find configurations which lead to error but it is harder. I don't know if it is really possible to prevent all errors with the method I used to manage floating point issues. But I want to keep the code as simple as possible.

The next steps are to set up the CI and improve the cmake file so that it is easy to install and use the library.

Best!

[Mr-Question]

Hello!

Thanks a lot, it has become much better. However, I have faced some use cases that prevent integration into target project. I could share some of them later if you are interested.

Best regards!

[pvigier]

Hello!

That would be very nice if you find some time! If it is not too complicated to fix, I will gladly fix these issues.

Best

[Mr-Question]

Hi, Pierre!

Thanks for your efforts, the component works quite fast on huge and dense point clouds, and I continue its evaluation.

I am not really sure whether I should post it in this tread or create a new one. The matter is that there could be holes in triangulation produced by MyGAL. Maybe I just use this functionality quite improperly.

Here are the code how I put the source points:

std::vector<mygal::Vector2<double>> points;
points.push_back(Vector2<double>(20, 9.927088740980544));
points.push_back(Vector2<double>(19.70941817426052, 12.32024538385613));
points.push_back(Vector2<double>(18.8545602565321, 14.57432046141823));
points.push_back(Vector2<double>(17.48510748171101, 16.5583153233885));
points.push_back(Vector2<double>(15.68064746731156, 18.15692739991711));
points.push_back(Vector2<double>(13.54604887042536, 19.2772511678347));
points.push_back(Vector2<double>(11.20536680255323, 19.85417748196109));
points.push_back(Vector2<double>(8.794633197446771, 19.85417748196109));
points.push_back(Vector2<double>(6.453951129574645, 19.2772511678347));
points.push_back(Vector2<double>(4.319352532688443, 18.15692739991711));
points.push_back(Vector2<double>(2.514892518288992, 16.5583153233885));
points.push_back(Vector2<double>(1.145439743467904, 14.57432046141824));
points.push_back(Vector2<double>(0.2905818257394817, 12.32024538385613));
points.push_back(Vector2<double>(0, 9.927088740980556));
points.push_back(Vector2<double>(0.2905818257394763, 7.53393209810498));
points.push_back(Vector2<double>(1.145439743467895, 5.279857020542872));
points.push_back(Vector2<double>(2.514892518288978, 3.295862158572605));
points.push_back(Vector2<double>(4.319352532688426, 1.69725008204399));
points.push_back(Vector2<double>(6.453951129574625, 0.5769263141264035));
points.push_back(Vector2<double>(8.794633197446746, 7.105427357601003e-15));
points.push_back(Vector2<double>(11.20536680255321, 0));
points.push_back(Vector2<double>(13.54604887042533, 0.5769263141263857));
points.push_back(Vector2<double>(15.68064746731153, 1.697250082043964));
points.push_back(Vector2<double>(17.48510748171099, 3.29586215857257));
points.push_back(Vector2<double>(18.85456025653208, 5.279857020542829));
points.push_back(Vector2<double>(19.70941817426051, 7.533932098104933));

Then I invoke MyGAL:

mygal::FortuneAlgorithm<double> aTool (points);
aTool.construct ();
mygal::Diagram<double> aDiagram = aTool.getDiagram ();
mygal::Triangulation aTriangulation = aDiagram.computeTriangulation ();

And here is the code that I use to produce OBJ file with triangulation:

for (size_t aNodeIt = 0; aNodeIt < points.size (); ++aNodeIt)
{
  std::cout.precision (16);
  std::cout << "v " << points[aNodeIt].x << " " << points[aNodeIt].y << " 0" << std::endl;
}

for (size_t i = 0; i < aDiagram.getNbSites (); ++i)
{
  int aNodes[3] = { -1, -1, -1 };
  const std::vector<mygal::Diagram<double>::Site>& aSites = aDiagram.getSites();
  aNodes[0] = static_cast<int> (aSites[i].index + 1);

  const std::vector<std::size_t>& aNeighbors = aTriangulation.getNeighbors (i);
  for (size_t j = 0; j < aNeighbors.size (); ++j)
  {
    const int aDstIdx = static_cast<int> (aSites[aNeighbors[j]].index + 1);

    if (aNodes[1] == -1)
    {
      aNodes[1] = aDstIdx;
    }
    else
    {
      aNodes[2] = aDstIdx;

      std::cout << "f " << aNodes[0] << " " << aNodes[1] << " " << aNodes[2] << std::endl;

      aNodes[1] = -1;
    }
  }
}

Finally, I have got the following result:

result.zip

Thank you in advance!

[Mr-Question]

Dear Pierre!

I have just came to that the code intended to collection of triangles is messy a bit... In fact, there is a part, that I intentionally excluded, which checks that the same triangle is not faced more than once.

Here is the question: how it would be possible to collect unique triangles in the most appropriate and fastest way?

What I actually would like to ask in addition: as the algorithm uses half-edge data structure, is it possible to use this feature to speed up collection of triangles (without maps) by introduction of some kind of iterator over half-edges and extending half-edge structure (or iterator itself) by bool flags that indicate visited ones?

Finally, I have found one more case failing production of correct diagram regular.txt :

[pvigier]

Hi,

Indeed, I think that in your first example (with the circle), the problem comes from the way you enumerate the triangles: it is not exhaustive. The Triangulation class is really simple and just returns the neighbors of sites in Delaunay triangulation.

If you want the triangles, I think that the best way is to use the vertices of the diagram. Each vertex is the center of the circumcircle of three sites. So for each vertex of the Voronoi diagram, there should be a triangle in the Delaunay triangulation.

The next question is how to get the three sites for each vertex. Once you have a half-edge that is incident with a vertex you can obtain the other half-edges and thus the other sites using twin, prev and next. The problem is that with the current implementation the Vertex class is really simple and just contains the coordinates of a point, there is no way to get an incident half-edge. We have to iterate the half-edges but we will process each vertex several times (at most three times, I think) which is not really efficient. The best way would be to keep track of one incident half-edge for each vertex (as it is already done for faces which keep track of one half-edge on the outer component). It should be possible to do this with no significant overhead.

I will try to find some time to implement this.

Concerning, the second issue, I am sorry but I do not think I can do anything more. I have already done as much as I can to prevent floating point issues from occurring. But with the method I used, I can not prevent them all. In this reddit comment, someone advises me to use Shewchuk's robust floating point primitives. But I think it is beyond my abilities to implement that.

Moreover, I am happy enough with the robustness of the current implementation as I only use it with random points, so there is really low probability to have aligned points. If you have data that contain a lot of aligned points, you should maybe consider another implementation that uses more advanced techniques, as the one mentioned above, to manage floating point issues.

Best, Pierre

[Mr-Question]

Well, in any case your implementation is a very good cross-platform solution. Hope, it will receive new features from time to time . Thank you for your time, Pierre.

Regards!