@@ -1952,6 +1952,181 @@ TEST(primal_clip, polygon_intersects_polygon)
19521952 EXPECT_NEAR (poly[9 ][0 ], 100 , EPS);
19531953 EXPECT_NEAR (poly[9 ][1 ], 250 , EPS);
19541954 }
1955+
1956+ /*
1957+ * Edge case where polygons overlap along diagonal of one polygon,
1958+ * and the edge of the other:
1959+ * +-----------+
1960+ * /| \ |
1961+ * / | \ |
1962+ * + | + |
1963+ * \ | / |
1964+ * \| / |
1965+ * +-----------+
1966+ */
1967+ {
1968+ Polygon2D subjectPolygon;
1969+ subjectPolygon.addVertex (Point2D {0 ., 0 .});
1970+ subjectPolygon.addVertex (Point2D {1 ., 0 .});
1971+ subjectPolygon.addVertex (Point2D {1 ., 1 .});
1972+ subjectPolygon.addVertex (Point2D {0 ., 1 .});
1973+
1974+ Polygon2D clipPolygon;
1975+ clipPolygon.addVertex (Point2D {0 ., 0 .});
1976+ clipPolygon.addVertex (Point2D {0.3 , 0.3 });
1977+ clipPolygon.addVertex (Point2D {0 ., 1 .});
1978+ clipPolygon.addVertex (Point2D {-0.3 , 0.7 });
1979+
1980+ Polygon2D poly = axom::primal::clip (subjectPolygon, clipPolygon, EPS);
1981+ EXPECT_NEAR (poly.signedArea (), 0.15 , EPS);
1982+ EXPECT_EQ (poly.numVertices (), 3 );
1983+
1984+ // Check vertices
1985+ EXPECT_NEAR (poly[0 ][0 ], 0 , EPS);
1986+ EXPECT_NEAR (poly[0 ][1 ], 1 , EPS);
1987+
1988+ EXPECT_NEAR (poly[1 ][0 ], 0 , EPS);
1989+ EXPECT_NEAR (poly[1 ][1 ], 0 , EPS);
1990+
1991+ EXPECT_NEAR (poly[2 ][0 ], 0.3 , EPS);
1992+ EXPECT_NEAR (poly[2 ][1 ], 0.3 , EPS);
1993+ }
1994+ }
1995+
1996+ /*
1997+ * 3x3 grid of quads that occupies the unit square,
1998+ * that clip against a lower or upper triangle:
1999+ *
2000+ * +---+---+---+
2001+ * | \ | | |
2002+ * | \| | |
2003+ * +---+---+---+
2004+ * | | \ | |
2005+ * | | \| |
2006+ * +---+---+---+
2007+ * | | | \ |
2008+ * | | | \|
2009+ * +---+---+---+
2010+ */
2011+ TEST (primal_clip, polygon_intersect_robustness)
2012+ {
2013+ using Polygon2D = axom::primal::Polygon<double , 2 >;
2014+ using Point2D = axom::primal::Point <double , 2 >;
2015+ constexpr double EPS = 1e-8 ;
2016+
2017+ const double x[] = {0 .,
2018+ 1 . / 3 .,
2019+ 2 . / 3 ,
2020+ 1 .,
2021+ 0 .,
2022+ 1 . / 3 .,
2023+ 2 . / 3 ,
2024+ 1 .,
2025+ 0 .,
2026+ 1 . / 3 .,
2027+ 2 . / 3 ,
2028+ 1 .,
2029+ 0 .,
2030+ 1 . / 3 .,
2031+ 2 . / 3 ,
2032+ 1 .};
2033+ const double y[] = {0 .,
2034+ 0 .,
2035+ 0 .,
2036+ 0 .,
2037+ 1 . / 3 .,
2038+ 1 . / 3 ,
2039+ 1 . / 3 .,
2040+ 1 . / 3 .,
2041+ 2 . / 3 .,
2042+ 2 . / 3 ,
2043+ 2 . / 3 .,
2044+ 2 . / 3 .,
2045+ 1 .,
2046+ 1 .,
2047+ 1 .,
2048+ 1 .};
2049+ const int conn[][4 ] = {{0 , 1 , 5 , 4 },
2050+ {1 , 2 , 6 , 5 },
2051+ {2 , 3 , 7 , 6 },
2052+ {4 , 5 , 9 , 8 },
2053+ {5 , 6 , 10 , 9 },
2054+ {6 , 7 , 11 , 10 },
2055+ {8 , 9 , 13 , 12 },
2056+ {9 , 10 , 14 , 13 },
2057+ {10 , 11 , 15 , 14 }};
2058+ Polygon2D fine[9 ];
2059+ for (int p = 0 ; p < 9 ; p++)
2060+ {
2061+ fine[p].addVertex (Point2D {x[conn[p][0 ]], y[conn[p][0 ]]});
2062+ fine[p].addVertex (Point2D {x[conn[p][1 ]], y[conn[p][1 ]]});
2063+ fine[p].addVertex (Point2D {x[conn[p][2 ]], y[conn[p][2 ]]});
2064+ fine[p].addVertex (Point2D {x[conn[p][3 ]], y[conn[p][3 ]]});
2065+ }
2066+
2067+ // Overlap polygons
2068+ Polygon2D clipLower;
2069+ clipLower.addVertex (Point2D {0 ., 0 .});
2070+ clipLower.addVertex (Point2D {1 ., 0 .});
2071+ clipLower.addVertex (Point2D {0 ., 1 .});
2072+
2073+ Polygon2D clipUpper;
2074+ clipUpper.addVertex (Point2D {1 ., 0 .});
2075+ clipUpper.addVertex (Point2D {1 ., 1 .});
2076+ clipUpper.addVertex (Point2D {0 ., 1 .});
2077+
2078+ // Expected VFs (volume fractions, areas) for overlaps.
2079+ const double lower_vf[] = {1 ., 1 ., 0.5 , 1 ., 0.5 , 0 ., 0.5 , 0 ., 0 .};
2080+ const double upper_vf[] = {0 ., 0 ., 0.5 , 0 ., 0.5 , 1 ., 0.5 , 1 ., 1 .};
2081+
2082+ // Each quad has an area of 1/9 (3x3 quads in the unit square).
2083+ constexpr double fine_area = (1 . / 3 .) * (1 . / 3 .);
2084+
2085+ for (int p : std::vector<int > {0 , 1 , 2 , 3 , 4 , 6 })
2086+ {
2087+ // Clip clipLower with overlapping fine quad
2088+ const auto overlapPoly = axom::primal::clip (clipLower, fine[p], EPS);
2089+ EXPECT_TRUE (overlapPoly.isValid ());
2090+
2091+ if (overlapPoly.isValid ())
2092+ {
2093+ // Area of overlapped polygon is some fraction of 1/9.
2094+ // Divide by 1/9 to get a volume fraction.
2095+ const double vf = overlapPoly.area () / fine_area;
2096+ EXPECT_NEAR (vf, lower_vf[p], EPS);
2097+ }
2098+ }
2099+
2100+ for (int p : std::vector<int > {5 , 7 , 8 })
2101+ {
2102+ // Clip clipLower with non-intersecting fine quad
2103+ const auto overlapPoly = axom::primal::clip (clipLower, fine[p], EPS);
2104+ EXPECT_FALSE (overlapPoly.isValid ());
2105+ EXPECT_EQ (overlapPoly.numVertices (), 0 );
2106+ }
2107+
2108+ for (int p : std::vector<int > {2 , 4 , 5 , 6 , 7 , 8 })
2109+ {
2110+ // Clip clipUpper with overlapping fine quad
2111+ const auto overlapPoly = axom::primal::clip (clipUpper, fine[p], EPS);
2112+ EXPECT_TRUE (overlapPoly.isValid ());
2113+
2114+ if (overlapPoly.isValid ())
2115+ {
2116+ // Area of overlapped polygon is some fraction of 1/9.
2117+ // Divide by 1/9 to get a volume fraction.
2118+ const double vf = overlapPoly.area () / fine_area;
2119+ EXPECT_NEAR (vf, upper_vf[p], EPS);
2120+ }
2121+ }
2122+
2123+ for (int p : std::vector<int > {0 , 1 , 3 })
2124+ {
2125+ // Clip clipUpper with non-intersecting fine quad
2126+ const auto overlapPoly = axom::primal::clip (clipUpper, fine[p], EPS);
2127+ EXPECT_FALSE (overlapPoly.isValid ());
2128+ EXPECT_EQ (overlapPoly.numVertices (), 0 );
2129+ }
19552130}
19562131
19572132// ------------------------------------------------------------------------------
0 commit comments