Modify

Opened 4 years ago

Closed 4 years ago

Last modified 4 years ago

#8652 closed Bugs (fixed)

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

Reported by: flo@… Owned by: Barend Gehrels
Milestone: Boost 1.55.0 Component: geometry
Version: Boost 1.53.0 Severity: Problem
Keywords: Cc:

Description

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)

Attachments (0)

Change History (5)

comment:1 Changed 4 years ago by 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.

comment:2 Changed 4 years ago by anonymous

Is there any workaround for this problem?

comment:3 Changed 4 years ago by 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;
	}
}

comment:4 Changed 4 years ago by Barend Gehrels

Resolution: fixed
Status: newclosed

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.

Last edited 4 years ago by Barend Gehrels (previous) (diff)

comment:5 Changed 4 years ago by Barend Gehrels

Milestone: To Be DeterminedBoost 1.55.0

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain Barend Gehrels.
The resolution will be deleted.

Add Comment


E-mail address and name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.