axom::primal::clip() for polygons is unreliable.

[BradWhitlock]

I was trying to clip 2 polygons (a quad and a triangle) to determine their overlap and the clip() routine sometimes makes invalid polygons, which can lead to wrong results or a crash.

#include <axom/primal.hpp>

#include <iostream>

#if 1
#include <cmath>
// Don't use these in real test
#define EXPECT_TRUE(value) if(!(value)) {std::cout << #value << " was false but should be true" << std::endl; }
#define EXPECT_NEAR(A,B,EPS) if(fabs((A)-(B))>=(EPS)) {std::cout << "Difference of " << #A << " (" << (A) << ") and " << #B << " (" << (B) << ") greater than " << (EPS) << std::endl; }
#endif

void run_clip_test()
{
  using Polygon = axom::primal::Polygon<double, 2>;
  using Point = axom::primal::Point<double, 2>;

  // 3x3 grid of quads that occupies unit square.
  const double x[] = {0., 1./3., 2./3, 1.,
                      0., 1./3., 2./3, 1.,
                      0., 1./3., 2./3, 1.,
                      0., 1./3., 2./3, 1.};
  const double y[] = {0., 0., 0., 0.,
                      1./3., 1./3, 1./3., 1./3.,
                      2./3., 2./3, 2./3., 2./3.,
                      1., 1., 1., 1.};
  const int conn[][4] = {{0,1,5,4}, {1,2,6,5}, {2,3,7,6},
                         {4,5,9,8}, {5,6,10,9}, {6,7,11,10},
                         {8,9,13,12}, {9,10,14,13}, {10,11,15,14}};
  Polygon fine[9];
  for(int p = 0; p < 9; p++)
  {
    fine[p].addVertex(Point{x[conn[p][0]], y[conn[p][0]]});
    fine[p].addVertex(Point{x[conn[p][1]], y[conn[p][1]]});
    fine[p].addVertex(Point{x[conn[p][2]], y[conn[p][2]]});
    fine[p].addVertex(Point{x[conn[p][3]], y[conn[p][3]]});
  }

  // Overlap polygons
  Polygon clipLower;
  clipLower.addVertex(Point{0., 0.});
  clipLower.addVertex(Point{1., 0.});
  clipLower.addVertex(Point{0., 1.});

  Polygon clipUpper;
  clipUpper.addVertex(Point{1., 0.});
  clipUpper.addVertex(Point{1., 1.});
  clipUpper.addVertex(Point{0., 1.});

  // Expected VFs for overlaps.
  const double lower_vf[] = {1., 1., 0.5,
                             1., 0.5, 0.,
                             0.5, 0., 0.};
  const double upper_vf[] = {0., 0., 0.5,
                             0., 0.5, 1.,
                             0.5, 1., 1.};

  constexpr double eps = 1.e-8;
  constexpr double fine_area = (1./3.) * (1./3.);
  for(int p : std::vector<int>{0,1,2,3,4,6})
  {
    // Clip clipLower with fine quad
    const auto overlapPoly = axom::primal::clip(clipLower, fine[p], eps);
    std::cout << "p=" << p
              << ", subject=" << clipLower
              << ", clip=" << fine[p]
              << ", overlap=" << overlapPoly << std::endl;
    EXPECT_TRUE(overlapPoly.isValid());

    if(overlapPoly.isValid())
    {
      const double vf = overlapPoly.area() / fine_area;
      EXPECT_NEAR(vf, lower_vf[p], eps);
    }
  }
  for(int p : std::vector<int>{2,4,5,6,7,8})
  {
    // Clip clipUpper with fine quad
    const auto overlapPoly = axom::primal::clip(clipUpper, fine[p], eps);
    std::cout << "p=" << p
              << ", subject=" << clipUpper
              << ", clip=" << fine[p]
              << ", overlap=" << overlapPoly << std::endl;
    EXPECT_TRUE(overlapPoly.isValid());

    if(overlapPoly.isValid())
    {
      const double vf = overlapPoly.area() / fine_area;
      EXPECT_NEAR(vf, lower_vf[p], eps);
    }
  }
}

int main(int argc, char *argv[])
{
  run_clip_test();
  return 0;
}

p=0, subject={3-gon:(0,0),(1,0),(0,1)}, clip={4-gon:(0,0),(0.333333,0),(0.333333,0.333333),(0,0.333333)}, overlap={4-gon:(0,0.333333),(0,0),(0.333333,0),(0.333333,0.333333)}
p=1, subject={3-gon:(0,0),(1,0),(0,1)}, clip={4-gon:(0.333333,0),(0.666667,0),(0.666667,0.333333),(0.333333,0.333333)}, overlap={4-gon:(0.333333,0.333333),(0.333333,0),(0.666667,0),(0.666667,0.333333)}
p=2, subject={3-gon:(0,0),(1,0),(0,1)}, clip={4-gon:(0.666667,0),(1,0),(1,0.333333),(0.666667,0.333333)}, overlap={3-gon:(0.666667,0),(1,0),(0,0)}
Difference of vf (0) and lower_vf[p] (0.5) greater than 1e-08
p=3, subject={3-gon:(0,0),(1,0),(0,1)}, clip={4-gon:(0,0.333333),(0.333333,0.333333),(0.333333,0.666667),(0,0.666667)}, overlap={4-gon:(0,0.666667),(0,0.333333),(0.333333,0.333333),(0.333333,0.666667)}
p=4, subject={3-gon:(0,0),(1,0),(0,1)}, clip={4-gon:(0.333333,0.333333),(0.666667,0.333333),(0.666667,0.666667),(0.333333,0.666667)}, overlap={4-gon:(0.333333,0.333333),(0.666667,0.333333),(0.666667,0.333333),(0,0)}
cliptest: ArrayBase.hpp:663: T &axom::ArrayBase<axom::primal::Point<double, 2>, 1, axom::Array<axom::primal::Point<double, 2>, 1, axom::MemorySpace::Dynamic>>::operator[](const axom::IndexType) [T = axom::primal::Point<double, 2>, DIM = 1, ArrayType = axom::Array<axom::primal::Point<double, 2>, 1, axom::MemorySpace::Dynamic>]: Assertion `inBounds(idx)' failed.
Abort (core dumped)