Skip to content

Commit

Permalink
Fix for divide by zero error in triangulatePolygon3D ("testDifference…
Browse files Browse the repository at this point in the history 
…Crash3" in issue Oslandia#190).
  • Loading branch information
@DRC1Spatial
DRC1Spatial committed Sep 16, 2019
1 parent 4369440 commit ecb0c76
Show file tree
Hide file tree
Showing 4 changed files with 191 additions and 15 deletions.
43 changes: 36 additions & 7 deletions src/algorithm/plane.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() );
}



/**
Expand Down
18 changes: 11 additions & 7 deletions src/triangulate/triangulatePolygon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/
Expand All @@ -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 );

Expand Down
123 changes: 122 additions & 1 deletion test/unit/SFCGAL/algorithm/DifferenceTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::string> typeNames;
Expand Down Expand Up @@ -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<Geometry> 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<Geometry> 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<Geometry> diff = algorithm::difference3D( *g1, *g2, algorithm::NoValidityCheck() );
}
catch(...)
{
}

try
{
std::auto_ptr<Geometry> diff = algorithm::difference3D( *g2, *g1, algorithm::NoValidityCheck() );
}
catch(...)
{
}
}
}

BOOST_AUTO_TEST_SUITE_END()
22 changes: 22 additions & 0 deletions test/unit/SFCGAL/algorithm/PlaneTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Polygon>() );

// See triangulatePolygon3D for this pattern
if (algorithm::hasPlane3D< Kernel >(degenerate_polygon->as<Polygon>()))
{
// Should not get here, OR plane3D with Plane3DInexactUnsafe should not divide by zero
auto div_by_zero_check = algorithm::plane3D< Kernel >( degenerate_polygon->as<Polygon>(), 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<Polygon>()) );

auto valid_plane = algorithm::plane3D< Kernel >( ok_polygon->as<Polygon>(), algorithm::Plane3DInexactUnsafe() );
}


BOOST_AUTO_TEST_SUITE_END()
Expand Down

0 comments on commit ecb0c76

Please # to comment.