Polygon clipping in primal is prone to making bad vertices in some cases. #1713
Description
I ran across some Polygon clipping problems in primal when looking at some 2D ElviraAlgorithm output. The output polygons contained an invalid (0,0) vertex well outside the input shape. That's BAD. I was able to work around this to some degree by passing smaller tolerances to axom::primal::clip() for the Polygon case. For some other polygons, tolerances did not help.
- Is there anything bad with this data (bad clipping planes, etc)?
- Would scaling really small polygons (and their clipping plane) and then scaling the vertices back down improve things numerically?
- Does the existing polygon clipping algorithm have unhandled corner cases?
Do the following
- Fix the polygon clipping algorithm
- Include these clipping cases in the test suite
I was able to reproduce it with Axom's Debug build using the clang compiler on rzwhippet and the following program:
#include <axom/primal.hpp>
#include <iostream>
#include <limits>
using FloatType = double;
using PlaneType = axom::primal::Plane<FloatType, 2>;
using VectorType = typename PlaneType::VectorType;
using PointType = typename PlaneType::PointType;
using PolygonType = axom::primal::Polygon<FloatType, 2>;
// This function produces a valid output polygon. (needed more precision)
PolygonType clip1()
{
// In the code, clipping this polygon results in a bad polygon. Note the (0,0) vertex.
// that was with the clang compiler in Debug mode.
// clip: P=normal: [ -0.0395524 -0.999217 ] offset: 0.0905541, before={4-gon:(0.28,-0.101708),(0.28,-0.101979),(0.32,-0.105354),(0.32,-0.102458)}
// clip: after={4-gon:(0,0),(0.28,-0.101979),(0.32,-0.105354),(0.32,-0.103292)}
const auto normal = VectorType{-0.0395524, -0.999217};
const FloatType offset = 0.0905541;
const auto P = PlaneType(normal, offset, false);
PolygonType shape;
shape.addVertex(PointType{0.28,-0.101708});
shape.addVertex(PointType{0.28,-0.101979});
shape.addVertex(PointType{0.32,-0.105354});
shape.addVertex(PointType{0.32,-0.102458});
auto clippedShape = axom::primal::clip(shape, P);
std::cout << "P=" << P << std::endl;
std::cout << "shape=" << shape << std::endl;
std::cout << "clippedShape=" << clippedShape << std::endl << std::endl;
// Makes
//P=normal: [ -0.0395524 -0.999217 ] offset: 0.0905541
//shape={4-gon:(0.28,-0.101708),(0.28,-0.101979),(0.32,-0.105354),(0.32,-0.102458)}
//clippedShape={4-gon:(0.28,-0.101708),(0.28,-0.101979),(0.32,-0.105354),(0.32,-0.103292)}
return clippedShape;
}
// This function is clip1's polygon with more precision and it DOES NOT make a
// valid output polygon (unless we make the tolerance really small).
PolygonType clip2(double eps)
{
// clip: before={4-gon:(0.28,-0.1017083319161494),(0.28,-0.1019791654925216),(0.32,-0.1053541654925216),(0.32,-0.1024583319161494)}, P=normal: [ -0.03955235934655436 -0.9992174992813733 ] offset: 0.09055408441368087
// clip: after={4-gon:(0,0),(0.28,-0.1019791654925216),(0.32,-0.1053541654925216),(0.32,-0.1032916652068308)}
const auto normal = VectorType{-0.03955235934655436, -0.9992174992813733};
const FloatType offset = 0.09055408441368087;
const auto P = PlaneType(normal, offset, false);
PolygonType shape;
shape.addVertex(PointType{0.28,-0.1017083319161494});
shape.addVertex(PointType{0.28,-0.1019791654925216});
shape.addVertex(PointType{0.32,-0.1053541654925216});
shape.addVertex(PointType{0.32,-0.1024583319161494});
auto clippedShape = axom::primal::clip(shape, P, eps);
std::cout << "P=" << P << std::endl;
std::cout << "shape=" << shape << std::endl;
std::cout << "clippedShape=" << clippedShape << std::endl << std::endl;
// Makes
// P=normal: [ -0.03955235934655436 -0.9992174992813733 ] offset: 0.09055408441368087
// shape={4-gon:(0.28,-0.1017083319161494),(0.28,-0.1019791654925216),(0.32,-0.1053541654925216),(0.32,-0.1024583319161494)}
// clippedShape={4-gon:(0,0),(0.28,-0.1019791654925216),(0.32,-0.1053541654925216),(0.32,-0.1032916652068308)}
// The error is the (0,0) vertex which is well outside the input shape.
return clippedShape;
}
// This is a different polygon where clip made an invalid output polygon. Basically no tolerance saves it.
PolygonType clip3(double eps)
{
// clip: before={3-gon:(0.46841046332037361566,-0.24999999999999977796),(0.4685774857907649138,-0.24999999999999977796),(0.46822188956348015365,-0.24640811891631614339)}, P=normal: [ 0.9999820004859855116 0.0059998920029159345801 ] offset: 0.46690205915894555933
// clip: after={3-gon:(0,0),(0.4685774857907649138,-0.24999999999999977796),(0.46839968767712247821,-0.24820405945815776638)}
const auto normal = VectorType{0.9999820004859855116, 0.0059998920029159345801 };
const FloatType offset = 0.46690205915894555933;
const auto P = PlaneType(normal, offset, false);
PolygonType shape;
shape.addVertex(PointType{0.46841046332037361566,-0.24999999999999977796});
shape.addVertex(PointType{0.4685774857907649138,-0.24999999999999977796});
shape.addVertex(PointType{0.46822188956348015365,-0.24640811891631614339});
auto clippedShape = axom::primal::clip(shape, P, eps);
std::cout << "P=" << P << std::endl;
std::cout << "shape=" << shape << std::endl;
std::cout << "clippedShape=" << clippedShape << std::endl << std::endl;
// Makes
// P=normal: [ 0.9999820004859855 0.005999892002915935 ] offset: 0.4669020591589456
// shape={3-gon:(0.4684104633203736,-0.2499999999999998),(0.4685774857907649,-0.2499999999999998),(0.4682218895634802,-0.2464081189163161)}
// clippedShape={3-gon:(0,0),(0.4685774857907649,-0.2499999999999998),(0.4683996876771225,-0.2482040594581578)}
// The error is the (0,0) vertex which is well outside the input shape.
return clippedShape;
}
void banner(const std::string &name)
{
std::cout << name << std::endl;
std::cout << "================================================" << std::endl;
}
int main(int argc, char *argv[])
{
banner("clip1");
clip1();
banner("clip2");
clip2(1.e-10); // default eps - fail
clip2(1.e-14); // smaller eps - pass
banner("clip3");
clip3(1.e-14); // small eps - fail
clip3(std::numeric_limits<double>::epsilon()); // smaller eps - fail
return 0;
}