From d5fbb739287d1bb23ef0e98f475326c852881b20 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 10 Feb 2025 12:47:33 -0800 Subject: [PATCH 1/5] Add OrientationIndexStressTest --- benchmarks/algorithm/CMakeLists.txt | 3 + .../algorithm/OrientationIndexStressTest.cpp | 160 ++++++++++++++++++ 2 files changed, 163 insertions(+) create mode 100644 benchmarks/algorithm/OrientationIndexStressTest.cpp diff --git a/benchmarks/algorithm/CMakeLists.txt b/benchmarks/algorithm/CMakeLists.txt index cf89a52e54..3b8f867a52 100644 --- a/benchmarks/algorithm/CMakeLists.txt +++ b/benchmarks/algorithm/CMakeLists.txt @@ -18,6 +18,9 @@ target_link_libraries(perf_unaryunion_segments geos) if (BUILD_BENCHMARKS) +add_executable(stress_orientation OrientationIndexStressTest.cpp) +target_link_libraries(stress_orientation geos) + if (benchmark_FOUND) add_executable(perf_orientation OrientationIndexPerfTest.cpp ${PROJECT_SOURCE_DIR}/src/algorithm/CGAlgorithmsDD.cpp diff --git a/benchmarks/algorithm/OrientationIndexStressTest.cpp b/benchmarks/algorithm/OrientationIndexStressTest.cpp new file mode 100644 index 0000000000..28a9e1ffbe --- /dev/null +++ b/benchmarks/algorithm/OrientationIndexStressTest.cpp @@ -0,0 +1,160 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2025 Martin Davis + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU Lesser General Public Licence as published + * by the Free Software Foundation. + * See the COPYING file for more information. + * + **********************************************************************/ + +/*-------------------------------------------------------------- +Stress test for the Orientation Index implementation. + +Usage: stress_orientation [ -v ] + +A robust orientation index implementation should be consistent +- i.e. it should produce the same result for the 3 possible +permutations of the input coordinates which have the same orientation: + +p0-p1 / p2 p1-p2 / p0 p2-p0 / p1 + +The robust implementation uses DoubleDouble arithmetic and a filter to improve computation time. +It is compared to +the simple floating-point orientation computation, which is not consistent. +--------------------------------------------------------------*/ + +#include +#include +#include +#include + +using namespace geos::algorithm; +using namespace geos::geom; +using namespace geos::io; + +int +orientationIndexFP(const Coordinate& p1, const Coordinate& p2,const Coordinate& q){ + double dx1 = p2.x-p1.x; + double dy1 = p2.y-p1.y; + double dx2 = q.x-p2.x; + double dy2 = q.y-p2.y; + double det = dx1 * dy2 - dx2 * dy1; + if (det > 0.0) return 1; + if (det < 0.0) return-1; + return 0; +} + +bool isVerbose = false; + +std::size_t failDD; +std::size_t failFP; + +void parseFlag(char* arg) { + char flag = arg[1]; + switch (flag) { + case 'v': + isVerbose = true; break; + } +} + +void parseArgs(int argc, char** argv) { + if (argc <= 1) return; + int i = 1; + //-- parse flags + while (i < argc && argv[i][0] == '-' && isalpha(argv[i][1])) { + parseFlag(argv[i]); + i++; + } +} + +char orientSym(int orientationIndex) { + if (orientationIndex < 0) return '-'; + if (orientationIndex > 0) return '+'; + return '0'; +} + +void report(std::string name, int orient0, int orient1, + int orient2) { + + std::string consistentInd = (orient0 == orient1 && orient0 == orient2) ? " " : " 1) { + parseArgs(argc, argv); + } + + srand (static_cast (time(0))); + + int i = 0; + while (true) { + runTest(); + i++; + + if (i % 1000 == 0) { + std::cout << "Num tests: " << i + << " DD fail = " << failDD << " (" << (100.0 * failDD / (double) i) << "%)" + << " FP fail = " << failFP << " (" << (100.0 * failFP / (double) i) << "%)" + << std::endl; + } + } +} \ No newline at end of file From ee0238b9a3ae4303b8e6d8985762fe3ee11eab1b Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 10 Feb 2025 14:16:15 -0800 Subject: [PATCH 2/5] Improve consistency testing --- .../algorithm/OrientationIndexStressTest.cpp | 73 +++++++++++-------- 1 file changed, 43 insertions(+), 30 deletions(-) diff --git a/benchmarks/algorithm/OrientationIndexStressTest.cpp b/benchmarks/algorithm/OrientationIndexStressTest.cpp index 28a9e1ffbe..a36fe7345a 100644 --- a/benchmarks/algorithm/OrientationIndexStressTest.cpp +++ b/benchmarks/algorithm/OrientationIndexStressTest.cpp @@ -16,16 +16,19 @@ Stress test for the Orientation Index implementation. Usage: stress_orientation [ -v ] + -v - displays the input and results from each test -A robust orientation index implementation should be consistent +A robust orientation index implementation should be internally consistent - i.e. it should produce the same result for the 3 possible permutations of the input coordinates which have the same orientation: p0-p1 / p2 p1-p2 / p0 p2-p0 / p1 +Also, the reverse orientations should themselves be consistent, +and be opposite in sign to the forward orientation. + The robust implementation uses DoubleDouble arithmetic and a filter to improve computation time. -It is compared to -the simple floating-point orientation computation, which is not consistent. +It is compared to the simple Floating-Point orientation computation, which is not robust. --------------------------------------------------------------*/ #include @@ -78,39 +81,49 @@ char orientSym(int orientationIndex) { return '0'; } -void report(std::string name, int orient0, int orient1, - int orient2) { - - std::string consistentInd = (orient0 == orient1 && orient0 == orient2) ? " " : " orientFunc ) { - int orient0 = Orientation::index(p0, p1, p2); - int orient1 = Orientation::index(p1, p2, p0); - int orient2 = Orientation::index(p2, p0, p1); + int orient0 = orientFunc(p0, p1, p2); + int orient1 = orientFunc(p1, p2, p0); + int orient2 = orientFunc(p2, p0, p1); + bool isConsistentForward = orient0 == orient1 && orient0 == orient2; + + int orientRev0 = orientFunc(p1, p0, p2); + int orientRev1 = orientFunc(p0, p2, p1); + int orientRev2 = orientFunc(p2, p1, p0); + bool isConsistentRev = orientRev0 == orientRev1 && orientRev0 == orientRev2; + + bool isConsistent = isConsistentForward && isConsistentRev; + if (isConsistent) { + bool isOpposite = orient0 == -orientRev0; + isConsistent &= isOpposite; + } if (isVerbose) { - report("DD", orient0, orient1, orient2); + std::string consistentInd = isConsistent ? " " : " int { + return Orientation::index(p0, p1, p2); + }); } -bool isAllOrientationsEqualFP(Coordinate p0, Coordinate p1, Coordinate p2) +bool isConsistentFP(Coordinate p0, Coordinate p1, Coordinate p2) { - int orient0 = orientationIndexFP(p0, p1, p2); - int orient1 = orientationIndexFP(p1, p2, p0); - int orient2 = orientationIndexFP(p2, p0, p1); - if (isVerbose) { - report("FP", orient0, orient1, orient2); - } - return orient0 == orient1 && orient0 == orient2; + return isConsistent("FP", p0, p1, p2, + [](Coordinate p0, Coordinate p1, Coordinate p2) -> int { + return orientationIndexFP(p0, p1, p2); + }); } Coordinate randomCoord() { @@ -126,8 +139,8 @@ void runTest() Coordinate p2 = Coordinate(LineSegment::midPoint(p0, p1)); - bool isCorrectDD = isAllOrientationsEqualDD(p0, p1, p2); - bool isCorrectFP = isAllOrientationsEqualFP(p0, p1, p2); + bool isCorrectDD = isConsistentDD(p0, p1, p2); + bool isCorrectFP = isConsistentFP(p0, p1, p2); if (isVerbose) { std::cout << " " << WKTWriter::toLineString(p0, p1) << " - " << WKTWriter::toPoint(p2) From 04ac61167bef386a9a33381a78f86e55cf89162f Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 10 Feb 2025 14:19:09 -0800 Subject: [PATCH 3/5] Fix incorrect test tag --- benchmarks/algorithm/OrientationIndexStressTest.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/algorithm/OrientationIndexStressTest.cpp b/benchmarks/algorithm/OrientationIndexStressTest.cpp index a36fe7345a..d37b7b6fc8 100644 --- a/benchmarks/algorithm/OrientationIndexStressTest.cpp +++ b/benchmarks/algorithm/OrientationIndexStressTest.cpp @@ -112,7 +112,7 @@ bool isConsistent(std::string tag, Coordinate p0, Coordinate p1, Coordinate p2, bool isConsistentDD(Coordinate p0, Coordinate p1, Coordinate p2) { - return isConsistent("FP", p0, p1, p2, + return isConsistent("DD", p0, p1, p2, [](Coordinate p0, Coordinate p1, Coordinate p2) -> int { return Orientation::index(p0, p1, p2); }); From 36351eb0ab03016fa8deee6e95e08e764b7bbb0e Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 10 Feb 2025 14:44:29 -0800 Subject: [PATCH 4/5] Improve input display precision --- benchmarks/algorithm/OrientationIndexStressTest.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/benchmarks/algorithm/OrientationIndexStressTest.cpp b/benchmarks/algorithm/OrientationIndexStressTest.cpp index d37b7b6fc8..41337b6dd9 100644 --- a/benchmarks/algorithm/OrientationIndexStressTest.cpp +++ b/benchmarks/algorithm/OrientationIndexStressTest.cpp @@ -36,6 +36,8 @@ It is compared to the simple Floating-Point orientation computation, which is no #include #include +#include + using namespace geos::algorithm; using namespace geos::geom; using namespace geos::io; @@ -143,7 +145,9 @@ void runTest() bool isCorrectFP = isConsistentFP(p0, p1, p2); if (isVerbose) { - std::cout << " " << WKTWriter::toLineString(p0, p1) << " - " << WKTWriter::toPoint(p2) + std::cout << std::setprecision(20) << " " + << "LINESTRING ( " << p0.x << " " << p0.y << ", " << p1.x << " " << p1.y << " )" + << " - " << "POINT ( " << p2.x << " " << p2.y << " )" << std::endl; } From 3d2dcbcdd8584691f3701d7e040ec2f2b9ad67e5 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 10 Feb 2025 16:42:48 -0800 Subject: [PATCH 5/5] Add diagonal line segment tests --- .../algorithm/OrientationIndexStressTest.cpp | 82 ++++++++++++++----- 1 file changed, 63 insertions(+), 19 deletions(-) diff --git a/benchmarks/algorithm/OrientationIndexStressTest.cpp b/benchmarks/algorithm/OrientationIndexStressTest.cpp index 41337b6dd9..af9e6be667 100644 --- a/benchmarks/algorithm/OrientationIndexStressTest.cpp +++ b/benchmarks/algorithm/OrientationIndexStressTest.cpp @@ -15,7 +15,8 @@ /*-------------------------------------------------------------- Stress test for the Orientation Index implementation. -Usage: stress_orientation [ -v ] +Usage: stress_orientation [ -v ] [ -d ] + -d - run diagonal line segment tests -v - displays the input and results from each test A robust orientation index implementation should be internally consistent @@ -29,6 +30,10 @@ and be opposite in sign to the forward orientation. The robust implementation uses DoubleDouble arithmetic and a filter to improve computation time. It is compared to the simple Floating-Point orientation computation, which is not robust. + +Two kinds of test generators are provided: +- random line segments with midpoints +- points at increasing decimal intervals along a diagonal line segment LINESTRING(2 0, 0 2) --------------------------------------------------------------*/ #include @@ -55,7 +60,9 @@ orientationIndexFP(const Coordinate& p1, const Coordinate& p2,const Coordinate& } bool isVerbose = false; +bool isDiagonalTest = false; +std::size_t count; std::size_t failDD; std::size_t failFP; @@ -64,6 +71,8 @@ void parseFlag(char* arg) { switch (flag) { case 'v': isVerbose = true; break; + case 'd': + isDiagonalTest = true; break; } } @@ -134,13 +143,8 @@ Coordinate randomCoord() { return Coordinate(x, y); } -void runTest() +void checkTest(Coordinate p0, Coordinate p1, Coordinate p2) { - Coordinate p0 = randomCoord(); - Coordinate p1 = randomCoord(); - - Coordinate p2 = Coordinate(LineSegment::midPoint(p0, p1)); - bool isCorrectDD = isConsistentDD(p0, p1, p2); bool isCorrectFP = isConsistentFP(p0, p1, p2); @@ -151,27 +155,67 @@ void runTest() << std::endl; } + count++; if (! isCorrectDD) failDD++; if (! isCorrectFP) failFP++; } -int main(int argc, char** argv) { - if (argc > 1) { - parseArgs(argc, argv); +void reportStats(std::string tag = "") { + std::cout << tag << "Num tests: " << count + << " DD fail = " << failDD << " (" << round((100.0 * failDD / (double) count)) << "%)" + << " FP fail = " << failFP << " (" << round((100.0 * failFP / (double) count)) << "%)" + << std::endl; +} + +void runDiagonalTests() +{ + const int DIAG_SIZE = 2; + Coordinate p0 = Coordinate(DIAG_SIZE, 0); + Coordinate p1 = Coordinate(0, DIAG_SIZE); + + const int MAX_PRECISION = 3; + int d = 10; + for (int i = 0; i < MAX_PRECISION; i++) { + int num_points = DIAG_SIZE * d; + for (int ix = 0; ix <= num_points; ix++) { + int iy = num_points - ix; + double x = ix / (double) d; + double y = iy / (double) d; + Coordinate p2 = Coordinate(x, y); + checkTest(p0, p1, p2); + } + d *= 10; + reportStats(); } +} +void runRandomTests() +{ srand (static_cast (time(0))); - int i = 0; - while (true) { - runTest(); - i++; + const size_t MAX_ITER = 10'000'000; + + for (size_t i = 1 ; i <= MAX_ITER; i++) { + Coordinate p0 = randomCoord(); + Coordinate p1 = randomCoord(); + Coordinate p2 = Coordinate(LineSegment::midPoint(p0, p1)); + checkTest(p0, p1, p2); - if (i % 1000 == 0) { - std::cout << "Num tests: " << i - << " DD fail = " << failDD << " (" << (100.0 * failDD / (double) i) << "%)" - << " FP fail = " << failFP << " (" << (100.0 * failFP / (double) i) << "%)" - << std::endl; + if (i % 10'000 == 0) { + reportStats(); } } +} + +int main(int argc, char** argv) { + if (argc > 1) { + parseArgs(argc, argv); + } + if (isDiagonalTest) { + runDiagonalTests(); + } + else { + runRandomTests(); + } + reportStats("Final: "); } \ No newline at end of file