Help regarding line intersection on sphere issue

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Help regarding line intersection on sphere issue

Stephen Hardy
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
Reply | Threaded
Open this post in threaded view
|

Re: Help regarding line intersection on sphere issue

Stephen Hardy

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
Reply | Threaded
Open this post in threaded view
|

Re: Help regarding line intersection on sphere issue

Adam Wulkiewicz
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