anonymous — Wed, 05 Jun 2013 11:35:54 GMT

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 — Wed, 05 Jun 2013 12:04:06 GMT

Is there any workaround for this problem?

anonymous — Wed, 05 Jun 2013 17:43:32 GMT

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;
	}
}

status changed; resolution set

Barend Gehrels — Tue, 27 Aug 2013 20:57:06 GMT

status new → closed
resolution → fixed

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.

milestone changed

Barend Gehrels — Tue, 27 Aug 2013 21:01:35 GMT

milestone To Be Determined → Boost 1.55.0

Boost C++ Libraries: Ticket #8652: boost::geometry::intersection fails for triangle-triangle intersection

status changed; resolution set

milestone changed