Falls on 10000 random points on Win32 on stack overflow:

[GoogleCodeExporter]

What steps will reproduce the problem?
1. Under windows 32 with just 10000 random points you get stack overflow

Please provide any additional information below.

Here is my code used:

// PolyDelaunay.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "poly2tri.h"

#include "stdafx.h"
#include "float.h"
#include "time.h"

#include <vector>
#include <cstdlib>
#include <iostream>
#include <algorithm>

using namespace std;

bool pointSortPredicate(const p2t::Point& a, const p2t::Point& b)
{  
  if (a.x < b.x)
    return true;
  else if (a.x > b.x)
    return false;
  else if (a.y < b.y)
    return true;
  else
    return false;  
};

bool pointComparisonPredicate(const p2t::Point& a, const p2t::Point& b)
{
  return a.x == b.x && a.y == b.y;
}

int _tmain(int argc, _TCHAR* argv[])
{
  size_t numPoints = 10000;
  std::vector<p2t::Point> points;

  for (size_t i = 0; i < numPoints; ++i)
  {
    p2t::Point p;
    p.set(rand(), rand());
    points.push_back(p);
  }

  std::cout << "Number of points: " << numPoints << std::endl;

  std::sort(points.begin(), points.end(), pointSortPredicate);
  std::vector<p2t::Point>::iterator newEnd = std::unique(points.begin(), points.end(), pointComparisonPredicate);
  points.resize(newEnd - points.begin());
  numPoints = points.size();

  std::cout << "Points after unique making: " << numPoints << std::endl;

  /*
  size_t numPoints2 = numPoints * 3;  

  double offset = 1.0e-5;
  for (size_t i = 0; i < numPoints; ++i)
  {
    GEO_Vrtx2D p = points.at(i);
    p.Set(p.GetX() + offset, p.GetY() + offset);
    points.push_back(p);
    p.Set(p.GetX() + offset, p.GetY() + offset);
    points.push_back(p);
  }
  numPoints = numPoints2;
  */

  // Determine min, max
  double minx = DBL_MAX;
  double miny = DBL_MAX;
  double maxx = -DBL_MAX;
  double maxy = -DBL_MAX;

  for (size_t i = 0; i < numPoints; ++i)
  { 
    if (points.at(i).x  < minx)
      minx = points.at(i).x;
    if (points.at(i).x > maxx)
      maxx = points.at(i).x;
    if (points.at(i).y < miny)
      miny = points.at(i).y;
    if (points.at(i).y > maxy)
      maxy = points.at(i).y;
  }

  // Outer rectangle 
  std::vector<p2t::Point*> polyline;
  polyline.push_back(new p2t::Point(minx, miny));
  polyline.push_back(new p2t::Point(maxx, miny));
  polyline.push_back(new p2t::Point(maxx, maxy));
  polyline.push_back(new p2t::Point(minx, maxy));

  clock_t start, end, diff;
  start = clock();

  p2t::CDT* delaunay = new p2t::CDT(polyline);

  // Add points
  for (size_t i = 0; i < numPoints; ++i)
  {
    delaunay->AddPoint(new p2t::Point(points.at(i).x, points.at(i).y));
  }

  delaunay->Triangulate();

  end = clock();

  diff = ((end - start) * 1000) / CLOCKS_PER_SEC;
  std::cout << "Run time: " << diff << std::endl;

  int a;
  cin >> a;

    return 0;
}

Original issue reported on code.google.com by pavel.ce...@celba.cz on 16 Nov 2011 at 2:10

[GoogleCodeExporter]

Falls in SweepContext::MeshClean() calling itself a lot of times.

If this gets fixed soon, I may test this more massively. Up to 1M of points for 
performance and also with data generating floating-point round-off errors like 
commented in code above.

Original comment by pavel.ce...@celba.cz on 16 Nov 2011 at 2:25

[GoogleCodeExporter]

just took a quick glance at the code.

When adding steiner points they must be inside the rectangle polygon.
Since you use the point set to find max an min x,y values there should be a 
chance that one of the points in the set will be on the edge of the rectangle 
you create, maybe this is your problem.

The triangulation algorithm assume that given polygon isn't self intersecting 
or touching and that steiner point are inside the polygon and not on the edge 
or outside.

Original comment by thahlen@gmail.com on 16 Nov 2011 at 9:12

[GoogleCodeExporter]

Well.. not only a chance, atleast 2 points would be on the corner or edge of 
the rectangle. Try growing the rectangle a tiny bit and see if it works better.

Original comment by thahlen@gmail.com on 16 Nov 2011 at 9:25

[GoogleCodeExporter]

Did this suggestion solve your problem? Want to know if I can close this issue 
:)

Original comment by thahlen@gmail.com on 19 Nov 2011 at 7:43

[GoogleCodeExporter]

Adding these lines after min-max finding:
  minx -= 100;
  miny -= 100;
  maxx += 100;
  maxy += 100;

Didn't solve the problem - I note that values generated are from 0 to approx. 
32K so adding 100 to boundary should be more than enough.

Original comment by pavel.ce...@celba.cz on 21 Nov 2011 at 8:23

[GoogleCodeExporter]

Ah I think I know the issue. 

The MeshClean algorithm is recursive. It traverse the triangle graph 
recursively in a polygon until it finds an edge. Since we only got a rectangle 
with 4 edges and a massive amount of internal triangles the recursive depth can 
get really big.

I know someone rewrote MeshClean to be nonrecursive: see bottom of this issue

http://code.google.com/p/poly2tri/issues/detail?id=12&can=1

this is not the c++ version but I guess you could just port it. I agree that 
current MeshClean is a weakness and probably should be rewritten.

Original comment by thahlen@gmail.com on 21 Nov 2011 at 1:36

[GoogleCodeExporter]

Well I think you should really fix that in c++ version. I have switched to use 
q-hull library for delaunay triangulation as even when you fix the issue, I'm 
not quite sure that your algorithm is really stable.

Original comment by pavel.ce...@celba.cz on 21 Nov 2011 at 4:36

[GoogleCodeExporter]

Well this lib is excellent for triangulating polygons with valid input. eg. no 
duplicate points, touching or intersecting edges. For that I'd say it is really 
stable. Thats what it was written for originally.

qhull is a nice lib, so unless you need constrained triangulation it is 
probably a better solution.

Original comment by thahlen@gmail.com on 21 Nov 2011 at 5:27

[GoogleCodeExporter]

Well yes, the qhull is really working, but the interface is terrible. Our 
company would however also use stable constrained triangulation (surface 
remeshing, etc...). But my point is that anyone who downloads the code will 
test it like I did and stack overflow errors on first run isn't probably the 
image of library you would like to have.

Original comment by pavel.ce...@celba.cz on 22 Nov 2011 at 12:18

[GoogleCodeExporter]

The c++ version is a simplified port of the java version I wrote. The java 
version separates polygon and pointset triangulations and have different 
MeshClean functions so this problem never arise. 

I have notified mason and hopefully he can implement a non recursive mesh clean 
for the c++ version. I think most have used this lib as a polygon tesselator 
and not used that many steiner points.

Original comment by thahlen@gmail.com on 22 Nov 2011 at 4:46

[GoogleCodeExporter]

I can fix this, although I need to make the time. This is a hobby project for 
me, not used professionally, and I have too many other side interests I'm 
pursing at the moment. 

I have a week off from work over the holidays, and will commit several hours to 
finding a good solution.

Original comment by mason.gr...@gmail.com on 14 Dec 2011 at 2:56

[GoogleCodeExporter]

With 10000 random points, you may have some duplicate points. You need to 
provide poly2tri with a valid input set to work correctly, and all the points 
should be contained within your constrained edges.

I will add a non-recursive version of MeshClean for good measure.

Original comment by mason.gr...@gmail.com on 4 Feb 2012 at 8:41

• Changed state: Invalid

[GoogleCodeExporter]

Well my code is obviously removing the duplicated points so that's not the 
problem.

Original comment by pavel.ce...@celba.cz on 5 Feb 2012 at 5:07

[GoogleCodeExporter]

I have the same problem: With around 20000 Steiner points (the limit varies), I 
get a stack overflow in SweepContext::MeshClean. I make sure that there are no 
duplicate points.

Could you please provide the promised non-recursive version of MeshClean? Would 
help me a lot.

Original comment by martin.h...@gmail.com on 25 Jun 2012 at 11:26

[GoogleCodeExporter]

Well I worked out a replacement for MeshClean, using a std::stack instead of 
the "real" stack. I tested it with 300,000 Steiner points, it worked. For 
10,000 Steiner points, the destruction of the CDT object takes exactly the same 
time with your and my version (both 750 msec). Your testbed also runs without 
errors.

I wonder if this could be implemented more efficiently - destruction of a CDT 
object with 300,000 Steiner points takes as long as 38 sec (but at least there 
is no overflow :) ).

Attached you'll find the code.
* added  #include <stack>
* new version of SweepContext::MeshClean
* (no other changes)

Original comment by martin.h...@gmail.com on 25 Jun 2012 at 12:40

Attachments:

• sweep_context.cc

[GoogleCodeExporter]

Well instead of calling destructor allocate on custom memory managed heap and 
throw the memory in heap away (without calling destructors).

Original comment by pavel.ce...@celba.cz on 1 Jul 2012 at 5:47