diff --git a/src/algorithm/plane.h b/src/algorithm/plane.h index 5a45a535..82ec5f57 100644 --- a/src/algorithm/plane.h +++ b/src/algorithm/plane.h @@ -107,21 +107,50 @@ void plane3D( /** * Returns the oriented 3D plane of a polygon (supposed to be planar). - * @warning result is rounded to double if exact is false (avoid huge expression tree) + * May return degenerate plane. */ template < typename Kernel > -CGAL::Plane_3< Kernel > plane3D( const Polygon& polygon, bool exact = true ) +CGAL::Plane_3< Kernel > plane3D( const Polygon& polygon ) { - CGAL::Vector_3< Kernel > nrml = normal3D< Kernel >( polygon, exact ); + CGAL::Vector_3< Kernel > nrml = normal3D< Kernel >( polygon, true ); - if ( !exact ) { - const double nrm = std::sqrt( CGAL::to_double( nrml.squared_length() ) ); - nrml = CGAL::Vector_3< Kernel >( nrml.x()/nrm, nrml.y()/nrm, nrml.z()/nrm ); - } + return CGAL::Plane_3< Kernel >( polygon.exteriorRing().pointN( 0 ).toPoint_3(), nrml ); +} + +struct Plane3DInexactUnsafe +{}; + +/** + * Returns the oriented 3D plane of a polygon (supposed to be planar) - inexact version. + * @warning Will divide by zero if polygon is degenerate. + * @warning result is rounded to double (avoid huge expression tree). + */ +template < typename Kernel > +CGAL::Plane_3< Kernel > plane3D( const Polygon& polygon, const Plane3DInexactUnsafe& ) +{ + CGAL::Vector_3< Kernel > nrml = normal3D< Kernel >( polygon, false ); + + const double nrm = std::sqrt( CGAL::to_double( nrml.squared_length() ) ); + nrml = CGAL::Vector_3< Kernel >( nrml.x()/nrm, nrml.y()/nrm, nrml.z()/nrm ); return CGAL::Plane_3< Kernel >( polygon.exteriorRing().pointN( 0 ).toPoint_3(), nrml ); } +/** + * Returns the oriented 3D plane of a polygon (supposed to be planar). + * This is legacy code for SFCGAL users and should be deprecated. + * @warning result is rounded to double if exact is false (avoid huge expression tree). + * @warning Will divide by zero if polygon is degenerate. This maintains the previous behaviour. + */ +template < typename Kernel > +CGAL::Plane_3< Kernel > plane3D( const Polygon& polygon, bool exact ) +{ + if ( exact ) + return plane3D< Kernel >( polygon ); + else + return plane3D< Kernel >( polygon, Plane3DInexactUnsafe() ); +} + /** diff --git a/src/triangulate/triangulatePolygon.cpp b/src/triangulate/triangulatePolygon.cpp index 3ed8d4ff..ff1e565d 100644 --- a/src/triangulate/triangulatePolygon.cpp +++ b/src/triangulate/triangulatePolygon.cpp @@ -145,6 +145,16 @@ void triangulatePolygon3D( return ; } + /* + * Make sure polygon is not degenerate + */ + if (!algorithm::hasPlane3D< Kernel >(polygon)) + { + BOOST_THROW_EXCEPTION( Exception( + ( boost::format( "can't find plane for polygon %s" ) % polygon.asText() ).str() + ) ); + } + /* * Prepare a Constraint Delaunay Triangulation */ @@ -153,13 +163,7 @@ void triangulatePolygon3D( /* * find polygon plane */ - Kernel::Plane_3 polygonPlane = algorithm::plane3D< Kernel >( polygon, false ) ; - - if ( polygonPlane.is_degenerate() ) { - BOOST_THROW_EXCEPTION( Exception( - ( boost::format( "can't find plane for polygon %s" ) % polygon.asText() ).str() - ) ); - } + Kernel::Plane_3 polygonPlane = algorithm::plane3D< Kernel >( polygon, algorithm::Plane3DInexactUnsafe() ) ; cdt.setProjectionPlane( polygonPlane ); diff --git a/test/unit/SFCGAL/algorithm/DifferenceTest.cpp b/test/unit/SFCGAL/algorithm/DifferenceTest.cpp index c20830c8..e095ccb1 100644 --- a/test/unit/SFCGAL/algorithm/DifferenceTest.cpp +++ b/test/unit/SFCGAL/algorithm/DifferenceTest.cpp @@ -41,6 +41,8 @@ BOOST_AUTO_TEST_CASE( testDifferenceXPoint ) BOOST_CHECK( algorithm::difference( Point( 0,0 ), Point( 0,0 ) )->isEmpty() ); // A different point BOOST_CHECK( *algorithm::difference( Point( 1,0 ), Point( 0,0 ) ) == Point( 1,0 ) ); + // A different point (reversed) + BOOST_CHECK( *algorithm::difference( Point( 0,0 ), Point( 1,0 ) ) == Point( 0,0 ) ); // check difference(X, point) == X std::vector typeNames; @@ -364,5 +366,124 @@ BOOST_AUTO_TEST_CASE( testDifferencePolygonVolume ) } -BOOST_AUTO_TEST_SUITE_END() +BOOST_AUTO_TEST_CASE( testDifference3DDivideByZeroCrash ) +{ + // Correct behaviour is to throw an exception about the degenerate polygon in g1. + // Incorrect behaviour is a divide by zero error + { + std::auto_ptr g1 = io::readWkt("SOLID((((1 1 0.5,0.5 1 0.5,0.5 1.5 0.5,1 1 0.5)),\ + ((0.5 1 0.5,1 1 0.5,1 1 -1,0.5 1 0.5)),\ + ((0.5 1.5 0.5,0.5 1 0.5,0.5 1 1,0.5 1.5 0.5)),\ + ((1 1 0.5,0.5 1.5 0.5,1.5 1.5 0.5,1 1 0.5)),\ + ((1 1 -1,1 1 0.5,1 0.5 0.5,1 1 -1)),\ + ((0.5 1 0.5,1 1 -1,-1 1 -1,0.5 1 0.5)),\ + ((0.5 1 1,0.5 1 0.5,-1 1 1,0.5 1 1)),\ + ((0.5 1.5 0.5,0.5 1 1,0.5 1.5 1.5,0.5 1.5 0.5)),\ + ((1.5 1.5 0.5,0.5 1.5 0.5,1.5 1.5 1.5,1.5 1.5 0.5)),\ + ((1 1 0.5,1.5 1.5 0.5,1 0.5 0.5,1 1 0.5)),\ + ((1 1 -1,1 0.5 0.5,1 -1 -1,1 1 -1)),\ + ((-1 1 -1,1 1 -1,-1 -1 -1,-1 1 -1)),\ + ((0.5 1 0.5,-1 1 -1,-1 1 1,0.5 1 0.5)),\ + ((0.5 1 1,-1 1 1,0.5 0.5 1,0.5 1 1)),\ + ((0.5 1.5 1.5,0.5 1 1,0.5 0.5 1,0.5 1.5 1.5)),\ + ((0.5 1.5 0.5,0.5 1.5 1.5,1.5 1.5 1.5,0.5 1.5 0.5)),\ + ((1.5 1.5 0.5,1.5 1.5 1.5,1.5 0.5 1.5,1.5 1.5 0.5)),\ + ((1 0.5 0.5,1.5 1.5 0.5,1.5 0.5 0.5,1 0.5 0.5)),\ + ((1 -1 -1,1 0.5 0.5,1 0.5 0.5,1 -1 -1)),\ + ((1 1 -1,1 -1 -1,-1 -1 -1,1 1 -1)),\ + ((-1 1 -1,-1 -1 -1,-1 1 1,-1 1 -1)),\ + ((0.5 0.5 1,-1 1 1,0.5 0.5 1,0.5 0.5 1)),\ + ((0.5 1.5 1.5,0.5 0.5 1,0.5 0.5 1.5,0.5 1.5 1.5)),\ + ((1.5 1.5 1.5,0.5 1.5 1.5,1.5 0.5 1.5,1.5 1.5 1.5)),\ + ((1.5 1.5 0.5,1.5 0.5 1.5,1.5 0.5 0.5,1.5 1.5 0.5)),\ + ((1 0.5 0.5,1.5 0.5 0.5,1 0.5 0.5,1 0.5 0.5)),\ + ((1 -1 -1,1 0.5 0.5,1 0.5 0,1 -1 -1)),\ + ((-1 -1 -1,1 -1 -1,-1 -1 1,-1 -1 -1)),\ + ((-1 1 1,-1 -1 -1,-1 -1 1,-1 1 1)),\ + ((0.5 0.5 1,-1 1 1,-1 -1 1,0.5 0.5 1)),\ + ((0.5 0.5 1,0.5 0.5 1,0.5 0.5 1.5,0.5 0.5 1)),\ + ((0.5 1.5 1.5,0.5 0.5 1.5,1.5 0.5 1.5,0.5 1.5 1.5)),\ + ((1.5 0.5 0.5,1.5 0.5 1.5,1 0.5 0.5,1.5 0.5 0.5)),\ + ((1 0.5 0.5,1.5 0.5 0.5,1 0.5 0,1 0.5 0.5)),\ + ((1 -1 -1,1 0.5 0,1 -1 1,1 -1 -1)),\ + ((-1 -1 1,1 -1 -1,1 -1 1,-1 -1 1)),\ + ((0.5 0.5 1,-1 -1 1,0 0.5 1,0.5 0.5 1)),\ + ((0.5 0.5 1.5,0.5 0.5 1,0 0.5 1,0.5 0.5 1.5)),\ + ((1.5 0.5 1.5,0.5 0.5 1.5,0.5 0.5 1,1.5 0.5 1.5)),\ + ((1 0.5 0.5,1.5 0.5 1.5,1 0.5 1,1 0.5 0.5)),\ + ((1.5 0.5 0.5,1 0.5 0.5,1 0.5 0,1.5 0.5 0.5)),\ + ((1 -1 1,1 0.5 0,1 0.5 0.5,1 -1 1)),\ + ((-1 -1 1,1 -1 1,0 0.5 1,-1 -1 1)),\ + ((0.5 0.5 1.5,0 0.5 1,0.5 0.5 1,0.5 0.5 1.5)),\ + ((1.5 0.5 1.5,0.5 0.5 1,1 0.5 1,1.5 0.5 1.5)),\ + ((1 0.5 0.5,1 0.5 1,1 -1 1,1 0.5 0.5)),\ + ((0 0.5 1,1 -1 1,0.5 0.5 1,0 0.5 1)),\ + ((1 0.5 1,0.5 0.5 1,1 -1 1,1 0.5 1))))"); + std::auto_ptr g2 = io::readWkt("SOLID((((-0.5 1 1,-0.5 1 -0.5,-1 1 1,-0.5 1 1)),\ + ((-0.5 1 -0.5,-0.5 1 1,-0.5 1.5 -0.5,-0.5 1 -0.5)),\ + ((-1 1 1,-0.5 1 -0.5,-1 1 -1,-1 1 1)),\ + ((-0.5 1 1,-1 1 1,-0.5 0.5 1,-0.5 1 1)),\ + ((-0.5 1.5 -0.5,-0.5 1 1,-0.5 1.5 1.5,-0.5 1.5 -0.5)),\ + ((-0.5 1 -0.5,-0.5 1.5 -0.5,1 1 -0.5,-0.5 1 -0.5)),\ + ((-1 1 -1,-0.5 1 -0.5,1 1 -1,-1 1 -1)),\ + ((-1 1 1,-1 1 -1,-1 -1 -1,-1 1 1)),\ + ((-0.5 0.5 1,-1 1 1,-0.5 -0.5 1,-0.5 0.5 1)),\ + ((-0.5 1 1,-0.5 0.5 1,-0.5 1.5 1.5,-0.5 1 1)),\ + ((-0.5 1.5 -0.5,-0.5 1.5 1.5,1.5 1.5 1.5,-0.5 1.5 -0.5)),\ + ((1 1 -0.5,-0.5 1.5 -0.5,1.5 1.5 -0.5,1 1 -0.5)),\ + ((-0.5 1 -0.5,1 1 -0.5,1 1 -1,-0.5 1 -0.5)),\ + ((-1 1 -1,1 1 -1,-1 -1 -1,-1 1 -1)),\ + ((-1 1 1,-1 -1 -1,-1 -1 1,-1 1 1)),\ + ((-0.5 -0.5 1,-1 1 1,-1 -1 1,-0.5 -0.5 1)),\ + ((-0.5 0.5 1,-0.5 -0.5 1,-0.5 -0.5 1.5,-0.5 0.5 1)),\ + ((-0.5 1.5 1.5,-0.5 0.5 1,-0.5 -0.5 1.5,-0.5 1.5 1.5)),\ + ((1.5 1.5 1.5,-0.5 1.5 1.5,1.5 -0.5 1.5,1.5 1.5 1.5)),\ + ((-0.5 1.5 -0.5,1.5 1.5 1.5,1.5 1.5 -0.5,-0.5 1.5 -0.5)),\ + ((1 1 -0.5,1.5 1.5 -0.5,1 0.5 -0.5,1 1 -0.5)),\ + ((1 1 -1,1 1 -0.5,1 0.5 -0.5,1 1 -1)),\ + ((-1 -1 -1,1 1 -1,1 -1 -1,-1 -1 -1)),\ + ((-1 -1 1,-1 -1 -1,1 -1 -1,-1 -1 1)),\ + ((-0.5 -0.5 1,-1 -1 1,0 -0.5 1,-0.5 -0.5 1)),\ + ((-0.5 -0.5 1.5,-0.5 -0.5 1,0 -0.5 1,-0.5 -0.5 1.5)),\ + ((-0.5 1.5 1.5,-0.5 -0.5 1.5,1.5 -0.5 1.5,-0.5 1.5 1.5)),\ + ((1.5 1.5 1.5,1.5 -0.5 1.5,1.5 1.5 -0.5,1.5 1.5 1.5)),\ + ((1 0.5 -0.5,1.5 1.5 -0.5,1.5 -0.5 -0.5,1 0.5 -0.5)),\ + ((1 1 -1,1 0.5 -0.5,1 -1 -1,1 1 -1)),\ + ((-1 -1 1,1 -1 -1,1 -1 1,-1 -1 1)),\ + ((0 -0.5 1,-1 -1 1,1 -1 1,0 -0.5 1)),\ + ((-0.5 -0.5 1.5,0 -0.5 1,0.5 -0.5 1,-0.5 -0.5 1.5)),\ + ((1.5 -0.5 1.5,-0.5 -0.5 1.5,0.5 -0.5 1,1.5 -0.5 1.5)),\ + ((1.5 1.5 -0.5,1.5 -0.5 1.5,1.5 -0.5 -0.5,1.5 1.5 -0.5)),\ + ((1 0.5 -0.5,1.5 -0.5 -0.5,1 -0.5 -0.5,1 0.5 -0.5)),\ + ((1 -1 -1,1 0.5 -0.5,1 -0.5 -0.5,1 -1 -1)),\ + ((1 -1 1,1 -1 -1,1 -0.5 -0,1 -1 1)),\ + ((0 -0.5 1,1 -1 1,0.5 -0.5 1,0 -0.5 1)),\ + ((1.5 -0.5 1.5,0.5 -0.5 1,1 -0.5 1,1.5 -0.5 1.5)),\ + ((1.5 -0.5 -0.5,1.5 -0.5 1.5,1 -0.5 0.5,1.5 -0.5 -0.5)),\ + ((1 -0.5 -0.5,1.5 -0.5 -0.5,1 -0.5 -0,1 -0.5 -0.5)),\ + ((1 -1 -1,1 -0.5 -0.5,1 -0.5 -0,1 -1 -1)),\ + ((1 -1 1,1 -0.5 -0,1 -0.5 0.5,1 -1 1)),\ + ((0.5 -0.5 1,1 -1 1,1 -0.5 1,0.5 -0.5 1)),\ + ((1.5 -0.5 1.5,1 -0.5 1,1 -0.5 0.5,1.5 -0.5 1.5)),\ + ((1.5 -0.5 -0.5,1 -0.5 0.5,1 -0.5 -0,1.5 -0.5 -0.5)),\ + ((1 -1 1,1 -0.5 0.5,1 -0.5 1,1 -1 1))))"); + + try + { + std::auto_ptr diff = algorithm::difference3D( *g1, *g2, algorithm::NoValidityCheck() ); + } + catch(...) + { + } + try + { + std::auto_ptr diff = algorithm::difference3D( *g2, *g1, algorithm::NoValidityCheck() ); + } + catch(...) + { + } + } +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/test/unit/SFCGAL/algorithm/PlaneTest.cpp b/test/unit/SFCGAL/algorithm/PlaneTest.cpp index 387785fa..d1d08674 100644 --- a/test/unit/SFCGAL/algorithm/PlaneTest.cpp +++ b/test/unit/SFCGAL/algorithm/PlaneTest.cpp @@ -78,6 +78,28 @@ BOOST_AUTO_TEST_CASE( testPlane ) } } +BOOST_AUTO_TEST_CASE( testPlane3DDivideByZeroCrash ) +{ + std::auto_ptr< Geometry > degenerate_polygon = io::readWkt("POLYGON((1 -1 -1,1 0.5 0.5,1 0.5 0.5,1 -1 -1))"); + BOOST_CHECK( degenerate_polygon->geometryTypeId() == TYPE_POLYGON ); + + // Should return degenerate plane without throwing + auto degenerate_plane = algorithm::plane3D< Kernel >( degenerate_polygon->as() ); + + // See triangulatePolygon3D for this pattern + if (algorithm::hasPlane3D< Kernel >(degenerate_polygon->as())) + { + // Should not get here, OR plane3D with Plane3DInexactUnsafe should not divide by zero + auto div_by_zero_check = algorithm::plane3D< Kernel >( degenerate_polygon->as(), algorithm::Plane3DInexactUnsafe() ); + } + + std::auto_ptr< Geometry > ok_polygon = io::readWkt("POLYGON((1 0.5 0.5,1.5 1.5 0.5,1.5 0.5 0.5,1 0.5 0.5))"); + BOOST_CHECK( ok_polygon->geometryTypeId() == TYPE_POLYGON ); + + BOOST_CHECK( algorithm::hasPlane3D< Kernel >(ok_polygon->as()) ); + + auto valid_plane = algorithm::plane3D< Kernel >( ok_polygon->as(), algorithm::Plane3DInexactUnsafe() ); +} BOOST_AUTO_TEST_SUITE_END()