The following intersection operation returns the second hole in the green polygon.
Since I've only just started playing with boost geometry, I presume I'm doing something wrong. typedef boost::geometry::model::polygon Thanks for any help ... Angus |
Solved with ... Correct(green); I evidently forgot to 'close' the ring. Sorry for bothering ... A. |
In reply to this post by angusj
Christoph, on rereading your post I can see that my previous reply wont be helpful as evidently you did 'close' the polygons (and the link I gave was to the wrong thread too).
Anyhow, I have compiled and tested your polygons and can confirm that there appears to be a small problem with the solution to an intersection operation (which returns a single point well beyond any potential intersection).
#include <boost/geometry.hpp> #include <boost/geometry/geometries/point_xy.hpp> #include <boost/geometry/geometries/polygon.hpp> #include <boost/geometry/multi/geometries/multi_polygon.hpp> #include <boost/geometry/domains/gis/io/wkt/wkt.hpp> typedef boost::geometry::model::polygon ggl_polygon; typedef boost::geometry::model::multi_polygon ggl_polygons; ggl_polygon blue, green; boost::geometry::read_wkt( "POLYGON((0.24806946917841693 26.015444246572663, 31.751930530821582 33.984555753427337, " "32.248069469178418 30.015444246572663, 0.24806946917841693 26.015444246572663))", green); boost::geometry::read_wkt( "POLYGON((17.763942722600319 32.236057277399681, 19.192448808558737 30.807551191441263, " "16.000000000000000 30.000000000000000, 17.763942722600319 32.236057277399681))", blue); ggl_polygons output; bg::intersection(green, blue, output); cout << output; |
thanks for providing the code. However if I do that using double precision, i get a triangle as output, that has little to do with the intersection. I made a mistake posting my Data. One element had a different value. I also used correct, as you suggested, but as I expected this changes nothing. I posted my complete code below. Note that the y-value 29.884635544079821 cannot belong to the first polygon. So it seems we found two errors. The question for me is if I should develop my own intersection code? I just have to calculate the difference/intersection of a lot of triangles in 2D. The library needs to handle these cases, they ALWAYS occur in practice. Greetings, Christoph #include <boost/geometry.hpp> #include <boost/geometry/geometries/polygon.hpp> #include <boost/geometry/geometries/adapted/boost_tuple.hpp> BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian) typedef boost::tuple<double, double> point; typedef boost::geometry::model::polygon<point> polygon; void boostGiveMeError() { std::vector<polygon> outDiff, outIntersect, outInnerDiff; polygon vecPoly, lpPoly; lpPoly.outer().push_back(boost::make_tuple(0.24806946917841693, 26.015444246572663)); lpPoly.outer().push_back(boost::make_tuple(31.751930530821582 , 33.984555753427337)); lpPoly.outer().push_back(boost::make_tuple(32.248069469178418 ,30.015444246572663)); lpPoly.outer().push_back(boost::make_tuple(0.24806946917841693, 26.015444246572663)); vecPoly.outer().push_back(boost::make_tuple(17.763942722600319, 32.236057277399681)); vecPoly.outer().push_back(boost::make_tuple(19.192448808558737, 30.807551191441263)); vecPoly.outer().push_back(boost::make_tuple(16.000000000000000, 30.000000000000000)); vecPoly.outer().push_back(boost::make_tuple(17.763942722600319, 32.236057277399681)); boost::geometry::correct(vecPoly); boost::geometry::correct(lpPoly); boost::geometry::difference(vecPoly, lpPoly, outDiff); for(unsigned i=0; i < outDiff.size(); i++) boost::geometry::unique(outDiff[i]); double debugMe[5][2]; for(int i=0; i < 5; i++) { debugMe[i][0] = outDiff[0].outer()[i].get<0>(); debugMe[i][1] = outDiff[0].outer()[i].get<1>(); } } THE DEBUGGER OUTPUT: - debugMe 0x0052eddc double [5][2] - [0] 0x0052eddc double [2] [0] 15.543935884491610 double [1] 29.884635544079821 double - [1] 0x0052edec double [2] [0] 16.000000000000000 double [1] 30.000000000000000 double - [2] 0x0052edfc double [2] [0] 17.763942722600319 double [1] 32.236057277399681 double - [3] 0x0052ee0c double [2] [0] 19.192448808558737 double [1] 30.807551191441263 double - [4] 0x0052ee1c double [2] [0] 15.543935884491610 double [1] 29.884635544079821 double On 08/19/2011 11:37 PM, angusj wrote: Christoph, on rereading your post I can see that my previous reply wont be helpful as evidently you did 'close' the polygons (and the link I gave was to the wrong thread too). Anyhow, I have compiled and tested your polygons and can confirm that there appears to be a small problem with the solution to an intersection operation (which returns a single point well beyond any potential intersection).#include <boost/geometry.hpp> #include <boost/geometry/geometries/point_xy.hpp> #include <boost/geometry/geometries/polygon.hpp> #include <boost/geometry/multi/geometries/multi_polygon.hpp> #include <boost/geometry/domains/gis/io/wkt/wkt.hpp> typedef boost::geometry::model::polygon ggl_polygon; typedef boost::geometry::model::multi_polygon ggl_polygons; ggl_polygon blue, green; boost::geometry::read_wkt( "POLYGON((0.24806946917841693 26.015444246572663, 31.751930530821582 33.984555753427337, " "32.248069469178418 30.015444246572663, 0.24806946917841693 26.015444246572663))", green); boost::geometry::read_wkt( "POLYGON((17.763942722600319 32.236057277399681, 19.192448808558737 30.807551191441263, " "16.000000000000000 30.000000000000000, 17.763942722600319 32.236057277399681))", blue); ggl_polygons output; bg::intersection(green, blue, output); cout << output; _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hi Christophe,
On 20-8-2011 0:13, Christoph Keller wrote:
Right, it indeed should not. Though the output triangle in general represents the (more or less) right difference, it is apparently not precise enough for your purposes. If you use ttmath, you will get the (more) correct output. I will see if it can be improved for double as well.
It is unclear to me which is the other error. The question for me is if I should develop my own intersection code? That is not the intention. Regards, Barend _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hello Barend,
Thanks for you answer. I think higher precision will not help me. Say I have two triangles A,B,C - where B=C. I substract B from A to get M. Then I intersect M and C. Because M and C may share a border i will get numeric problems with any precision. This happens quite often in my use case. I think I will avoid this by checking first if the triangles have a parrallel border - and if they do I will take measures to avoid these problems. I have found it is extremely complicated to write code that works with FP for this kind of problem. The best approach may be adaptive precision like in CGal. If that is not an option, i would just define an epsilon (maybe about 10E-6) and if a point is less than epsilon away from a line, then we assume it is on the line. Maybe the epsilon could be adapted, before including the .hpp file. Greetings, Christoph On 08/28/2011 10:22 PM, Barend Gehrels wrote: Hi Christophe, _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hi Christophe,
Personally I've seen several cases where ttmath did the job, and double didn't. But yes, if you are going up to the limits, in the end you will get problems with ttmath as well (but note that you can define precision by a template parameter). Besides this, I did some research on your case and have to do more - it seems that it is not unsolveable (carefully stated...)
Yes, in these kind of cases this might be better. But in general, I did have problems by merging intersection points on distance - this is why that part is rewritten (about a year ago). The best approach may be adaptive precision like in CGal. We selected high precision, using an external library (now ttmath but it should be able to use another as well - in the past I used GMP and CLN as well)
See above - this is a solution in some cases but will give problems in other cases. However, we might think about this further - it could be in a strategy or policy, which could be selected at compile time (this has the same effect but is better then a required epsilon). Regards, Barend _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Maybe i could use integers (32 bit long) instead of double? I have to cut out streets (rectangles) out of a 2d mesh. The result is fed into a nav mesh generator, so it does not matter if the results are not absolutely precise, but they should not be off by a very wide margin. Thanks, Christoph On 08/29/2011 10:44 PM, Barend Gehrels wrote: Hi Christophe, _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Christoph,
As it happens the Boost.Polygon library provide polygon intersection operations that use 32 bit long instead of double and are 100% numerically robust. The library uses snap
rounding and bounds the maximum error to 1 integer unit square (you can move sqrt(2) distance diagonally at maximum due to integer snapping error). To get high precision you can scale your data. You might consider giving it a try.
Regards,
Luke
From: [hidden email] [mailto:[hidden email]] On Behalf Of Christoph Keller Sent: Monday, August 29, 2011 8:51 PM To: Generic Geometry Library Discussion Subject: Re: [ggl] Re: Help needed with Intersection operation Maybe i could use integers (32 bit long) instead of double? I have to cut out streets (rectangles) out of a 2d mesh. The result is fed into a nav mesh generator, so it does not matter if the results are not absolutely precise, but they should not be off by a very wide margin. Thanks, Christoph On 08/29/2011 10:44 PM, Barend Gehrels wrote: Hi Christophe, _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
On 30-8-2011 18:47, Simonson, Lucanus J wrote:
Hi Christoph, If integer is an option, I don't exactly understand why ttmath is not. Anyway, if you want to use integers, Boost.Polygon, as Luke suggests, is a good alternative and tested extensively with integer coordinates. Besides that, I looked further to the problem and in the end it is (of course) very simple and fixed, at least temporary helped (I want to work further on it but don't know if I've the time coming month). It was a robustness issue and also influenced the testcase of Enrico and that if Brandon (submitted long ago). I should have been more defensive at that point. So thanks for the report, this again helped to make the library better, and it is possible that you encounter more problems, if so the reports are welcome. Regards, Barend _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
On 08/30/2011 11:33 PM, Barend Gehrels wrote:
On 30-8-2011 18:47, Simonson, Lucanus J wrote:Hello Barend, The difference between integer and ttmath is that if we look at the original test case and substract B from A (lets call the result C) and then intersect C with B, then no precision would have avoided this kind of numerical problem. (The problem is that a point from B is very near the border of C in the order of epsilon). With integer we get some rounding error in each operation(I can live with that). It is ok, if I get a very small triangle from the intersecton operation, that can be thrown away. This whole series of operations is a bit strange anyway. Thanks for all your effort. Greetings, Christoph _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hi Christoph,
> The difference between integer and ttmath is that if we look at the > original test case and substract B from A (lets call the result C) and > then intersect C with B, then no precision would have avoided this > kind of numerical problem. (The problem is that a point from B is very > near the border of C in the order of epsilon). With integer we get > some rounding error in each operation(I can live with that). It is ok, > if I get a very small triangle from the intersecton operation, that > can be thrown away. > > This whole series of operations is a bit strange anyway. > > Thanks for all your effort. You are welcome ;-) Just curious if you tested the fix I submitted, is your problem now solved with that, or did you select another way and did not test (I can imagine that as well). Regards, Barend _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hello Barend,
Thanks for that fix. It works for me and I am pretty happy the Geometry module is in boost. I did not have other errors and the one I had did not occur again. Thanks! Christoph On 09/06/2011 07:17 PM, Barend Gehrels wrote: > Hi Christoph, > >> The difference between integer and ttmath is that if we look at the >> original test case and substract B from A (lets call the result C) >> and then intersect C with B, then no precision would have avoided >> this kind of numerical problem. (The problem is that a point from B >> is very near the border of C in the order of epsilon). With integer >> we get some rounding error in each operation(I can live with that). >> It is ok, if I get a very small triangle from the intersecton >> operation, that can be thrown away. >> >> This whole series of operations is a bit strange anyway. >> >> Thanks for all your effort. > > You are welcome ;-) > > Just curious if you tested the fix I submitted, is your problem now > solved with that, or did you select another way and did not test (I > can imagine that as well). > > Regards, Barend > > > _______________________________________________ > ggl mailing list > [hidden email] > http://lists.osgeo.org/mailman/listinfo/ggl > _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Hi Christoph,
OK, great to hear. Regards, Barend On 12-10-2011 2:22, Christoph Keller wrote: > Hello Barend, > > Thanks for that fix. > > It works for me and I am pretty happy the Geometry module is in boost. > I did not have other errors and the one I had did not occur again. > > Thanks! > Christoph > > On 09/06/2011 07:17 PM, Barend Gehrels wrote: >> Hi Christoph, >> >>> The difference between integer and ttmath is that if we look at the >>> original test case and substract B from A (lets call the result C) >>> and then intersect C with B, then no precision would have avoided >>> this kind of numerical problem. (The problem is that a point from B >>> is very near the border of C in the order of epsilon). With integer >>> we get some rounding error in each operation(I can live with that). >>> It is ok, if I get a very small triangle from the intersecton >>> operation, that can be thrown away. >>> >>> This whole series of operations is a bit strange anyway. >>> >>> Thanks for all your effort. >> >> You are welcome ;-) >> >> Just curious if you tested the fix I submitted, is your problem now >> solved with that, or did you select another way and did not test (I >> can imagine that as well). >> >> Regards, Barend >> >> >> _______________________________________________ >> ggl mailing list >> [hidden email] >> http://lists.osgeo.org/mailman/listinfo/ggl >> > > _______________________________________________ > ggl mailing list > [hidden email] > http://lists.osgeo.org/mailman/listinfo/ggl > -- Barend Gehrels http://about.me/barendgehrels _______________________________________________ ggl mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/ggl |
Free forum by Nabble | Edit this page |