boost::geometry::intersection fails for triangle-triangle intersection

[flo@…]

The following example fails. The computed intersection polygon is supposed to be identical to the smaller of the two triangles. However, the computed intersection polygon is more or less the outline of the larger of the two triangles with an extra collinear point.

#include <iostream>
#include <deque>

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

#include <boost/foreach.hpp>


int main()
{
    typedef
boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> >
polygon;

    polygon green, blue;

	 #define SMALL_VAL "1e-8"

    boost::geometry::read_wkt( "POLYGON((0 0, 0.05 0.04, 0.05 0, 0 0))", green);
  
    boost::geometry::read_wkt( "POLYGON((0.02 -2.77556e-17, 0.05 0.02, 0.05 -2.77556e-17, 0.02 -2.77556e-17))", blue);

	 //#define SMALL_VAL "1e-8"
	 //boost::geometry::read_wkt( "POLYGON((0.02 " SMALL_VAL ", 0.05 0.02, 0.05 " SMALL_VAL ", 0.02 " SMALL_VAL "))", blue);

    std::deque<polygon> output;
    boost::geometry::intersection(green, blue, output);

    int i = 0;
    std::cout << "The computed difference is:" << std::endl;
    BOOST_FOREACH(polygon const& p, output)
    {
		std::cout << boost::geometry::dsv(p) << std::endl;
    }

	 std::cout << std::endl;
	 std::cout << "The expected output would be: " << std::endl;
	 std::cout << "(((0.05, 0.02), (0.05, 0), (0.02, 0), (0.05, 0.02)))" << std::endl;
    
    return 0;
}

The output (compiled with gcc version 4.7.2 (Debian 4.7.2-5) on a 64bit debian wheezy system) is:

The computed difference is: (((0.05, 0.02), (0.05, -2.77556e-17), (0, 0), (0.05, 0.04), (0.05, 0.02)))

The expected output would be: (((0.05, 0.02), (0.05, 0), (0.02, 0), (0.05, 0.02)))

g++ -v

Using built-in specs. COLLECT_GCC=g++ COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/4.7/lto-wrapper Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion='Debian 4.7.2-5' --with-bugurl=​file:///usr/share/doc/gcc-4.7/README.Bugs --enable-languages=c,c++,go,fortran,objc,obj-c++ --prefix=/usr --program-suffix=-4.7 --enable-shared --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.7 --libdir=/usr/lib --enable-nls --with-sysroot=/ --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --enable-gnu-unique-object --enable-plugin --enable-objc-gc --with-arch-32=i586 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu Thread model: posix gcc version 4.7.2 (Debian 4.7.2-5)

[anonymous]

Another funny effect: Uncomment the two SMALL_VAL lines and use them as a replacement for polygon blue. If the value of SMALL_VAL is set to 1e-15, no intersection is detected at all.

[anonymous]

Is there any workaround for this problem?

[anonymous]

Mapping the floating point coordinates to integer values does not solve the problem either (see the following test program). The intersection remains empty.

#include <iostream>
#include <deque>
#include <limits>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/assign/list_of.hpp>
#include <boost/geometry/geometries/adapted/boost_tuple.hpp>
#include <boost/geometry/views/identity_view.hpp>
#include <boost/geometry/core/closure.hpp> 
#include <boost/geometry/core/exterior_ring.hpp>

#include <boost/foreach.hpp>

#define INF__ (std::numeric_limits<double>::infinity())

inline
void mapPointToInt(const double &x, const double &xOffset, const double &xSize,
int &xi)
{
	xi = ((x - xOffset)/xSize)*INT_MAX;
}

inline
void mapPointToFloatingPoint(int xI, const double &xOffset, const double &xSize, double &x) 
{
	x = (double(xI)/double(INT_MAX))*xSize + xOffset;
}

int main()
{
    typedef
boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<int> >
polygon;
	typedef boost::geometry::model::d2::point_xy<int> Point;

    polygon green, blue;

	 double greenPoints[4][2] = {{ 0, 0}, {0.05, 0.04}, {0.05, 0}, {0, 0}};
	 double bluePoints[4][2] = {{ 0.02, -2.77556e-17}, {0.05, 0.02}, {0.05 -2.77556e-17}, {0.02 -2.77556e-17}};

	 // Get the maximum extension
	 //
	 double min[2] = { INF__, INF__ };
	 double max[2] = { -INF__, -INF__ };

	 for(int i = 0; i < 3; ++i ) {
		for(int j = 0; j < 2; ++j) {
			min[j] = (greenPoints[i][j] < min[j]) ? greenPoints[i][j] : min[j];
			min[j] = (bluePoints[i][j] < min[j]) ? bluePoints[i][j] : min[j];

			max[j] = (greenPoints[i][j] > max[j]) ? greenPoints[i][j] : max[j];
			max[j] = (bluePoints[i][j] > max[j]) ? bluePoints[i][j] : max[j];
		}
	}

	std::cout << "x-range: " << min[0] << " - " << max[0] << std::endl;
	std::cout << "y-range: " << min[1] << " - " << max[1] << std::endl;

	double size[2] = { max[0] - min[0], max[1] - min[1] };

	for(int i = 0; i < 4; ++i) {

		{
			int xI, yI;
			mapPointToInt(greenPoints[i][0], min[0], size[0], xI);
			mapPointToInt(greenPoints[i][1], min[1], size[1], yI);

			boost::geometry::append(green, boost::geometry::make<Point>(xI, yI));
		}

		{
			int xI, yI;
			mapPointToInt(bluePoints[i][0], min[0], size[0], xI);
			mapPointToInt(bluePoints[i][1], min[1], size[1], yI);

			boost::geometry::append(blue, boost::geometry::make<Point>(xI, yI));
		}
	}

	boost::geometry::correct(green);
	boost::geometry::correct(blue);

	std::cout << "green poly: " <<
		boost::geometry::dsv(green) << std::endl;
	std::cout << "blue poly: " <<
		boost::geometry::dsv(blue) << std::endl;

    std::deque<polygon> output;
    boost::geometry::intersection(green, blue, output);

    int i = 0;
    std::cout << "The computed difference is:" << std::endl;
    BOOST_FOREACH(polygon const& p, output)
    {
		std::cout << boost::geometry::dsv(p) << std::endl;
    }
	typedef boost::geometry::identity_view<
 				typename boost::geometry::ring_type<polygon>::type//, 
 				//boost::geometry::open
 	> view_type;

	view_type view(boost::geometry::exterior_ring(output[0]));
 	for(auto it = boost::begin(view); it != boost::end(view); ++it) {

		double x, y;
		
		mapPointToFloatingPoint(
				boost::geometry::get<0>(*it),
				min[0], 
				size[0], 
				x);
		mapPointToFloatingPoint(
				boost::geometry::get<1>(*it),
				min[1],
				size[1],
				y);
		std::cout << "   (" << x << ", " << y << ")" << std::endl;
	}
}

[Barend Gehrels]

The boundaries of double are reached by using coordinates as specified. This results in rounding errors and imprecise results.

One of the imprecise points which was generated is fixed in commit 85451 (will be released in Boost 1.55).

The output was: POLYGON((0.05 0.02,0.05 -2.77556e-17,0.02 -2.77556e-17,0.05 0.02)) Area: 0.0003

The new output is: POLYGON((0.05 0.02,0.05 0,0.02 -2.77556e-17,0.05 0.02)) Area: 0.0003

So still not as your expectation but it comes closer.

If you use "long double" (better suited for cases like this), the output is correct: The computed intersection is: POLYGON((0.05 0.02,0.05 0,0.02 0,0.05 0.02)) Area: 0.0003

(I changed output to WKT and added area in your testprogram).

I will add a testcase for this ticket to avoid future regressions.

SQL Server gives: POLYGON ((0.019999999552965164 0, 0.049999999813735485 0, 0.049999999813735485 0.019999999552965164, 0.019999999552965164 0)) Area: 0.000299999995902181

I did not check your integer-sample, if that is a separate problem, please create a separate ticket.

Thanks for your report.