Modify

Ticket #8652 (closed Bugs: fixed)

Opened 11 months ago

Last modified 8 months ago

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

Reported by: flo@… Owned by: barendgehrels
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

Change History

comment:1 Changed 11 months 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 11 months ago by anonymous

Is there any workaround for this problem?

comment:3 Changed 11 months 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 8 months ago by barendgehrels

  • Status changed from new to closed
  • Resolution set to 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.

Last edited 8 months ago by barendgehrels (previous) (diff)

comment:5 Changed 8 months ago by barendgehrels

  • Milestone changed from To Be Determined to Boost 1.55.0
View

Add a comment

Modify Ticket

Change Properties
<Author field>
Action
as closed
The resolution will be deleted. Next status will be 'reopened'
Author


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

 
Note: See TracTickets for help on using tickets.