Hi,
I have been working on a code deriving plate tectonic motions on the earth’s surface based on geophysical evidence. I recently tried replacing my spherical geometry code with boost::geometry, using the spherical_equatorial coordinate system. The code became much smaller - this is great. However, I had a few regression issues, and I’ve traced them down to some peculiarities around the poles. The following is as simple as I can get it - intersect two lines, both crossing the pole, intersecting at the pole. === https://gist.github.com/sjhardy/cf3d61c0bbdd334eaf4f === #include <iostream> #include <boost/geometry.hpp> int main(int argc, const char *argv[]) { typedef boost::geometry::model::point<double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree> > spoint; typedef boost::geometry::model::linestring<spoint> sline; sline sline1, sline2; boost::geometry::append(sline1, spoint(90.0, 80.0)); boost::geometry::append(sline1, spoint(-90.0, 80.0)); boost::geometry::append(sline2, spoint(180.0, 80.0)); boost::geometry::append(sline2, spoint(0.0, 80.0)); typedef std::vector<spoint> pointvec; pointvec res3; boost::geometry::intersection(sline1, sline2, res3); for(pointvec::iterator p = res3.begin(); p != res3.end(); p++) std::cout << p->get<0>() << " " << p->get<1>() << "\n"; return 1; } === This returns $ ./boostgeomtestcase 90 80 0 80 rather than (0 90) or equivalent, as expected. I’m using boost 1.58.0. * Is there something that I am doing wrong above that is giving me the incorrect result? * Is this a known issue (I could not find a bug for it)? Many thanks, Stephen ________________________________ The information in this e-mail may be confidential and subject to legal professional privilege and/or copyright. National ICT Australia Limited accepts no liability for any damage caused by this email or its attachments. _______________________________________________ Geometry mailing list [hidden email] http://lists.boost.org/mailman/listinfo.cgi/geometry |
On 17 May 2015, at 10:50 pm, Stephen Hardy <[hidden email]> wrote: > Hi, > > I have been working on a code deriving plate tectonic motions on the earth’s surface based on geophysical evidence. > > I recently tried replacing my spherical geometry code with boost::geometry, using the spherical_equatorial coordinate system. The code became much smaller - this is great. > > However, I had a few regression issues, and I’ve traced them down to some peculiarities around the poles. > > The following is as simple as I can get it - intersect two lines, both crossing the pole, intersecting at the pole. > > === > https://gist.github.com/sjhardy/cf3d61c0bbdd334eaf4f > === > > #include <iostream> > #include <boost/geometry.hpp> > > int main(int argc, const char *argv[]) > { > > typedef boost::geometry::model::point<double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree> > spoint; > typedef boost::geometry::model::linestring<spoint> sline; > > sline sline1, sline2; > > boost::geometry::append(sline1, spoint(90.0, 80.0)); > boost::geometry::append(sline1, spoint(-90.0, 80.0)); > boost::geometry::append(sline2, spoint(180.0, 80.0)); > boost::geometry::append(sline2, spoint(0.0, 80.0)); > > typedef std::vector<spoint> pointvec; > pointvec res3; > boost::geometry::intersection(sline1, sline2, res3); > for(pointvec::iterator p = res3.begin(); p != res3.end(); p++) > std::cout << p->get<0>() << " " << p->get<1>() << "\n"; > > return 1; > } > > === > > This returns > > $ ./boostgeomtestcase > 90 80 > 0 80 > > rather than (0 90) or equivalent, as expected. I’m using boost 1.58.0. > > * Is there something that I am doing wrong above that is giving me the incorrect result? > * Is this a known issue (I could not find a bug for it)? > > Many thanks, > Stephen > > Hi, Further to the above, I’ve now found Ticket #9162 ( https://svn.boost.org/trac/boost/ticket/9162 ) from Sep ’13 - "intersects returns incorrect result in spherical equatorial". This relates to point in polygon tests failing. I can confirm that this test still fails on boost 1.58.0. regards, Stephen ________________________________ The information in this e-mail may be confidential and subject to legal professional privilege and/or copyright. National ICT Australia Limited accepts no liability for any damage caused by this email or its attachments. _______________________________________________ Geometry mailing list [hidden email] http://lists.boost.org/mailman/listinfo.cgi/geometry |
Hi Stephen,
Stephen Hardy wrote: > On 17 May 2015, at 10:50 pm, Stephen Hardy <[hidden email]> wrote: > >> Hi, >> >> I have been working on a code deriving plate tectonic motions on the earth’s surface based on geophysical evidence. >> >> I recently tried replacing my spherical geometry code with boost::geometry, using the spherical_equatorial coordinate system. The code became much smaller - this is great. >> >> However, I had a few regression issues, and I’ve traced them down to some peculiarities around the poles. >> >> The following is as simple as I can get it - intersect two lines, both crossing the pole, intersecting at the pole. >> >> === >> https://gist.github.com/sjhardy/cf3d61c0bbdd334eaf4f >> === >> >> #include <iostream> >> #include <boost/geometry.hpp> >> >> int main(int argc, const char *argv[]) >> { >> >> typedef boost::geometry::model::point<double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree> > spoint; >> typedef boost::geometry::model::linestring<spoint> sline; >> >> sline sline1, sline2; >> >> boost::geometry::append(sline1, spoint(90.0, 80.0)); >> boost::geometry::append(sline1, spoint(-90.0, 80.0)); >> boost::geometry::append(sline2, spoint(180.0, 80.0)); >> boost::geometry::append(sline2, spoint(0.0, 80.0)); >> >> typedef std::vector<spoint> pointvec; >> pointvec res3; >> boost::geometry::intersection(sline1, sline2, res3); >> for(pointvec::iterator p = res3.begin(); p != res3.end(); p++) >> std::cout << p->get<0>() << " " << p->get<1>() << "\n"; >> >> return 1; >> } >> >> === >> >> This returns >> >> $ ./boostgeomtestcase >> 90 80 >> 0 80 >> >> rather than (0 90) or equivalent, as expected. I’m using boost 1.58.0. >> >> * Is there something that I am doing wrong above that is giving me the incorrect result? >> * Is this a known issue (I could not find a bug for it)? >> >> Many thanks, >> Stephen >> >> > > Hi, > > Further to the above, I’ve now found Ticket #9162 ( https://svn.boost.org/trac/boost/ticket/9162 ) from Sep ’13 - "intersects returns incorrect result in spherical equatorial". This relates to point in polygon tests failing. > > I can confirm that this test still fails on boost 1.58.0. Thanks for the info about the failure! I think an issue mentioned in this ticket is a different one than the one you're experiencing. The ticket is related to the spherical side and/or winding strategies. But in your case I think that the segment/segment intersection strategy should be blamed. In particular, in this case as in all other cases cart_intersect strategy is used in get_turns(). So the cartesian intersection of the segments of linestrings is found. AFAIR non-cartesian intersection of segments isn't properly supported right now. It's even probable that it's wrong that this code compiles. Anyway, could you please create a separate ticket? Or maybe would you like to contribute to the library and add the spherical CS support for intersection of 2 segments? Regards, Adam _______________________________________________ Geometry mailing list [hidden email] http://lists.boost.org/mailman/listinfo.cgi/geometry |
Free forum by Nabble | Edit this page |