[geos-commits] [SCM] GEOS branch master updated. bd4e378a3e0adec55a4e36fab001039fb72f3ce2

git at osgeo.org git at osgeo.org
Thu Jul 25 19:26:41 PDT 2019


This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".

The branch, master has been updated
       via  bd4e378a3e0adec55a4e36fab001039fb72f3ce2 (commit)
       via  56417ce550e64eb4c1bf90f9959d74f34b6a296c (commit)
       via  9520ab2f0fb036990271bcb70e0fcea3d4c671ea (commit)
       via  92d9e66a58aaa2ef44412986d6fb83b5840d972e (commit)
       via  8343c2cd19f5ce7bc729d97c3e79aea46c3de85c (commit)
       via  cc91e6dd19c9fb1884caa3ae052e02e0e2bdfcf8 (commit)
       via  ac9f3d5a6a669af6e49f21a7f5d2940488e8b1b8 (commit)
       via  c95d2b06c15ef81d94a2e8a782bd6875b7a69f50 (commit)
       via  3e0ac6270eb5786b646efec5e2b7c85cb6de9633 (commit)
      from  58fcec4990e79be7ab4a4300a3035cc73f203876 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit bd4e378a3e0adec55a4e36fab001039fb72f3ce2
Merge: 56417ce cc91e6d
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 15:03:12 2019 -0700

    Merge branch 'jts-448' of github.com:pramsey/geos into jts-448


commit 56417ce550e64eb4c1bf90f9959d74f34b6a296c
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 13:43:02 2019 -0700

    Add filter to only use DD math is we are in the "danger zone" for circumcenter, namely that we have a potential equally-sided right triangle on our hands.

diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index fb93f8d..77cc51d 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -67,6 +67,8 @@ public:
     void circumcentre(Coordinate& resultPoint);
     void circumcentreDD(Coordinate& resultPoint);
 
+    bool isIsoceles();
+
     static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
 
 private:
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index c495815..008fc28 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -20,6 +20,18 @@ namespace geos {
 namespace geom { // geos::geom
 
 
+bool
+Triangle::isIsoceles()
+{
+    double len0 = p1.distance(p2);
+    double len1 = p0.distance(p2);
+    double len2 = p0.distance(p1);
+    if (len0 == len1 || len1 == len2 || len2 == len0)
+        return true;
+    else
+        return false;
+}
+
 void
 Triangle::inCentre(Coordinate& result)
 {
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index a72c6ea..c4e32c0 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
@@ -421,7 +421,10 @@ public:
         Triangle triangle(triEdges[0]->orig().getCoordinate(),
                           triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
         Coordinate cc;
-        triangle.circumcentreDD(cc);
+        if (triangle.isIsoceles())
+            triangle.circumcentreDD(cc);
+        else
+            triangle.circumcentre(cc);
 
         Vertex ccVertex(cc);
 

commit 9520ab2f0fb036990271bcb70e0fcea3d4c671ea
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 13:16:34 2019 -0700

    Add actual test case

diff --git a/tests/unit/triangulate/VoronoiTest.cpp b/tests/unit/triangulate/VoronoiTest.cpp
index a7539db..3a6948e 100644
--- a/tests/unit/triangulate/VoronoiTest.cpp
+++ b/tests/unit/triangulate/VoronoiTest.cpp
@@ -12,6 +12,7 @@
 
 #include <geos/io/WKTWriter.h>
 #include <geos/io/WKTReader.h>
+#include <geos/io/WKBReader.h>
 #include <geos/geom/GeometryCollection.h>
 #include <geos/geom/GeometryFactory.h>
 
@@ -41,15 +42,29 @@ typedef group::object object;
 
 group test_voronoidiag_group("geos::triangulate::Voronoi");
 
+std::unique_ptr<Geometry>
+readTextOrHex(const char* wkt)
+{
+    WKTReader wktreader;
+    WKBReader wkbreader;
+    std::string wktstr(wkt);
+    if ("01" == wktstr.substr(0, 2) || "00" == wktstr.substr(0, 2))
+    {
+        std::stringstream is(wkt);
+        return wkbreader.readHEX(is);
+    }
+    else
+        return wktreader.read(wktstr);
+}
+
 //helper function for funning triangulation
 void
 runVoronoi(const char* sitesWkt, const char* expectedWkt, const double tolerance)
 {
-    WKTReader reader;
     WKTWriter writer;
     geos::triangulate::VoronoiDiagramBuilder builder;
-    std::unique_ptr<Geometry> sites(reader.read(sitesWkt));
-    std::unique_ptr<Geometry> expected(reader.read(expectedWkt));
+    std::unique_ptr<Geometry> sites(readTextOrHex(sitesWkt));
+    std::unique_ptr<Geometry> expected(readTextOrHex(expectedWkt));
     std::unique_ptr<GeometryCollection> results;
     const GeometryFactory& geomFact(*GeometryFactory::getDefaultInstance());
     builder.setSites(*sites);
@@ -204,5 +219,19 @@ void object::test<9>
     runVoronoi(wkt, expected, 100);
 }
 
+//Validity issue and DD calculation
+template<>
+template<>
+void object::test<10>
+()
+{
+    const char* wkt = "0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B8
 1E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540";
+
+    const char* expected =
+        "GEOMETRYCOLLECTION EMPTY";
+
+    runVoronoi(wkt, expected, 100);
+}
+
 } // namespace tut
 

commit 92d9e66a58aaa2ef44412986d6fb83b5840d972e
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 12:17:25 2019 -0700

    Fix up mistakes in drafting and add simple unit test

diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index 26f8f2d..fb93f8d 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -65,6 +65,7 @@ public:
      * @param resultPoint the point into which to write the inCentre of the triangle
      */
     void circumcentre(Coordinate& resultPoint);
+    void circumcentreDD(Coordinate& resultPoint);
 
     static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
 
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index f477285..7d884a9 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -196,18 +196,18 @@ CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
 Coordinate
 CGAlgorithmsDD::circumcentreDD(const Coordinate& a, const Coordinate& b, const Coordinate& c)
 {
-    DD ax = DD(a.x) + DD(-c.x);
-    DD ay = DD(a.y) + DD(-c.y);
-    DD bx = DD(b.x) + DD(-c.x);
-    DD by = DD(b.y) + DD(-c.y);
+    DD ax = DD(a.x) - DD(c.x);
+    DD ay = DD(a.y) - DD(c.y);
+    DD bx = DD(b.x) - DD(c.x);
+    DD by = DD(b.y) - DD(c.y);
 
-    DD denom = detDD(ax, ay, bx, by);
+    DD denom = DD(2) * detDD(ax, ay, bx, by);
     DD asqr = (ax * ax) + (ay * ay);
     DD bsqr = (bx * bx) + (by * by);
     DD numx = detDD(ay, asqr, by, bsqr);
     DD numy = detDD(ax, asqr, bx, bsqr);
     double ccx = (DD(c.x) - (numx / denom)).ToDouble();
-    double ccy = (DD(c.y) - (numy / denom)).ToDouble();
+    double ccy = (DD(c.y) + (numy / denom)).ToDouble();
     Coordinate cc(ccx, ccy);
     return cc;
 }
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index 6552745..c495815 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -14,6 +14,7 @@
 
 #include <geos/geom/Triangle.h>
 #include <geos/geom/Coordinate.h>
+#include <geos/algorithm/CGAlgorithmsDD.h>
 
 namespace geos {
 namespace geom { // geos::geom
@@ -53,6 +54,12 @@ Triangle::circumcentre(Coordinate& result)
     result = Coordinate(ccx, ccy);
 }
 
+void
+Triangle::circumcentreDD(Coordinate& result)
+{
+    result = algorithm::CGAlgorithmsDD::circumcentreDD(p0, p1, p2);
+}
+
 /* public static */
 const Coordinate
 Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2)
@@ -63,7 +70,7 @@ Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordin
     return c;
 }
 
-
+/* private */
 double
 Triangle::det(double m00, double m01, double m10, double m11) const
 {
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index 6ea08b2..a72c6ea 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
@@ -421,7 +421,7 @@ public:
         Triangle triangle(triEdges[0]->orig().getCoordinate(),
                           triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
         Coordinate cc;
-        triangle.circumcentre(cc);
+        triangle.circumcentreDD(cc);
 
         Vertex ccVertex(cc);
 
diff --git a/src/util/math.cpp b/src/util/math.cpp
index 3931dc3..a154df6 100644
--- a/src/util/math.cpp
+++ b/src/util/math.cpp
@@ -12,7 +12,7 @@
  *
  **********************************************************************/
 
-#include "geos/util.h"
+#include <geos/util.h>
 #include <cmath>
 
 namespace geos {
@@ -115,6 +115,7 @@ rint_vc(double val)
     }
 }
 
+
 } // namespace geos.util
 } // namespace geos
 
diff --git a/tests/unit/geom/TriangleTest.cpp b/tests/unit/geom/TriangleTest.cpp
index 0ab6511..c375819 100644
--- a/tests/unit/geom/TriangleTest.cpp
+++ b/tests/unit/geom/TriangleTest.cpp
@@ -6,6 +6,8 @@
 // geos
 #include <geos/geom/Triangle.h>
 #include <geos/geom/Coordinate.h>
+#include <geos/algorithm/CGAlgorithmsDD.h>
+
 // std
 #include <cmath>
 
@@ -28,6 +30,7 @@ struct test_triangle_data {
     test_triangle_data()
         : a(3, 3), b(9, 3), c(6, 6), d(-4, -2), e(-8, -2), f(-4, -4)
     {}
+
 };
 
 typedef test_group<test_triangle_data> group;
@@ -173,4 +176,29 @@ void object::test<6>
     //  std::cout << "CicumCenter of triangle DEF:: " << c2.x << " " << c2.y << std::endl;
 }
 
+// Test circumcentreDD()
+template<>
+template<>
+void object::test<7>
+()
+{
+    using geos::geom::Triangle;
+    using geos::geom::Coordinate;
+    using geos::algorithm::CGAlgorithmsDD;
+
+    Coordinate x1(193600.80333333334, 469345.355);
+    Coordinate x2(193600.80333333334, 469345.0175);
+    Coordinate x3(193601.10666666666, 469345.0175);
+
+    Coordinate y1(193600.80333333334, 469345.355);
+    Coordinate y2(193601.10666666666, 469345.0175);
+    Coordinate y3(193601.10666666666, 469345.355);
+
+    Coordinate cc1 = CGAlgorithmsDD::circumcentreDD(x1, x2, x3);
+    Coordinate cc2 = CGAlgorithmsDD::circumcentreDD(y1, y2, y3);
+
+    ensure(cc1 == cc2);
+}
+
+
 } // namespace tut

commit 8343c2cd19f5ce7bc729d97c3e79aea46c3de85c
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Wed Jul 24 15:05:53 2019 -0700

    First draft of circumcenterDD support
    from JTS https://github.com/locationtech/jts/pull/448

diff --git a/include/geos/algorithm/CGAlgorithmsDD.h b/include/geos/algorithm/CGAlgorithmsDD.h
index 8b2cb50..7f586f9 100644
--- a/include/geos/algorithm/CGAlgorithmsDD.h
+++ b/include/geos/algorithm/CGAlgorithmsDD.h
@@ -116,10 +116,35 @@ public:
 
     static int signOfDet2x2(double dx1, double dy1, double dx2, double dy2);
 
+    static DD detDD(double x1, double y1, double x2, double y2);
+    static DD detDD(const DD& x1, const DD& y1, const DD& x2, const DD& y2);
+
+    /**
+    * Computes the circumcentre of a triangle. The circumcentre is the centre of
+    * the circumcircle, the smallest circle which encloses the triangle. It is
+    * also the common intersection point of the perpendicular bisectors of the
+    * sides of the triangle, and is the only point which has equal distance to
+    * all three vertices of the triangle.
+    * <p>
+    * The circumcentre does not necessarily lie within the triangle. For example,
+    * the circumcentre of an obtuse isosceles triangle lies outside the triangle.
+    * <p>
+    * This method uses {@link DD} extended-precision arithmetic to
+    * provide more accurate results than {@link #circumcentre(Coordinate, Coordinate, Coordinate)}
+    *
+    * @param a
+    *          a vertex of the triangle
+    * @param b
+    *          a vertex of the triangle
+    * @param c
+    *          a vertex of the triangle
+    * @return the circumcentre of the triangle
+    */
+    static geom::Coordinate circumcentreDD(const geom::Coordinate& a, const geom::Coordinate& b, const geom::Coordinate& c);
 
 protected:
 
-    static int signOfDet2x2(DD& x1, DD& y1, DD& x2, DD& y2);
+    static int signOfDet2x2(const DD& x1, const DD& y1, const DD& x2, const DD& y2);
 
 };
 
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index 2fe8a55..f477285 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -87,7 +87,7 @@ CGAlgorithmsDD::orientationIndex(const Coordinate& p1,
 }
 
 int
-CGAlgorithmsDD::signOfDet2x2(DD& x1, DD& y1, DD& x2, DD& y2)
+CGAlgorithmsDD::signOfDet2x2(const DD& x1, const DD& y1, const DD& x2, const DD& y2)
 {
     DD mx1y2(x1 * y2);
     DD my1x2(y1 * x2);
@@ -192,6 +192,39 @@ CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
     rv.y = y.ToDouble();
 }
 
+/* public static */
+Coordinate
+CGAlgorithmsDD::circumcentreDD(const Coordinate& a, const Coordinate& b, const Coordinate& c)
+{
+    DD ax = DD(a.x) + DD(-c.x);
+    DD ay = DD(a.y) + DD(-c.y);
+    DD bx = DD(b.x) + DD(-c.x);
+    DD by = DD(b.y) + DD(-c.y);
+
+    DD denom = detDD(ax, ay, bx, by);
+    DD asqr = (ax * ax) + (ay * ay);
+    DD bsqr = (bx * bx) + (by * by);
+    DD numx = detDD(ay, asqr, by, bsqr);
+    DD numy = detDD(ax, asqr, bx, bsqr);
+    double ccx = (DD(c.x) - (numx / denom)).ToDouble();
+    double ccy = (DD(c.y) - (numy / denom)).ToDouble();
+    Coordinate cc(ccx, ccy);
+    return cc;
+}
+
+/* public static */
+DD
+CGAlgorithmsDD::detDD(double x1, double y1, double x2, double y2)
+{
+    return detDD(DD(x1), DD(y1), DD(x2), DD(y2));
+}
+
+/* public static */
+DD
+CGAlgorithmsDD::detDD(const DD& x1, const DD& y1, const DD& x2, const DD& y2)
+{
+    return (x1 * y2) - (y1 * x2);
+}
 
 } // namespace geos::algorithm
 } // namespace geos

commit cc91e6dd19c9fb1884caa3ae052e02e0e2bdfcf8
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 13:43:02 2019 -0700

    Add filter to only use DD math is we are in the "danger zone" for circumcenter, namely that we have a potential equally-sided right triangle on our hands.

diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index fb93f8d..77cc51d 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -67,6 +67,8 @@ public:
     void circumcentre(Coordinate& resultPoint);
     void circumcentreDD(Coordinate& resultPoint);
 
+    bool isIsoceles();
+
     static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
 
 private:
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index c495815..008fc28 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -20,6 +20,18 @@ namespace geos {
 namespace geom { // geos::geom
 
 
+bool
+Triangle::isIsoceles()
+{
+    double len0 = p1.distance(p2);
+    double len1 = p0.distance(p2);
+    double len2 = p0.distance(p1);
+    if (len0 == len1 || len1 == len2 || len2 == len0)
+        return true;
+    else
+        return false;
+}
+
 void
 Triangle::inCentre(Coordinate& result)
 {
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index a72c6ea..c4e32c0 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
@@ -421,7 +421,10 @@ public:
         Triangle triangle(triEdges[0]->orig().getCoordinate(),
                           triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
         Coordinate cc;
-        triangle.circumcentreDD(cc);
+        if (triangle.isIsoceles())
+            triangle.circumcentreDD(cc);
+        else
+            triangle.circumcentre(cc);
 
         Vertex ccVertex(cc);
 

commit ac9f3d5a6a669af6e49f21a7f5d2940488e8b1b8
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 13:16:34 2019 -0700

    Add actual test case

diff --git a/tests/unit/triangulate/VoronoiTest.cpp b/tests/unit/triangulate/VoronoiTest.cpp
index a7539db..3a6948e 100644
--- a/tests/unit/triangulate/VoronoiTest.cpp
+++ b/tests/unit/triangulate/VoronoiTest.cpp
@@ -12,6 +12,7 @@
 
 #include <geos/io/WKTWriter.h>
 #include <geos/io/WKTReader.h>
+#include <geos/io/WKBReader.h>
 #include <geos/geom/GeometryCollection.h>
 #include <geos/geom/GeometryFactory.h>
 
@@ -41,15 +42,29 @@ typedef group::object object;
 
 group test_voronoidiag_group("geos::triangulate::Voronoi");
 
+std::unique_ptr<Geometry>
+readTextOrHex(const char* wkt)
+{
+    WKTReader wktreader;
+    WKBReader wkbreader;
+    std::string wktstr(wkt);
+    if ("01" == wktstr.substr(0, 2) || "00" == wktstr.substr(0, 2))
+    {
+        std::stringstream is(wkt);
+        return wkbreader.readHEX(is);
+    }
+    else
+        return wktreader.read(wktstr);
+}
+
 //helper function for funning triangulation
 void
 runVoronoi(const char* sitesWkt, const char* expectedWkt, const double tolerance)
 {
-    WKTReader reader;
     WKTWriter writer;
     geos::triangulate::VoronoiDiagramBuilder builder;
-    std::unique_ptr<Geometry> sites(reader.read(sitesWkt));
-    std::unique_ptr<Geometry> expected(reader.read(expectedWkt));
+    std::unique_ptr<Geometry> sites(readTextOrHex(sitesWkt));
+    std::unique_ptr<Geometry> expected(readTextOrHex(expectedWkt));
     std::unique_ptr<GeometryCollection> results;
     const GeometryFactory& geomFact(*GeometryFactory::getDefaultInstance());
     builder.setSites(*sites);
@@ -204,5 +219,19 @@ void object::test<9>
     runVoronoi(wkt, expected, 100);
 }
 
+//Validity issue and DD calculation
+template<>
+template<>
+void object::test<10>
+()
+{
+    const char* wkt = "0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B8
 1E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540";
+
+    const char* expected =
+        "GEOMETRYCOLLECTION EMPTY";
+
+    runVoronoi(wkt, expected, 100);
+}
+
 } // namespace tut
 

commit c95d2b06c15ef81d94a2e8a782bd6875b7a69f50
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 12:17:25 2019 -0700

    Fix up mistakes in drafting and add simple unit test

diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index 26f8f2d..fb93f8d 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -65,6 +65,7 @@ public:
      * @param resultPoint the point into which to write the inCentre of the triangle
      */
     void circumcentre(Coordinate& resultPoint);
+    void circumcentreDD(Coordinate& resultPoint);
 
     static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
 
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index f477285..7d884a9 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -196,18 +196,18 @@ CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
 Coordinate
 CGAlgorithmsDD::circumcentreDD(const Coordinate& a, const Coordinate& b, const Coordinate& c)
 {
-    DD ax = DD(a.x) + DD(-c.x);
-    DD ay = DD(a.y) + DD(-c.y);
-    DD bx = DD(b.x) + DD(-c.x);
-    DD by = DD(b.y) + DD(-c.y);
+    DD ax = DD(a.x) - DD(c.x);
+    DD ay = DD(a.y) - DD(c.y);
+    DD bx = DD(b.x) - DD(c.x);
+    DD by = DD(b.y) - DD(c.y);
 
-    DD denom = detDD(ax, ay, bx, by);
+    DD denom = DD(2) * detDD(ax, ay, bx, by);
     DD asqr = (ax * ax) + (ay * ay);
     DD bsqr = (bx * bx) + (by * by);
     DD numx = detDD(ay, asqr, by, bsqr);
     DD numy = detDD(ax, asqr, bx, bsqr);
     double ccx = (DD(c.x) - (numx / denom)).ToDouble();
-    double ccy = (DD(c.y) - (numy / denom)).ToDouble();
+    double ccy = (DD(c.y) + (numy / denom)).ToDouble();
     Coordinate cc(ccx, ccy);
     return cc;
 }
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index 6552745..c495815 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -14,6 +14,7 @@
 
 #include <geos/geom/Triangle.h>
 #include <geos/geom/Coordinate.h>
+#include <geos/algorithm/CGAlgorithmsDD.h>
 
 namespace geos {
 namespace geom { // geos::geom
@@ -53,6 +54,12 @@ Triangle::circumcentre(Coordinate& result)
     result = Coordinate(ccx, ccy);
 }
 
+void
+Triangle::circumcentreDD(Coordinate& result)
+{
+    result = algorithm::CGAlgorithmsDD::circumcentreDD(p0, p1, p2);
+}
+
 /* public static */
 const Coordinate
 Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2)
@@ -63,7 +70,7 @@ Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordin
     return c;
 }
 
-
+/* private */
 double
 Triangle::det(double m00, double m01, double m10, double m11) const
 {
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index 6ea08b2..a72c6ea 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
@@ -421,7 +421,7 @@ public:
         Triangle triangle(triEdges[0]->orig().getCoordinate(),
                           triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
         Coordinate cc;
-        triangle.circumcentre(cc);
+        triangle.circumcentreDD(cc);
 
         Vertex ccVertex(cc);
 
diff --git a/src/util/math.cpp b/src/util/math.cpp
index 3931dc3..a154df6 100644
--- a/src/util/math.cpp
+++ b/src/util/math.cpp
@@ -12,7 +12,7 @@
  *
  **********************************************************************/
 
-#include "geos/util.h"
+#include <geos/util.h>
 #include <cmath>
 
 namespace geos {
@@ -115,6 +115,7 @@ rint_vc(double val)
     }
 }
 
+
 } // namespace geos.util
 } // namespace geos
 
diff --git a/tests/unit/geom/TriangleTest.cpp b/tests/unit/geom/TriangleTest.cpp
index 0ab6511..c375819 100644
--- a/tests/unit/geom/TriangleTest.cpp
+++ b/tests/unit/geom/TriangleTest.cpp
@@ -6,6 +6,8 @@
 // geos
 #include <geos/geom/Triangle.h>
 #include <geos/geom/Coordinate.h>
+#include <geos/algorithm/CGAlgorithmsDD.h>
+
 // std
 #include <cmath>
 
@@ -28,6 +30,7 @@ struct test_triangle_data {
     test_triangle_data()
         : a(3, 3), b(9, 3), c(6, 6), d(-4, -2), e(-8, -2), f(-4, -4)
     {}
+
 };
 
 typedef test_group<test_triangle_data> group;
@@ -173,4 +176,29 @@ void object::test<6>
     //  std::cout << "CicumCenter of triangle DEF:: " << c2.x << " " << c2.y << std::endl;
 }
 
+// Test circumcentreDD()
+template<>
+template<>
+void object::test<7>
+()
+{
+    using geos::geom::Triangle;
+    using geos::geom::Coordinate;
+    using geos::algorithm::CGAlgorithmsDD;
+
+    Coordinate x1(193600.80333333334, 469345.355);
+    Coordinate x2(193600.80333333334, 469345.0175);
+    Coordinate x3(193601.10666666666, 469345.0175);
+
+    Coordinate y1(193600.80333333334, 469345.355);
+    Coordinate y2(193601.10666666666, 469345.0175);
+    Coordinate y3(193601.10666666666, 469345.355);
+
+    Coordinate cc1 = CGAlgorithmsDD::circumcentreDD(x1, x2, x3);
+    Coordinate cc2 = CGAlgorithmsDD::circumcentreDD(y1, y2, y3);
+
+    ensure(cc1 == cc2);
+}
+
+
 } // namespace tut

commit 3e0ac6270eb5786b646efec5e2b7c85cb6de9633
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Wed Jul 24 15:05:53 2019 -0700

    First draft of circumcenterDD support
    from JTS https://github.com/locationtech/jts/pull/448

diff --git a/include/geos/algorithm/CGAlgorithmsDD.h b/include/geos/algorithm/CGAlgorithmsDD.h
index 8b2cb50..7f586f9 100644
--- a/include/geos/algorithm/CGAlgorithmsDD.h
+++ b/include/geos/algorithm/CGAlgorithmsDD.h
@@ -116,10 +116,35 @@ public:
 
     static int signOfDet2x2(double dx1, double dy1, double dx2, double dy2);
 
+    static DD detDD(double x1, double y1, double x2, double y2);
+    static DD detDD(const DD& x1, const DD& y1, const DD& x2, const DD& y2);
+
+    /**
+    * Computes the circumcentre of a triangle. The circumcentre is the centre of
+    * the circumcircle, the smallest circle which encloses the triangle. It is
+    * also the common intersection point of the perpendicular bisectors of the
+    * sides of the triangle, and is the only point which has equal distance to
+    * all three vertices of the triangle.
+    * <p>
+    * The circumcentre does not necessarily lie within the triangle. For example,
+    * the circumcentre of an obtuse isosceles triangle lies outside the triangle.
+    * <p>
+    * This method uses {@link DD} extended-precision arithmetic to
+    * provide more accurate results than {@link #circumcentre(Coordinate, Coordinate, Coordinate)}
+    *
+    * @param a
+    *          a vertex of the triangle
+    * @param b
+    *          a vertex of the triangle
+    * @param c
+    *          a vertex of the triangle
+    * @return the circumcentre of the triangle
+    */
+    static geom::Coordinate circumcentreDD(const geom::Coordinate& a, const geom::Coordinate& b, const geom::Coordinate& c);
 
 protected:
 
-    static int signOfDet2x2(DD& x1, DD& y1, DD& x2, DD& y2);
+    static int signOfDet2x2(const DD& x1, const DD& y1, const DD& x2, const DD& y2);
 
 };
 
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index 2fe8a55..f477285 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -87,7 +87,7 @@ CGAlgorithmsDD::orientationIndex(const Coordinate& p1,
 }
 
 int
-CGAlgorithmsDD::signOfDet2x2(DD& x1, DD& y1, DD& x2, DD& y2)
+CGAlgorithmsDD::signOfDet2x2(const DD& x1, const DD& y1, const DD& x2, const DD& y2)
 {
     DD mx1y2(x1 * y2);
     DD my1x2(y1 * x2);
@@ -192,6 +192,39 @@ CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
     rv.y = y.ToDouble();
 }
 
+/* public static */
+Coordinate
+CGAlgorithmsDD::circumcentreDD(const Coordinate& a, const Coordinate& b, const Coordinate& c)
+{
+    DD ax = DD(a.x) + DD(-c.x);
+    DD ay = DD(a.y) + DD(-c.y);
+    DD bx = DD(b.x) + DD(-c.x);
+    DD by = DD(b.y) + DD(-c.y);
+
+    DD denom = detDD(ax, ay, bx, by);
+    DD asqr = (ax * ax) + (ay * ay);
+    DD bsqr = (bx * bx) + (by * by);
+    DD numx = detDD(ay, asqr, by, bsqr);
+    DD numy = detDD(ax, asqr, bx, bsqr);
+    double ccx = (DD(c.x) - (numx / denom)).ToDouble();
+    double ccy = (DD(c.y) - (numy / denom)).ToDouble();
+    Coordinate cc(ccx, ccy);
+    return cc;
+}
+
+/* public static */
+DD
+CGAlgorithmsDD::detDD(double x1, double y1, double x2, double y2)
+{
+    return detDD(DD(x1), DD(y1), DD(x2), DD(y2));
+}
+
+/* public static */
+DD
+CGAlgorithmsDD::detDD(const DD& x1, const DD& y1, const DD& x2, const DD& y2)
+{
+    return (x1 * y2) - (y1 * x2);
+}
 
 } // namespace geos::algorithm
 } // namespace geos

-----------------------------------------------------------------------

Summary of changes:
 include/geos/algorithm/CGAlgorithmsDD.h          | 27 +++++++++++++++++-
 include/geos/geom/Triangle.h                     |  3 ++
 src/algorithm/CGAlgorithmsDD.cpp                 | 35 +++++++++++++++++++++++-
 src/geom/Triangle.cpp                            | 21 +++++++++++++-
 src/triangulate/quadedge/QuadEdgeSubdivision.cpp |  5 +++-
 src/util/math.cpp                                |  3 +-
 tests/unit/geom/TriangleTest.cpp                 | 28 +++++++++++++++++++
 tests/unit/triangulate/VoronoiTest.cpp           | 35 ++++++++++++++++++++++--
 8 files changed, 149 insertions(+), 8 deletions(-)


hooks/post-receive
-- 
GEOS


More information about the geos-commits mailing list