Boost Geometry / intersection() seems to return inconsistent result

Question

I have a 3D vector adapted for Boost Geometry as a 2D point, and as a ring:

BOOST_GEOMETRY_REGISTER_POINT_2D(Vector3, float, cs::cartesian, x, y)
BOOST_GEOMETRY_REGISTER_RING( std::vector< Vector3 > )

Then:

1. Draw some non-convex polygon (ring)
2. Draw line segment, that cuts the non-convex polygon and divides it into 2 (the smaller one will be most often a triangle)
3. Mirror smaller of the new 2 polygons over the line-segment

Resulting are two polygons, that are overlapping and have 1 tangent edge.

I then check for intersection of the two polygons. In 15% cases, the intersection result is empty, which is a surprise (the smaller polygon can have area 1.0f..10.f, so it's not a corner case)

std::deque< Polygon > output;
bg::intersection(bigger_Polygon, mirrored_over_cutting_lineseg_Polygon, output);
// output.size() == 0 in 15% of cases

What can be the reason? I tried doing boost::geometry::correct() on each polygon before calling intersection(), but it did not help. I am running out of ideas

EDIT::

I have tested if creating new rings, with Boost Geometry types and double storage type will help:

void my_intersection( std::vector<Vector3>& polyA, std::vector<Vector3>& polyB, std::deque< ... > & output ) {
    typedef bg::model::d2::point_xy<double> point_type;
    bg::model::ring< point_type > ringA;
    bg::model::ring< point_type > ringB;

    for( int i = 0; i < (int) polyA.size(); i ++ ) {
        bg::append( ringA, bg::make< point_type >( polyA[i].x, polyA[i].y ) );
    }

    for( int i = 0; i < (int) polyB.size(); i ++ ) {
        bg::append( ringB, bg::make< point_type >( polyB[i].x, polyB[i].y ) );
    }
    ...
}

I do two intersection() calls, for polyA, polyB (my initial float Vector3), and for ringA, ringB. Then, inconsistency appears:

A[6]( 58.20822143554688 100.0000076293945 , 89.18041229248047 100.0000076293945 , 100.0000076293945 93.08255767822266 , 100 80 , 64.98564147949219 80 , 58.20822143554688 100.0000076293945 )
B[4]( 89.18040466308594 100 , 100 93.08255004882812 , 93.72125244140625 90.17939758300781 , 89.18040466308594 100 )
INFO: ------ 1 vs 0 ------ INCONSISTENCY

The "1" means: output deque has size() == 1, so intersection occurs (this is for ringA/ringB intersection). "0" is for Vector3 -- empty result.

EDIT2:

Using boost models with float storage type causes incorrect results being returned also for the ringA & ringB calls. I have this confirmed. I got once confused that doubles do not change error's "logic", but it was because accidental removal of correct() calls. With the correct() calls and double storage type for ringA/ringB redundant rings I was unable to obtain empty intersection.

EDIT3:

Here are 5 cases in which intersection() returns:

• empty result for the initial two rings (std::vector< Vector3 >),
• correct result of size() == 1, when first creating double-typed copy of the std::vector<> rings (using boost::geometry models).

Case 1:

A[6]( 58.20822143554688 100.0000076293945 , 89.18041229248047 100.0000076293945 , 100.0000076293945 93.08255767822266 , 100 80 , 64.98564147949219 80 , 58.20822143554688 100.0000076293945 )

B[4]( 89.18040466308594 100 , 100 93.08255004882812 , 93.72125244140625 90.17939758300781 , 89.18040466308594 100 )

Case 2:

A[10]( 0 100 , 66.90238189697266 99.99999237060547 , 70.97279357910156 80 , 40 80 , 40 60 , 28.31221580505371 60 , 20 67.16078948974609 , 20 80 , 0 80 , 0 100 )

B[4]( 28.31221961975098 60.00000381469727 , 20.00000762939453 67.16079711914062 , 27.08192825317383 68.22066497802734 , 28.31221961975098 60.00000381469727 )

Case 3:

A[10]( 0 100 , 72.89675903320312 100 , 73.80842590332031 80 , 40 80 , 40 60 , 26.65167617797852 60 , 20 65.58068084716797 , 20 80 , 0 80 , 0 100 )

B[4]( 26.65167999267578 60.00000381469727 , 20.00000381469727 65.5806884765625 , 25.49577522277832 66.55047607421875 , 26.65167999267578 60.00000381469727 )

Case 4:

A[6]( 47.28099060058594 99.99999237060547 , 95.71660614013672 100 , 100 97.21295166015625 , 100 80 , 68.72442626953125 80.00000762939453 , 47.28099060058594 99.99999237060547 )

B[4]( 95.71659851074219 99.99999237060547 , 99.99998474121094 97.21293640136719 , 97.45189666748047 96.08384704589844 , 95.71659851074219 99.99999237060547 )

Case 5:

A[6]( 57.69097518920898 100 , 91.16551208496094 100 , 99.99999237060547 92.9193115234375 , 100 80 , 64.8609619140625 80 , 57.69097518920898 100 )

B[4]( 91.16550445556641 99.99999237060547 , 99.99998474121094 92.9193115234375 , 93.08920288085938 91.37748718261719 , 91.16550445556641 99.99999237060547 )

EDIT4:

Here is a function, that I use to mirror polygon over the crossing line (x0,y0)-(x1,y1). The tangent edge is created using this function -- after mirroring, point lands in the same place.

Vector3 mirror_point( Vector3 p, float x0, float y0, float x1, float y1 ) {
    float dx = x1 - x0;
    float dy = y1 - y0;

    float a = ( dx * dx - dy * dy ) / ( dx * dx + dy * dy );
    float b = 2.0f * dx * dy / ( dx * dx + dy * dy );

    float x2 = a * ( p.x - x0 ) + b * ( p.y - y0 ) + x0;
    float y2 = b * ( p.x - x0 ) - a * ( p.y - y0 ) + y0;

    return Vector3( x2, y2, p.z );
}

Answer 1

My analysis on your input:

The second polygon (starting with 24.57) is counter clockwise. Also the second polygon of the second set (starting with 90.61) is counter clockwise. So boost::geometry::correct should certainly be called. And it makes difference.

So, if I use geometry::correct, I get the following results:

1) first combination, using double: intersection area=12.3854, one geometry, 4 points 2) first combination, using float: intersection area=12.3854, one geometry, 4 points (the same) This result is identical to the results of SQL Server

3) second combination, using double: intersection area=34.7862, one geometry, 4 points 4) second combination, using float: intersection area=34.7862, one geometry, 4 points This result is identical to the results of SQL Server.

Note that in both cases the second polygon is within the first polygon (in the second case non touching, in the first case it is touching - according to SQL Server).

So all output seems correct. You mention: "eliminate the empty intersections", that has been solved recently in Boost.Geometry. That fix is not yet released in 1.52, but will be in 1.53. So if that is the specific problem, you have to use the version of Boost.Trunk.

However, that would not cause empty output.