Geometry Boost List: Re: [geometry] Help regarding line intersection on sphere issue

Subject: Re: [geometry] Help regarding line intersection on sphere issue
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2015-05-26 10:27:11

Hi Stephen,

Stephen Hardy wrote:
> On 17 May 2015, at 10:50 pm, Stephen Hardy <stephen.hardy_at_[hidden]> 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