[Request] Provide an argument to correct giving a vector defining the interior or exterior of a spherical polygon.

[rconde01]

Discussed in #1061

^{Originally posted by rconde01 October 4, 2022}
I'm trying to find a way to specify the inside of a polygon so that the results are correct when calling correct. For example, you might have a "cone" where the half angle is greater than 90 degrees. Currently this is corrected such that the axis of the code is flipped to the other side. Specifying a vector in the interior of the polygon should be able to resolve the ambiguity...but so far I don't see a way to do this.

[rconde01]

I guess one way to implement is do a "contains" query, and then reverse if the answer is the opposite of the expected answer. But maybe there are more efficient ways.

[vissarion]

Could you please provide a minimal code example and explain the issue in more details. I am afraid I cannot fully understand your problem with the current description.

[rconde01]

I think I can explain more without code - but if you still want code later I can provide it. First let's start with boolean operations on a 2d plane. For a polygon, there's a clear inside and outside. Therefore based on the template argument for whether rings are specified clockwise or counter clockwise, correct can easily fix outer and inner rings to be "correct". However, things are not so simple for a spherical polygon. For example, an outer ring could valid define the interior of the polygon no matter which direction it's vertices are specified in. Therefore, in order to properly perform a "correct" a point defining the intended interior is required.

[vissarion]

OK, I see your point now, thanks for clarifications. I think a simpler choice from what you propose is to change correct to do nothing for the ordering of non-cartesian polygons/rings since both orderings are correct. Now a user can easily change the ordering to contain a given point so I would not include this functionality in correct. What do you think?

[rconde01]

Here's a code example btw:
https://godbolt.org/z/6Pxbco356

A couple thoughts:

• On the one hand it's reasonable...but it seems to break user expectation because correcting vertex order is one of the main things correct is used for
• So yes, you could take your polygon, do a contains operation, and then optionally invert based on the expectation - but is that the most efficient way to achieve this?

[vissarion]

OK, my proposal above was not correct. The main issue here is that Boost.Geometry does not support polygons that cover more than half of the spheroid. So if one gives an order that is opposite from the one in polygon declaration (default: clockwise) then the polygon is invalid and correct reverse the ordering to make it valid. Note that for invalid geometries the results of operations (such as within) does not make sense and the output can be anything.

@awulkiew @barendgehrels I cannot find a point in documentation that explicitly states that those polygons are nor supported. Am I missing something here?

[rconde01]

Ouch - I definitely don't like that outcome. I admit I'm not an expert in the internals, but I don't see any reason that should be the case. Already the operations work as expected if the vertex order matches the intended definition - it's just correct which doesn't have the information to work in a meaningful way. I guess you could say something like "correct will choose the direction of minimum area". And then users could do a contains query to solve the original problem I posed. Again, I don't know if that's the most efficient way. But generally I think either representation is valid, and the user has to provide more information to disambiguate.

[vissarion]

OK I think you have to face this feature (Boost.Geometry does not support polygons that cover more than half of the spheroid) as a limitation. I agree that one can use the opposite (to the one defined by the polygon) vertex orientation to represent the complement of a polygon but this is not supported in Boost.Geometry now. Indeed if this is implemented you proposal with the extra point could be used for correct. But I do not see this as a straightforward task since all operations and algorithms should be enhanced to work with that kind of polygons.

[rconde01]

I guess I'm confused when you say it doesn't support it. Maybe correct doesn't support it, but otherwise everything appears to work as expected (although I admit I haven't thoroughly tested everything). Also, "does not support polygons that cover more than half of the spheroid"...I wonder what that means exactly. Is that a polygon whose area is greater than 50% of the total area...or a polygon whose vertices all lie on a hemisphere opposed from the center, or something else?

[vissarion]

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <cmath>
#include <iostream>
#include <vector>

namespace bg = boost::geometry;

using SphericalPoint =
    bg::model::d2::point_xy<double, bg::cs::spherical_equatorial<bg::degree>>;

using SphericalPolygon =
    boost::geometry::model::polygon<SphericalPoint,
                                    false,  // clockwise if true
                                    true,   // closed
                                    std::vector, std::vector, std::allocator,
                                    std::allocator>;

int main() 
{
    SphericalPolygon quad  {{{0,0}, {0,1}, {1,1}, {1,0}, {0,0}}};
    SphericalPoint p {0.5,0.5};

    std::cout << "within=" << bg::within(p,quad) << std::endl;
    std::cout << bg::wkt(quad) << std::endl;
    std::cout << "is valid=" << bg::is_valid(quad) << std::endl;
    std::cout << bg::area(quad) << std::endl;

    bg::correct(quad);

    std::cout << "within=" << bg::within(p,quad) << std::endl;
    std::cout << bg::wkt(quad) << std::endl;
    std::cout << "is valid=" << bg::is_valid(quad) << std::endl;
    std::cout << bg::area(quad) << std::endl;

    return 0;
}

ΟΚ let me share a simple example based on your code where I use CW ordering while the definition of the polygon was CCW. Then quad is invalid and the area result is wrong (is the negative of the correct polygon). Also if the quad is the complement of POLYGON((0 0,1 0,1 1,0 1,0 0)) then within with POINT(0.5,0.5) should return false but this is not the case. After applying correct the results are as expected.

[rconde01]

sigh - ok...that's a good example. Though, this is a very clear case - I'm still not sure for a much more complex case how the determination is made. Maybe "if all points of a polygon fall within a (freely oriented) hemisphere, then the all of the interior of that polygon must also fall within this hemisphere." It would be nice if this could be fixed, but you don't seem very amenable to that priority.

[rconde01]

Say you have a bunch of polygons which all seem valid by your definition - if you union them they might create a polygon which violates the definition of valid. It seems like this issue has to be fixed for general usage.

[vissarion]

I agree this is a interesting non-trivial feature. I would be interested to work on it but I cannot promise timelines. Also I do not know if other libraries support this feature (AFAIK Microsoft's SQL Server does; at least for certain operations).

[vissarion]

Say you have a bunch of polygons which all seem valid by your definition - if you union them they might create a polygon which violates the definition of valid. It seems like this issue has to be fixed for general usage.

In Boost Geometry union (as any other) function should not generate invalid geometries from valid ones, if you have an example that this happens please file an (different) issue.

[rconde01]

SQL Server supports boolean operations on spherical polygons?

TIL: https://learn.microsoft.com/en-us/sql/relational-databases/spatial/spatial-data-types-overview?view=sql-server-ver16

[rconde01]

In Boost Geometry union (as any other) function should not generate invalid geometries from valid ones, if you have an example that this happens please file an (different) issue.

#1097

This case just follows from the fact that the "greater than half" polygons are not handled correctly.

[rconde01]

btw if you can provide guidance, I would be interested in collaborating on fixing this

[vissarion]

sorry for the late reply, sure, I can help!

[rconde01]

great - if you can give me an outline of the tasks that need to be completed i can start digging in...would that be possible?

[rconde01]

@vissarion just pinging to remind you of this

[vissarion]

sorry for the slow replies here... I think it make sense to start with some simple algorithm like area since the whole project seems complicated (i.e. to provide full support for polygons larger than half of the globe). Then there is some issue of consistency since we do not want to change the behaviour of correct() before providing full support.
One choice would be to let some functions return right results and declare in the documentation about the status. That is, area is treating "large spherical" polygons correctly even if applying correct will change the input polygon. @barendgehrels what do you think?

[rconde01]

That seems like a good place to start to me.

[awulkiew]

@vissarion AFAIR polygons that cover less than half of the globe are there because either in setops or buffer the algorithm has to detect the order of points in rings (I don't remember why) and based on that do something with them. @barendgehrels should know more about it.

[rconde01]

@vissarion So one approach to area is pretty simple. If you assume the area should be positive, then you check for a negative result and add 4*pi. This does seem to break stuff downstream though (for example, correct).

[rconde01]

Just wanted to note that I think potentially this enhancement can apply to 2d polygons as well. You could imagine defining a "pure" hole in the 2d plane such that everything outside the exterior ring is "inside the polygon"...although i guess you could do the same by just allowing an empty exterior ring.