Angus Johnson wrote:

I believe Clipper is currently substantially more reliable that GPC because I don't think GPC properly deals with complex intersections. I'm not aware of the tree data structure recommended by Vatti et al. (I couldn't find any mention of it in Vatti's original paper.) My implementation uses linked and double linked list data structures which are described in Vatti's paper. I would be very interested if you could provide a link to Vatti's and Bentley-Ottman's discussion of a tree structured approach.

Concerning robustness - I believe numerically robust problems only arise when complex intersections are encountered. Specifically, how does an algorithm cope with mulitple intersections occuring in very close proximity. Simply relying on 'tolerance' when using floating point values will guarantee failure at some point. On failure, intersection operations are performed in an incorrect order resulting in incorrect edge sorts. This fortunately is easily detected after each sweep. While, not a final solution, I have an update that's close to release that tests the edge sorts before actually processing the intersections. If the test sort order is wrong, the 'tolerance' value is adjusted (both up and down) and retested. Fortunately, this testing has only a minor impact on performance (%10). In my stress tests*, and using this 'tuned' tolerance approach, the incedence of test failure drops to about one in a million tests. (*The test involves generating random subject and clip polygons with 100 vertices each rotated at varying degrees where the vertices are specifically tuned to generate may complex intersections.) Given that complex intersections in clipping operations are infrequent events in normal use (which is why GPC rarely blinks in normal use), I would argue that the likelihood of users encountering a failure with Clipper would be very much lower than one in a million. You also mentioned scaling as an issue with floating point implementations. I agree that it is an issue and that input polygons with very small or very large coordinate values should be scaled. Anyhow, a fully numerically robust solution is still my final goal.

That's a worthwhile goal to be sure. When I discussed scalability I meant algorithmic complexity and ability to process large data volumes in input geometry in a single test case. An algorithm that scales near linearly will be able to process inputs on the order of one hundred million verticies in reasonable run time, for example. If your sweepline data structure is a linked list and not a tree or bucketed list then your algorithmic complexity must be worse than the O(n log n) wrt input verticies plus intersection vertices described in Vatti's paper. In that case I wouldn't expect it to scale past tens of millions of vertices in reasonable runtime. In my testing the empirical scaling factor for GPC was about O(n^1.75) so you could be substantially faster than that but still worse than O(n log n). While you mention that odds of failure may be one in a million tests for tests of a given size, if you increase the size of the input by a factor of one million your odds of failure get close to one in one. That was the point I was trying to make about algorithmic scalability being undermined by lack of robustness. It doesn't matter if your implementation could process one hundred million vertex test cases in reasonable runtime if the odds of it succeeding are remote. I just wanted to clarify this point because it isn't always intuitively obvious. Also, I believe GPC does handle high order intersections because I'm quite certain I was generating large numbers of them in the tests I was running that GPC seemed to handle just fine. I could be wrong of course, but I suggest you double check that GPC doesn't properly deal with complex intersection. I know it doesn't deal with overlapping or self intersecting polygons at all, but that is a separate issue from polygons that have vertices in common or polygons that abut along edges but do not overlap. GPC is used quite extensively and paid for by many people, we should be responsible about making claims about its short comings. Regards, Luke

Angus Johnson wrote:

...

That's a worthwhile goal to be sure. When I discussed scalability I meant algorithmic complexity and ability to process large data volumes in input geometry in a single test case. Ahh, thanks for the clarification. I haven't yet done any significant testing of processing very large data volumes, nor am I competent to determine an O value, so I'll defer to your estimation there.

What you have already done is significantly harder than determining O values. Don't sell yourself short. If you plot number of vertices in input plus output on the x axis and runtime on the y axis in a excel spreadsheet you can have excel create an exponential trend line for the data and display the equation on the graph. The equation will look like 2.79235 * x^1.34984. Your empirical scaling factor (empirical O value) is then 1.35. In this hypothetical example 1.35 is faster than O(n * sqrt n) but slower than O(n log n). You should be able to estimate your empirical scaling factor very quickly by generating progressively large random test cases with this method. Read the experiment section of my boostcon paper http://www.boost.org/doc/libs/1_44_0/libs/polygon/doc/GTL_boostcon2009.pdf to get an idea how to structure the experiment and what the graph should look like (page 9.) I can share the code for the test harness with you if you like. The tricky part is making sure that the geometry doesn't qualitatively change as you make successively larger test cases because this can introduce a systematic bias in you empirical scaling factor. I'd be interested in seeing this kind of performance analysis on your clipper implementation. If you are using a linked list as your sweepline data structure your worst case runtime complexity is O(n^2) wrt. input plus intersection vertices and your expected runtime complexity is O(n^1.5) in the common case. In the worst case there is O(n) line segments in your sweepline at each sweepline "stop" (unique x coordinate) while in the common case there is O(n^0.5) line segments in your sweepline at each sweepline stop and O(n) stops.

...

Also, I believe GPC does handle high order intersections because I'm quite certain I was generating large numbers of them in the tests I was running that GPC seemed to handle just fine. I could be wrong of course, but I suggest you double check that GPC doesn't properly deal with complex intersection. I know it doesn't deal with overlapping or self intersecting polygons at all, but that is a separate issue from polygons that have vertices in common or polygons that abut along edges but do not overlap. GPC is used quite extensively and paid for by many people, we should be responsible about making claims about its short comings.

I was careful to qualify my statements about GPC with "I don't think" and "I believe". However, I accept that this still isn't fair treatment and I should be sure before making these claims. Let me simply state that under my stress-testing, GPC fails about 4% of the time (and with extreme stress testing - rounding random vertices to nearest 100 instead of nearest 10 - it fails far more frequently). The full source code for my stress-tests can be downloaded here: http://www.angusj.com/delphi/clipper.php

What is GPC's failure mode that you observe when it fails your stress test? Does it crash or does it just produce garbage output polygons. In my testing I found that given self intersecting or overlapping polygons in an input GPC produced garbage output. I would prefer a crash since it is as much work to check the result of a boolean is correct as to perform the boolean. An algorithm that produces the wrong result and reports success is worse than an algorithm that crashes. If you are substantially more robust than GPC they you are substantially robust. Epsilon methods won't get you to 100%, but there are other methods that you can use that get you all the way. You are asymptotically approaching them with your iterative epsilon idea. The really tough problem wrt. robustness is that numerically error move line segments laterally and can introduce spurious intersections in the output. These intersections become fatal errors for down stream processing of the output of your algorithm. There is a paper on floating point snap rounding by Milenkovic that could help you with 100% numerical robustness that avoids spurious intersections. Search for Polygon Set Operation and or Milenkovic and use cite seer to walk the reference links. I can send you my archive of papers if you don't have access to an academic database. I greatly approve of your stress testing approach. Try turning off rounding of the random vertices and then we'll really see where epsilon methods fall short. It isn't really fair to round the inputs with your knowledge of which epsilon you use and compare to someone elses algorithm. The kind stress testing you have been doing is exactly the right way to stress test numerical robustness. I wish others would do the same. Perhaps you can include the Boost.Geometry clipping implementation in your scaling and stress testing. If your implementation is both faster and more reliable than what they have I would think we would use yours to replace theirs in their library. I don't think users of the library would find that confusing. Also boost geometry is something of a test bed for development of better algorithms in addition to a vehicle for releasing them. From that standpoint multiple algorithms that do the same thing are part in parcel of developing better algorithms. Regards, Luke

Yesterday afternoon I did just that and I get a trendline with an equation of y = 3E-08x2.8413 which suggests scaling will be a problem.

Can't sleep :(. Just realised I was using the wrong values in my chart. The number of edges was actually the square of the values I put into the chart. I'd designed my test polygon by creating a square filled with horizontal squiggly lines designed to maximally stress the vatti algorithm (since it results in a near maximal number of LocalMinima with left and right bounds all of length 1). Anyhow, reworking the data with the correct number of edges, the equation for this (near worst case) test polygon returns an equation of y = 1E-08x^1.4821 . If I generate the squiggly lines vertically in my test polygon (and destress the vatti algorithm) which is likely to be closer to normal use, I suspect this equation will drop further. Will leave that till the morning.

Angus Johnson wrote:

...

Anyhow, reworking the data with the correct number of edges, the equation for this (near worst case) test polygon returns an equation of y = 1E-08x^1.4821 . If I generate the squiggly lines vertically in my test polygon (and destress the vatti algorithm) which is likely to be closer to normal use, I suspect this equation will drop further.

OK, when I rotated the test polygon 90 deg the trendline equation becomes : y = 2E-07x^1.0841 I would suggest that in normal use, my clipping library would perform somewhere between O(n^1.0841) and O(n^1.4821) though I suspect nearer to the lower value.

Here's the function that produces my test polygons (generating what I believe would be near worst case vertices for the Vatti algorithm):

//struct TDoublePoint { double X; double Y; }; //typedef std::vector< TDoublePoint > TPolygon;

TPolygon MakeTestPoly(int startX, int startY, int boxSize){ int k,x,y;

TPolygon result; result.resize(boxSize * boxSize +2); x = startX+4; y = startY; k = 0; for (int i = 1; i <= boxSize; ++i) { for (int j = 1; j <= boxSize; ++j) { x = x +(i % 2)*4 -2; y = y + (j % 2)*2 -1; result[k].X = x; result[k].Y = y; k++; } y = y + 4; x = x +(i % 2)*4 -2; } result[k].X = startX; result[k].Y = y -4; result[k+1].X = startX; result[k+1].Y = startY; return result; }

Just a C++ style pointer, we have what we call resource acquisition is initialization (RAII) policy where we initialize all variables at the time they are declared rather than deferring that to later. The rationale for this is that it avoids bugs introduced later when a variable is attempted to be used between where it was declared and where it is initialized resulting in undefined behavior and different behavior on different platforms/compiler optimization flags. RAII becomes much more important when memory allocation, locks, file handles, sockets and other tricky to handle resources are thrown into the mix because it allows you to encapsulate the acquisition and reliquishment of those resources in the declaration and destructor call for your objects which is much less error prone than the alternative. TPolygon MakeTestPoly(int startX, int startY, int boxSize){ int k = 0; int x = startX+4 int y = startY; TPolygon result; //it would be better if the TPolygon constructor accepted the size result.resize(boxSize * boxSize +2); //making this line uneeded //and to help prevent people from using the TPolygon but forgetting to resize it first. for (int i = 1; i <= boxSize; ++i) { for (int j = 1; j <= boxSize; ++j) { x = x +(i % 2)*4 -2; y = y + (j % 2)*2 -1; //squiggles right twices as fast as it squiggles up result[k].X = x; result[k].Y = y; k++; } y = y + 4; x = x +(i % 2)*4 -2; } result[k].X = startX; result[k].Y = y -4; result[k+1].X = startX; result[k+1].Y = startY; return result; } This case is actually a common case and not a worst case for your algoirthm because the average number of line segments intersecting your sweepline is more like O(n^1.5) because the line segments are very short on average. If you instead make n horizontal rectangles with random start and end coordinates and a small height within separate ranges of X for the start and end you will get O(n) unique x coordinates and on average O(n) (roughly 1/2 n) line segments intersecting your sweepline data structure at each unique x. This will be the worst case input and have expected runtime complexity of O(n^2) using a linked list as the sweepline data structure while an optimal Vatti using a tree would have O(n log n) complexity and empirical scale factor of about O(n^1.05). Regards, Luke

Angus Johnson wrote:

...

Just a C++ style pointer, we have what we call resource acquisition is initialization (RAII) policy where we initialize all variables at the time they are declared rather than deferring that to later. The rationale for this is that it avoids bugs introduced later when a variable is attempted to be used between where it was declared and where it is initialized resulting in undefined behavior and different behavior on different platforms/compiler optimization flags.

Thanks. I'll try and remember that. Delphi enforces local variable declarations before any logical operations can be performed (IIRC to allow a one pass compile), so I'm going to have to break a long standing habit.

...

TPolygon result; //it would be better if the TPolygon constructor accepted the size

As you can see I'm still coming to grips with the fundamentals of C++. I forgot that the std vector's constructor accepted a size parameter.

...

y = y + (j % 2)*2 -1; //squiggles right twices as fast as it squiggles up

That was so I could see it was really doing what I expected, otherwise it just looked like a fat blurry line :).

...

This will be the worst case input and have expected runtime complexity of O(n^2) using a linked list as the sweepline data structure while an optimal Vatti using a tree would have O(n log n) complexity and empirical scale factor of about O(n^1.05).

I've given my data structure a little more thought and believe it isn't a simple link-list structure as I stated earlier. The main linked list, the active edge list (AEL) contains a list of active edges which also link to their adjacent edges in an edge bound (LML). I imagine that this is the tree structure you been mentioning.

It doesn't sound like a tree. I don't quite grasp what you mean by adjacent edges in an edge bound. The empirical scaling factor of an optimal sweepline that uses a tree should be about n^1.05 for all test cases. That's what we see with CGAL, for example. I use the std::map (internally a red-black tree) as the sweep line data structure to maintain a sorted sequence of line segments that intersect the sweepline. I can lookup, insert and remove line segments from the map in O(log n) (amortized) runtime. Each line segment end point scanned during the progress of the sweepline gives rise to a constant factor number of lookup, insertion and removals from the sweepline. There are no more than O(n) elements in the map at any given time and there are O(n) line segment end points to process. Intersection events add to the number of events to be processed but not the max number of elements in the map. Therefore the runtime complexity is O((n+s) log n) for n line segments and s intersection points (Vatti's reported runtime complexity.) I usually simplfy to O(n log n) wrt line segments plus intersections because the number of intersections is at worst n^2 and log(n^2) = 2 log(n) so differs by only a constant factor. If the sweepline data structure requires linear lookup then you have O((n + s)n) runtime complexity. If the number of unique x values is smaller than O(n) then you could get a significant optimization with the list data structure by keeping your place in the sweepline data structure from the previous y value and a runtime complexity of O((n + s)x). There is another, more general, data dependent optimization where instead of line segments it is x-monotone chains of line segments that are processed by the sweepline, giving O(n + (x+s)log x) runtime complexity where x is the number of x montone chains. In the worst case x is n, but if you have nice smooth arcs in your data formed of large numbers of small line segments the speedup is dramatic. Implementing the x-monotone chain sweepline is similar in principle to implementing the sweepline over curves. There is a pretty significant difference between making the sub-optimal sweepline algorithm numerically robust and the optimal tree based sweepline algorithm numerically robust, so you might want to step back and decide whether you are looking to improve the robustness of your current algorithm or to improve the algorithm then come back to robustness. I have been mentoring a successful Google Summer of Code project this year that implemented an optimal and robust sweepline algorithm for the vornoi diagram of points and soon line segments. The robust and optimal sweepline for line segment intersection and polygon clipping is also achievable, but very tricky. Your work shows that you know how to test your code well, which counts for a lot. I'd like to encourage you to persevere in pursuing robustness as well as algorithmic optimality. Regards, Luke

Angus Johnson wrote:

...

There is a pretty significant difference between making the sub-optimal sweepline algorithm numerically robust and the optimal tree based sweepline algorithm numerically robust, so you might want to step back and decide whether you are looking to improve the robustness of your current algorithm or to improve the algorithm then come back to robustness.

I believe my current algorithm (which very closely adhers to Vatti's data structure) is close to maximally efficient. The data structure at least doesn't need touching but no doubt the implementation can always be improved.

Your empirical measurement of y = 1E-08x^1.4821 scaling is almost the exact scaling factor of 1.5 I would expect a sub-optimal sweepline algorithm to exhibit for the test case you described. Like I said, O(n log n) corresponds to scaling factor of less than 1.1. 1.48 is equivalent to an algorithmic complexity of O(n * sqrt(n)). 1.5 is a lot slower than 1.05 once you get over one million (like 500X slower) so this is a big issue. No matter how efficient your implementation it can't compete with an more optimal algorithm on larger tests cases. The data structure used for delivering the sorted sequence of edges to the sweepline is less important than the data structre used for the active edge list (which should be a tree and not a list). Please read the references of my boostcon paper http://www.boost.org/doc/libs/1_44_0/libs/polygon/doc/GTL_boostcon2009.pdf, in particular Ruud, B. Building a Mutable Set, 2003. Retrieved March 3, 2009, from Dr. Dobb's Portal: http://www.ddj.com/cpp/184401664 gives a clear description of how to use the std set or map to implement the sweepline data structure. Regards, Luke

Thanks for the links. Both interesting reads. A quick read has already given me a couple of ideas for improvements (though I think they will only be marginal).

I've now had a much better look at the Bentley-Ottmann algorithm and can see now why it could make a significant performance improvement. I'm also wondering if the event queue may also be manipulated to better manage complex intersection management too. Anyhow, all this would require a very major rewrite, so maybe it'll be something to look at in a few months time.

Jonathan Klein wrote:

Hi, i'm tring to serialize a std::map but get strange compiler errors.

This comes up all the time. You can't serialize a pointer to a prmiitive. Wrap the primitive in a class and serialize that. Also see BOOST_SERIALIZATION_STRONG_TYPEDEF

Also i tried to serialize a **class, which also failed.

This is not supported. There is an enhancement request to do so though. Robert Ramey

I can review your code if you like.

Luke, just a very quick follow-up to my last email. The comments in the C++ code is much briefer than that in the Delphi code so you'll probably find it helpful to look at both together.