difference returns wrong results if double is the data type

[yurik42]

The sample code:

static const char p1_wkt[] = "POLYGON((-124.19999999845255 51.901455507812500, -124.19999999460376 51.935823093966235, -123.99999999648789 51.935823093966235, -123.99999999317222 51.902564290309876, -124.19999999845255 51.901455507812500))";
static const char p2_wkt[] = "POLYGON((-123.99999999367975 51.907655109375000, -123.99999999291659 51.900000006653443, -124.19999999861555 51.900000005468293, -124.19999999792353 51.906179355468751, -123.99999999367975 51.907655109375000))";

TEST_F( SpatialOperationF, t3_double )
{
    using point = boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian>;
    using polygon = boost::geometry::model::polygon<point>;
    using polygon_vector = std::vector<polygon>;

    auto left = boost::geometry::from_wkt<polygon>( p1_wkt );
    auto right = boost::geometry::from_wkt<polygon>( p2_wkt );

    {
        THIS_IS_A_BUG();

        polygon_vector actual;
        boost::geometry::difference( left, right, actual );
        EXPECT_EQ( 0, actual.size() );
    }
    {
        polygon_vector actual;
        boost::geometry::intersection( left, right, actual );
        EXPECT_EQ( 1, actual.size() );
    }
    {
        polygon_vector actual;
        /* swap arguments */
        boost::geometry::difference( right, left, actual );
        EXPECT_EQ( 1, actual.size() );
    }
}

The same code ran with float data type produces expected results (1 polygon in difference results)

[vissarion]

I can reproduce it. For double this is a bug and should be fixed. It is probably created by inaccuracies in computation (maybe in azimuth) for almost colinear points.

Interestingly with long double it returns slightly different results i.e.
MULTIPOLYGON(((-124.2 51.9015,-124.2 51.9015,-124.2 51.9358,-124 51.9358,-124 51.9026,-124 51.9026,-124 51.9077,-124.2 51.9062,-124.2 51.9015)))

looks like a spike but if you zoom in it isn't

This is caused because the three points with longitude around -124.2 are not colinear.

EDIT: ubuntu 22.04, develop branch, gcc-9

[JohanDore]

Hi All,

I believe I hit this same issue and thought I may help by creating a few test cases showing to help:

BOOST_AUTO_TEST_CASE(Boost_Diffence_ClearAll)
{
	using PointXY = boost::geometry::model::d2::point_xy<double>;
	using Polygon = boost::geometry::model::polygon<PointXY>;
	using MultiPolygon = boost::geometry::model::multi_polygon<Polygon>;

	MultiPolygon polyA;
	read_wkt("MULTIPOLYGON((("
	    "-2.0785311235613415 -0.6304193410175202,"
      "-2.0534946127981359 -0.6304193410175202,"
      "-2.0534946127981359 -0.8867932112327471,"
      "-2.3098684830133629 -0.8867932112327471,"
      "-2.3098684830133629 -0.6554558517807265,"
      "-2.2848319722501573 -0.6554558517807265,"
      "-2.0785311235613415 -0.6554558517807265,"
      "-2.0785311235613415 -0.6304193410175202)))", polyA);

	MultiPolygon polyB;
	read_wkt("MULTIPOLYGON((("
			"-2.0785311235613420 -0.6304193410175202,"
			"-2.0534946127981359 -0.6304193410175202,"
			"-2.0534946127981359 -0.6554558517807265,"
			"-2.0785311235613420 -0.6554558517807265,"
			"-2.0785311235613420 -0.6304193410175202)))", polyB);

	MultiPolygon polyAB;
	//Clip off the small polyB appendix top/left in part polyA
	difference(polyA, polyB, polyAB);

	//Without an error we should get:
	//BOOST_CHECK_CLOSE(area(polyAB), area(polyA) - area(polyB), 0.0001);
	//BOOST_CHECK_EQUAL(polyAB.size(), 1u);
	//BOOST_CHECK(polyAB[0].inners().empty());

	//But we get this:
	BOOST_CHECK(polyAB.empty());
}

BOOST_AUTO_TEST_CASE(Boost_Diffence_ClearsTooMuch)
{
	using PointXY = boost::geometry::model::d2::point_xy<double>;
	using Polygon = boost::geometry::model::polygon<PointXY>;
	using MultiPolygon = boost::geometry::model::multi_polygon<Polygon>;

	MultiPolygon polyA;
	read_wkt("MULTIPOLYGON((("
		"-6.1723606999999996 3.2834021000000000,"
		"-6.1723606999999996 2.8006724999999992,"
		"-5.7133718999999994 2.8006724999999992,"
		"-5.7133718999999994 3.2834021000000000,"
		"-5.6896310999999997 3.2834021000000000,"
		"-5.6896310999999997 2.7769316999999996,"
		"-6.1961014999999993 2.7769316999999996,"
		"-6.1961014999999993 3.2834021000000000,"
		"-6.1723606999999996 3.2834021000000000)))", polyA);

	MultiPolygon polyB;
	read_wkt("MULTIPOLYGON((("
		"-6.1723606999999996 2.8006724999999997,"
		"-5.7133718999999994 2.8006724999999997,"
		"-5.7133718999999994 2.7769316999999996,"
		"-6.1723606999999996 2.7769316999999996,"
		"-6.1723606999999996 2.8006724999999997)))", polyB);

	MultiPolygon polyAB;
	//Clip off the small polyB at the bottom center part of polyA
	difference(polyA, polyB, polyAB);

	//Without an error we should get:
	//BOOST_CHECK_CLOSE(area(polyAB), area(polyA) - area(polyB), 0.0001);
	//BOOST_CHECK_EQUAL(polyAB.size(), 2u);
	//BOOST_CHECK(polyAB[0].inners().empty());
	//BOOST_CHECK(polyAB[1].inners().empty());

	//But we get this:
	BOOST_CHECK_CLOSE(area(polyAB), 0.5*(area(polyA) - area(polyB)), 0.0001);
	BOOST_CHECK(polyAB.size() == 1u);
	BOOST_CHECK(polyAB[0].inners().empty());
}

[barendgehrels]

Not fixed by my concept fix for #1294

And also Johan's two cases are (unfortunately) not fixed by that.

[barendgehrels]

Not fixed by my concept fix for #1226 and #1326
Neither are Johan's two cases.

[barendgehrels]

Not fixed by my concept fix for #1342 in PR #1343
Neither are Johan's two cases.

[barendgehrels]

@JohanDore I moved your two cases to #1345 because they have a different root cause.
I will create a PR with a fix for it today.

[JohanDore]

Thanks a lot @barend ***@***.***> Johan Doré CTO

[barendgehrels]

Though the case looks similar, it has a different cause than the new #1345.
In the reported case, two turns are just missed. In the cloned case, a wrong decision was made.
To be continued.