Port locationtech/jts#812, Closes libgeosGH-527

include/geos/operation/overlayng/OverlayUtil.h

@@ -60,6 +60,7 @@ class GEOS_DLL OverlayUtil {
 
     static constexpr double SAFE_ENV_BUFFER_FACTOR = 0.1;
     static constexpr int SAFE_ENV_GRID_FACTOR = 3;
+    static constexpr double AREA_HEURISTIC_TOLERANCE = 0.1;
 
     /**
     * Computes an envelope which covers the extent of the result of
@@ -85,6 +86,14 @@ class GEOS_DLL OverlayUtil {
     */
     static bool isDisjoint(const Envelope* envA, const Envelope* envB, const PrecisionModel* pm);
 
+    static bool isLess(double v1, double v2, double tol) {
+        return v1 <= v2 * (1 + tol);
+    };
+
+    static bool isGreater(double v1, double v2, double tol) {
+        return v1 >= v2 * (1 - tol);
+    }
+
 
 public:
 
@@ -157,6 +166,25 @@ class GEOS_DLL OverlayUtil {
 
     static std::unique_ptr<Geometry> toLines(OverlayGraph* graph, bool isOutputEdges, const GeometryFactory* geomFact);
 
+    /**
+    * A heuristic check for overlay result correctness
+    * comparing the areas of the input and result.
+    * The heuristic is necessarily coarse, but it detects some obvious issues.
+    * (e.g. https://github.com/locationtech/jts/issues/798)
+    * <p>
+    * <b>Note:</b> - this check is only safe if the precision model is floating.
+    * It should also be safe for snapping noding if the distance tolerance is reasonably small.
+    * (Fixed precision models can lead to collapse causing result area to expand.)
+    *
+    * @param geom0 input geometry 0
+    * @param geom1 input geometry 1
+    * @param opCode the overlay opcode
+    * @param result the overlay result
+    * @return true if the result area is consistent
+    */
+    static bool isResultAreaConsistent(
+        const Geometry* geom0, const Geometry* geom1,
+        int opCode, const Geometry* result);
 
     /**
     * Round the key point if precision model is fixed.

src/operation/overlayng/OverlayNG.cpp

@@ -36,6 +36,7 @@
 #include <geos/geom/Location.h>
 #include <geos/geom/Geometry.h>
 #include <geos/util/Interrupt.h>
+#include <geos/util/TopologyException.h>
 
 #include <algorithm>
 
@@ -266,7 +267,22 @@ OverlayNG::computeEdgeOverlay()
     }
 
     GEOS_CHECK_FOR_INTERRUPTS();
-    return extractResult(opCode, &graph);
+    std::unique_ptr<Geometry> result = extractResult(opCode, &graph);
+
+    /**
+     * Heuristic check on result area.
+     * Catches cases where noding causes vertex to move
+     * and make topology graph area "invert".
+     */
+    if (OverlayUtil::isFloating(pm)) {
+        bool isAreaConsistent = OverlayUtil::isResultAreaConsistent(
+            inputGeom.getGeometry(0),
+            inputGeom.getGeometry(1),
+            opCode, result.get());
+        if (! isAreaConsistent)
+            throw util::TopologyException("Result area inconsistent with overlay operation");
+    }
+    return result;
 }
 
 /*private*/

src/operation/overlayng/OverlayUtil.cpp

@@ -224,6 +224,41 @@ OverlayUtil::resultDimension(int opCode, int dim0, int dim1)
     return resultDimension;
 }
 
+/* public static */
+bool
+OverlayUtil::isResultAreaConsistent(
+    const Geometry* geom0, const Geometry* geom1,
+    int opCode, const Geometry* result)
+{
+    if (geom0 == nullptr || geom1 == nullptr)
+        return true;
+
+    double areaResult = result->getArea();
+    double areaA = geom0->getArea();
+    double areaB = geom1->getArea();
+    bool isConsistent = true;
+
+    switch (opCode) {
+    case OverlayNG::INTERSECTION:
+        isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
+                    && isLess(areaResult, areaB, AREA_HEURISTIC_TOLERANCE);
+        break;
+    case OverlayNG::DIFFERENCE:
+        isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
+                    && isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
+        break;
+    case OverlayNG::SYMDIFFERENCE:
+        isConsistent = isLess(areaResult, areaA + areaB, AREA_HEURISTIC_TOLERANCE);
+        break;
+    case OverlayNG::UNION:
+        isConsistent = isLess(areaA, areaResult, AREA_HEURISTIC_TOLERANCE)
+                    && isLess(areaB, areaResult, AREA_HEURISTIC_TOLERANCE)
+                    && isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
+        break;
+    }
+    return isConsistent;
+}
+
 
 /*public static*/
 std::unique_ptr<Geometry>

tests/xmltester/tests/robust/overlay/TestOverlay-jts-798.xml

@@ -0,0 +1,46 @@
+<run>
+<desc>JTS #798 - intersection fails with original OverlayNG.
+Problem is caused by Fast Noding causing an area to "invert"
+(but still producing correct noding).
+In this case Snapping Noding works.
+Fixed by adding result area heuristic check.
+
+https://github.com/locationtech/jts/issues/798
+Also https://github.com/Toblerity/Shapely/issues/1216
+</desc>
+<precisionModel type="FLOATING"/>
+
+<case>
+<desc> JTS #798 : 1- intersection fails with original OverlayNG. </desc>
+  <a>
+POLYGON ((66697.40120137333 185279.95469107336, 66698.375 185273.625, 66697.375 185280.125, 66697.40120137333 185279.95469107336))
+    </a>
+  <b>
+POLYGON ((66710 185280, 66710 185260, 66690 185260, 66690 185280, 66710 185280))
+    </b>
+<test><op name="overlayAreaTest" arg1="A" arg2="B">true</op></test>
+</case>
+
+<case>
+<desc> JTS #798 : 2 - intersection fails with original OverlayNG. </desc>
+  <a>
+POLYGON ((61607.62679190189 190194.1555478014, 61620.5389918983 190197.4990478009, 61619.64380966499 190197.26724827816, 61607.62679190189 190194.1555478014))
+    </a>
+  <b>
+POLYGON ((61620.04309999943 190199.41420000046, 61620.41290000081 190197.98570000008, 61620.53899999708 190197.4989, 61607.62680000067 190194.15549999848, 61607.07620000094 190196.27259999886, 61620.04309999943 190199.41420000046))
+    </b>
+<test><op name="overlayAreaTest" arg1="A" arg2="B">true</op></test>
+</case>
+
+<case>
+<desc> JTS #798 : 3 - intersection fails with original OverlayNG. </desc>
+  <a>
+POLYGON ((181093.43788657014 172099.78376270586, 181093.375 172099.375, 181093.57688359552 172100.6872433709, 181093.43788657014 172099.78376270586))
+    </a>
+  <b>
+POLYGON ((181118.78476615425 172110.76716212876, 181097.49019999802 172097.5095999986, 181081.35400000215 172105.30889999866, 181118.78476615425 172110.76716212876))
+    </b>
+<test><op name="overlayAreaTest" arg1="A" arg2="B">true</op></test>
+</case>
+
+</run>

tests/xmltester/tests/robust/overlay/TestOverlay-jts-808.xml

@@ -0,0 +1,21 @@
+<run>
+<desc>JTS #808 - intersection fails with original OverlayNG.
+Fast noding fails, and then problem is caused in Snapping Noding.
+Fixed by adding result area heuristic check.
+
+https://github.com/locationtech/jts/issues/808
+</desc>
+<precisionModel type="FLOATING"/>
+
+<case>
+<desc> JTS #300 - intersection fails with original OverlayNG </desc>
+  <a>
+POLYGON ((5.171796850040617 61.63109289654953, 5.171796676245829 61.63109295594216, 5.171765430643211 61.63109513643147, 5.17171972428688 61.63110275399046, 5.171693714982271 61.63111146271694, 5.17165467365684 61.631121693785616, 5.171647369933198 61.631123368287206, 5.171585767657087 61.631155078277516, 5.17157056467294 61.63116346060828, 5.171554021565499 61.63117504109626, 5.171510591367186 61.63121921397677, 5.171472849060388 61.631287634284774, 5.171431175423897 61.631360348585645, 5.171410832239342 61.63138889993328, 5.171369882904709 61.631423704808434, 5.17134097067367 61.631445185087415, 5.171322078939662 61.63146505865625, 5.171300404098109 61.631515054192896, 5.171288682554027 61.631544777638794, 5.171272603528746 61.631585550662514, 5.171264482517947 61.63162370220873, 5.171263616533105 61.631641384465254, 5.17126323358138 61.63164860015275, 5.171270591021064 61.63164912989441, 5.171294515991664 61.63165085584459, 5.171305224690445 61.63165162836994, 5.171350771945397 61.63165336444806, 5.171357855375412 61.631631252842645, 5.171413619065784 61.63163529798867, 5.171422119715523 61.63160876219859, 5.171533647053681 61.63161685241401, 5.171542147594096 61.631590316614194, 5.17176520228535 61.63160649675632, 5.171773702623785 61.631579960937884, 5.1723871038916736 61.63162445441739, 5.172404103499401 61.63157138268005, 5.172348339756998 61.63156733794006, 5.172365339390634 61.63151426620851, 5.17242110304039 61.63151831093958, 5.1724126032782385 61.6315448468102, 5.1725241307204435 61.63155293621354, 5.172532630389954 61.63152640033399, 5.172588394089763 61.63153044499734, 5.17260539328609 61.63147737322696, 5.1726611569099346 61.631481417858794, 5.172678155946876 61.631428346076376, 5.172733919494761 61.631432390676736, 5.1727594177860885 61.63135278298381, 5.172815181211687 61.63135682754822, 5.172840679213901 61.63127721983486, 5.1730079691739075 61.63128935335251, 5.172999469995483 61.631315889271136, 5.173055233395172 61.63131993373632, 5.173046734246392 61.6313464696586, 5.173158261188462 61.631354558530205, 5.17313276391995 61.631434166319124, 5.173077000301604 61.63143012188125, 5.173060001946544 61.63148319372768, 5.173004238252239 61.6314791492583, 5.172987239737798 61.631532221092684, 5.17293147596753 61.63152817659181, 5.172905977931779 61.63160778432411, 5.172961741841018 61.631611828838345, 5.172970241156623 61.63158529292391, 5.173026005036225 61.631589337411135, 5.173034504288828 61.63156280149145, 5.173090268138785 61.631566845951625, 5.173098767328383 61.63154031002671, 5.1731545311487 61.63154435445985, 5.173163030275293 61.631517818529666, 5.173218794065968 61.631521862935784, 5.173227293129557 61.631495327000366, 5.1732830568905905 61.63149937137947, 5.173291555891169 61.63147283543881, 5.1733331657906385 61.631475853262, 5.17334589872798 61.631469721527424, 5.173350581238179 61.63146669613823, 5.173372816385254 61.63139727195077, 5.173428579994363 61.63140131626687, 5.173446141375093 61.631346484215555, 5.173450425130651 61.631321443589385, 5.1734522570350405 61.6313107352209, 5.1734524972988885 61.631185832191626, 5.173440807018409 61.631184984342674, 5.173449305772492 61.63115844838812, 5.173393542596964 61.63115440408958, 5.1734020413806805 61.6311278681387, 5.173123225984876 61.631107646329674, 5.17313172498352 61.63108111040032, 5.17268562158828 61.63104875436774, 5.172677122219062 61.63107529026145, 5.172565596444439 61.631067201018645, 5.172548597470678 61.63112027278586, 5.172492834515744 61.63111622812168, 5.1724673357910245 61.63119583575317, 5.172411572713812 61.631191791053055, 5.172352074781295 61.63137754213452, 5.172240548028423 61.63136945260412, 5.172257547706421 61.631316380888315, 5.172201784447651 61.63131233609816, 5.1722442835822795 61.63117965681724, 5.1721885205718054 61.63117561202681, 5.172197020394999 61.63114907617274, 5.172141257447527 61.63114503136419, 5.172149757300357 61.631118495513796, 5.1720939944158895 61.63111445068714, 5.172102494298358 61.631087914840414, 5.171935205884021 61.631075780238405, 5.171926705862584 61.63110231607178, 5.171796850040617 61.63109289654953))
+    </a>
+  <b>
+MULTIPOLYGON (((5.17165467365684 61.631121693785616, 5.171647369933198 61.631123368287206, 5.171585767657087 61.631155078277516, 5.17157056467294 61.63116346060828, 5.171584599055308 61.63115572259883, 5.171585767657088 61.631155078277516, 5.1716449403479 61.631124618924865, 5.171647369933206 61.631123368287206, 5.17165467365684 61.631121693785616)), ((5.1714526857206184 61.63132281631538, 5.171431175423897 61.631360348585645, 5.171410832239342 61.63138889993328, 5.171431175423889 61.631360348585666, 5.1714526857206184 61.63132281631538)), ((5.171369882904709 61.631423704808434, 5.17134097067367 61.631445185087415, 5.171322078939662 61.63146505865625, 5.171322078939656 61.63146505865626, 5.171340970673664 61.63144518508743, 5.171355403249634 61.63143446243659, 5.171369882904709 61.631423704808434)), ((5.171320993973893 61.63146756125566, 5.171300404098109 61.631515054192896, 5.171300404024736 61.63151505437896, 5.171300404098101 61.63151505419292, 5.171308477044493 61.63149643300548, 5.171320993973893 61.63146756125566)), ((5.171263616533105 61.631641384465254, 5.17126323358138 61.63164860015275, 5.171270591021064 61.63164912989441, 5.171265335707341 61.631648751507534, 5.1712632335813815 61.63164860015275, 5.1712634250572425 61.631644992309006, 5.171263616533105 61.631641384465254)), ((5.171270591021064 61.63164912989441, 5.171294515991664 61.63165085584459, 5.171305224690445 61.63165162836994, 5.171270598747127 61.63164913045177, 5.171270591021064 61.63164912989441)), ((5.171305224690445 61.63165162836994, 5.171350771945397 61.63165336444806, 5.171357855375412 61.631631252842645, 5.1713507719482 61.63165336443931, 5.171305224690445 61.63165162836994)))
+    </b>
+<test><op name="overlayAreaTest" arg1="A" arg2="B">true</op></test>
+</case>
+
+</run>