From 928379cec5321855125b470061b8b30e6f40a0e0 Mon Sep 17 00:00:00 2001 From: angusj Date: Wed, 13 Dec 2023 15:15:27 +1000 Subject: [PATCH] Added CLIPPER2_HI_PRECISION option in CMakeLists.txt --- CPP/BenchMark/GetIntersectPtBenchmark.cpp | 83 +++++++++---------- CPP/CMakeLists.txt | 26 ++++-- .../include/clipper2/clipper.core.h | 70 +++++++++++++--- .../include/clipper2/clipper.engine.h | 4 +- CPP/Clipper2Lib/include/clipper2/clipper.h | 16 ++-- 5 files changed, 124 insertions(+), 75 deletions(-) diff --git a/CPP/BenchMark/GetIntersectPtBenchmark.cpp b/CPP/BenchMark/GetIntersectPtBenchmark.cpp index ea9697d4..3d4e3852 100644 --- a/CPP/BenchMark/GetIntersectPtBenchmark.cpp +++ b/CPP/BenchMark/GetIntersectPtBenchmark.cpp @@ -38,6 +38,40 @@ struct SetConsoleTextColor ////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////// +// Miscellaneous functions +///////////////////////////////////////////////////////// + +double GetSineAbc(const Point64& a, const Point64& b, const Point64& c) +{ + double dpB = DotProduct(a, b, c); + double SqrCosB = dpB * dpB / (DistanceSqr(a, b) * DistanceSqr(b, c)); + double cos2B = SqrCosB * 2 - 1; // trig. itentity + return std::sqrt(1 - cos2B); // sin(B) = Sqrt(1-cos(2B)) +} + +static inline Point64 ReflectPoint(const Point64& pt, const Point64& pivot) +{ + return Point64(pivot.x + (pivot.x - pt.x), pivot.y + (pivot.y - pt.y)); +} + +static inline Point64 MidPoint(const Point64& p1, const Point64& p2) +{ + Point64 result; + result.x = (p1.x + p2.x) / 2; + result.y = (p1.y + p2.y) / 2; + return result; +} + +static inline Point64 MakeRandomPoint(int64_t min_val, int64_t max_val) +{ + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution x(min_val, max_val); + std::uniform_int_distribution y(min_val, max_val); + return Point64(x(gen), y(gen)); +} + typedef std::function GipFunction; @@ -106,7 +140,9 @@ static bool GIP_Func_F(const Point64& ln1a, const Point64& ln1b, ///////////////////////////////////////////////////////// // GIP_F_Mod: Modified GIP_Func_F (see above). -// Uses static_cast instead of nearbyint. +// Uses static_cast instead of nearbyint. Surprisingly, +// while faster here, this is much slower than GIP_Func_F +// when replacing GetIntersectPoint() in clipper.core.h. ///////////////////////////////////////////////////////// #define CC_MIN(x,y) ((x)>(y)?(y):(x)) #define CC_MAX(x,y) ((x)<(y)?(y):(x)) @@ -160,46 +196,6 @@ typedef std::vector::iterator test_iter; // global data std::vector tests; -///////////////////////////////////////////////////////// -// Miscellaneous functions -///////////////////////////////////////////////////////// - -double GetSine(const Point64& a, const Point64& b, const Point64& c) -{ - Point64 ab = Point64(b.x - a.x, b.y - a.y); - Point64 bc = Point64(b.x - c.x, b.y - c.y); - double dp = DotProduct(ab,bc); - double distSqr1 = DistanceSqr(a, b); - double distSqr2 = DistanceSqr(b, c); - // cos(A)^2 = dp^2/DistSqr(ab)/DistSqr(bc) - double cosSqr = dp * dp / distSqr1 / distSqr2; - // cos2A = cos(A)^2 * 2 - 1 - double cos2 = cosSqr * 2 - 1; - // sinA = Sqrt(1-cos2A) - return std::sqrt(1 - cos2); -} - -static inline Point64 ReflectPoint(const Point64& pt, const Point64& pivot) -{ - return Point64(pivot.x + (pivot.x - pt.x), pivot.y + (pivot.y - pt.y)); -} - -static inline Point64 MidPoint(const Point64& p1, const Point64& p2) -{ - Point64 result; - result.x = (p1.x + p2.x) / 2; - result.y = (p1.y + p2.y) / 2; - return result; -} - -static inline Point64 MakeRandomPoint(int64_t min_val, int64_t max_val) -{ - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_int_distribution x(min_val, max_val); - std::uniform_int_distribution y(min_val, max_val); - return Point64(x(gen), y(gen)); -} static inline GipFunction GetGipFunc(int index) { @@ -285,7 +281,6 @@ int main(int argc, char** argv) BENCHMARK(BM_GIP_Current); BENCHMARK(BM_GIP_Func_F); BENCHMARK(BM_GIP_F_Mod); - benchmark::Counter(); bool first_pass = true; for (int current_pow10 = 12; current_pow10 <= 18; ++current_pow10) @@ -294,7 +289,7 @@ int main(int argc, char** argv) // at their midpoints, while using random coordinates that are // restricted to the specified power of 10 range int64_t max_coord = static_cast(pow(10, current_pow10)); - for (int64_t i = 0; i < 10000; ++i) + for (int64_t i = 0; i < 100000; ++i) { Point64 ip1 = MakeRandomPoint(-max_coord, max_coord); Point64 ip2 = MakeRandomPoint(-max_coord, max_coord); @@ -306,7 +301,7 @@ int main(int argc, char** argv) // the closer segments are to collinear, the less accurate are // calculations that determine intersection points. // So exclude segments that are **almost** collinear. eg. sin(1deg) ~= 0.017 - if (std::abs(GetSine(ip1, actual, ip3)) < 0.017) continue; + if (std::abs(GetSineAbc(ip1, actual, ip3)) < 0.017) continue; // alternatively, only exclude segments that are collinear //if (!CrossProduct(ip1, actual, ip3)) continue; tests.push_back(TestRecord(participants, actual, ip1, ip2, ip3, ip4)); diff --git a/CPP/CMakeLists.txt b/CPP/CMakeLists.txt index 5fab727d..f95dd295 100644 --- a/CPP/CMakeLists.txt +++ b/CPP/CMakeLists.txt @@ -9,21 +9,27 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) set_property(GLOBAL PROPERTY USE_FOLDERS ON) +# CLIPPER2_HI_PRECISION: See GetIntersectPoint() in clipper.core.h +option(CLIPPER2_HI_PRECISION "Caution: enabling this will compromise performance" OFF) + option(CLIPPER2_UTILS "Build utilities" ON) option(CLIPPER2_EXAMPLES "Build examples" ON) option(CLIPPER2_TESTS "Build tests" ON) option(USE_EXTERNAL_GTEST "Use system-wide installed GoogleTest" OFF) -option(USE_EXTERNAL_GBENCHMARK "Use the googlebenchmark" OFF) +option(USE_EXTERNAL_GBENCHMARK "Use the googlebenchmark" OFF) option(BUILD_SHARED_LIBS "Build shared libs" OFF) set(CLIPPER2_USINGZ "ON" CACHE STRING "Build Clipper2Z, either \"ON\" or \"OFF\" or \"ONLY\"") -set(CLIPPER2_MAX_PRECISION 8 CACHE STRING "Maximum precision allowed for double to int64 scaling") + +# CLIPPER2_MAX_DECIMAL_PRECISION: maximum decimal precision when scaling PathsD to Paths64. +# Caution: excessive scaling will increase the likelihood of integer overflow errors. +set(CLIPPER2_MAX_DECIMAL_PRECISION 8 CACHE STRING "Maximum decimal precision range") + if (APPLE) set(CMAKE_SHARED_LIBRARY_SUFFIX ".dylib") endif () include(GNUInstallDirs) - set(CLIPPER2_INC_FOLDER ${PROJECT_SOURCE_DIR}/Clipper2Lib/include/clipper2) configure_file(clipper.version.in ${CLIPPER2_INC_FOLDER}/clipper.version.h) @@ -52,10 +58,12 @@ if (NOT (CLIPPER2_USINGZ STREQUAL "ONLY")) add_library(Clipper2 ${CLIPPER2_INC} ${CLIPPER2_SRC}) target_compile_definitions( - Clipper2 - PUBLIC CLIPPER2_MAX_PRECISION=${CLIPPER2_MAX_PRECISION} + Clipper2 PUBLIC + CLIPPER2_MAX_DECIMAL_PRECISION=${CLIPPER2_MAX_DECIMAL_PRECISION} + $<$:CLIPPER2_HI_PRECISION> ) + target_include_directories(Clipper2 PUBLIC Clipper2Lib/include ) @@ -74,9 +82,11 @@ if (NOT (CLIPPER2_USINGZ STREQUAL "OFF")) add_library(Clipper2Z ${CLIPPER2_INC} ${CLIPPER2_SRC}) target_compile_definitions( - Clipper2Z - PUBLIC USINGZ CLIPPER2_MAX_PRECISION=${CLIPPER2_MAX_PRECISION} - ) + Clipper2Z PUBLIC + USINGZ + CLIPPER2_MAX_DECIMAL_PRECISION=${CLIPPER2_MAX_DECIMAL_PRECISION} + $<$:CLIPPER2_HI_PRECISION> + ) target_include_directories(Clipper2Z PUBLIC Clipper2Lib/include diff --git a/CPP/Clipper2Lib/include/clipper2/clipper.core.h b/CPP/Clipper2Lib/include/clipper2/clipper.core.h index 141e522e..69910789 100644 --- a/CPP/Clipper2Lib/include/clipper2/clipper.core.h +++ b/CPP/Clipper2Lib/include/clipper2/clipper.core.h @@ -1,6 +1,6 @@ /******************************************************************************* * Author : Angus Johnson * -* Date : 9 December 2023 * +* Date : 13 December 2023 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2023 * * Purpose : Core Clipper Library structures and functions * @@ -58,10 +58,10 @@ namespace Clipper2Lib static const double PI = 3.141592653589793238; #endif -#ifdef CLIPPER2_MAX_PRECISION - const int MAX_DECIMAL_PRECISION = CLIPPER2_MAX_PRECISION; +#ifdef CLIPPER2_MAX_DECIMAL_PRECISION + const int CLIPPER2_MAX_DEC_PRECISION = CLIPPER2_MAX_DECIMAL_PRECISION; #else - const int MAX_DECIMAL_PRECISION = 8; // see Discussions #564 + const int CLIPPER2_MAX_DEC_PRECISION = 8; // see Discussions #564 #endif static const int64_t MAX_COORD = INT64_MAX >> 2; @@ -242,6 +242,14 @@ namespace Clipper2Lib (std::numeric_limits::max)(), (std::numeric_limits::max)()); + template + static inline Point MidPoint(const Point& p1, const Point& p2) + { + Point result; + result.x = (p1.x + p2.x) / 2; + result.y = (p1.y + p2.y) / 2; + return result; + } // Rect ------------------------------------------------------------------------ @@ -623,18 +631,19 @@ namespace Clipper2Lib // Miscellaneous ------------------------------------------------------------ - inline void CheckPrecision(int& precision, int& error_code) + inline void CheckPrecisionRange(int& precision, int& error_code) { - if (precision >= -MAX_DECIMAL_PRECISION && precision <= MAX_DECIMAL_PRECISION) return; + if (precision >= -CLIPPER2_MAX_DEC_PRECISION && + precision <= CLIPPER2_MAX_DEC_PRECISION) return; error_code |= precision_error_i; // non-fatal error - DoError(precision_error_i); // does nothing unless exceptions enabled - precision = precision > 0 ? MAX_DECIMAL_PRECISION : -MAX_DECIMAL_PRECISION; + DoError(precision_error_i); // does nothing when exceptions are disabled + precision = precision > 0 ? CLIPPER2_MAX_DEC_PRECISION : -CLIPPER2_MAX_DEC_PRECISION; } - inline void CheckPrecision(int& precision) + inline void CheckPrecisionRange(int& precision) { int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); } template @@ -721,6 +730,40 @@ namespace Clipper2Lib return Area(poly) >= 0; } +#if CLIPPER2_HI_PRECISION + // caution: this will compromise performance + // https://github.com/AngusJohnson/Clipper2/issues/317#issuecomment-1314023253 + // See also CPP/BenchMark/GetIntersectPtBenchmark.cpp + #define CC_MIN(x,y) ((x)>(y)?(y):(x)) + #define CC_MAX(x,y) ((x)<(y)?(y):(x)) + inline bool GetIntersectPoint(const Point64& ln1a, const Point64& ln1b, + const Point64& ln2a, const Point64& ln2b, Point64& ip) + { + double ln1dy = (double)(ln1b.y - ln1a.y); + double ln1dx = (double)(ln1a.x - ln1b.x); + double ln2dy = (double)(ln2b.y - ln2a.y); + double ln2dx = (double)(ln2a.x - ln2b.x); + double det = (ln2dy * ln1dx) - (ln1dy * ln2dx); + if (det == 0.0) return false; + int64_t bb0minx = CC_MIN(ln1a.x, ln1b.x); + int64_t bb0miny = CC_MIN(ln1a.y, ln1b.y); + int64_t bb0maxx = CC_MAX(ln1a.x, ln1b.x); + int64_t bb0maxy = CC_MAX(ln1a.y, ln1b.y); + int64_t bb1minx = CC_MIN(ln2a.x, ln2b.x); + int64_t bb1miny = CC_MIN(ln2a.y, ln2b.y); + int64_t bb1maxx = CC_MAX(ln2a.x, ln2b.x); + int64_t bb1maxy = CC_MAX(ln2a.y, ln2b.y); + int64_t originx = (CC_MIN(bb0maxx, bb1maxx) + CC_MAX(bb0minx, bb1minx)) >> 1; + int64_t originy = (CC_MIN(bb0maxy, bb1maxy) + CC_MAX(bb0miny, bb1miny)) >> 1; + double ln0c = (ln1dy * (double)(ln1a.x - originx)) + (ln1dx * (double)(ln1a.y - originy)); + double ln1c = (ln2dy * (double)(ln2a.x - originx)) + (ln2dx * (double)(ln2a.y - originy)); + double hitx = ((ln1dx * ln1c) - (ln2dx * ln0c)) / det; + double hity = ((ln2dy * ln0c) - (ln1dy * ln1c)) / det; + ip.x = originx + (int64_t)nearbyint(hitx); + ip.y = originy + (int64_t)nearbyint(hity); + return true; +} +#else inline bool GetIntersectPoint(const Point64& ln1a, const Point64& ln1b, const Point64& ln2a, const Point64& ln2b, Point64& ip) { @@ -733,15 +776,16 @@ namespace Clipper2Lib double det = dy1 * dx2 - dy2 * dx1; if (det == 0.0) return false; double t = ((ln1a.x - ln2a.x) * dy2 - (ln1a.y - ln2a.y) * dx2) / det; - if (t <= 0.0) ip = ln1a; // ?? check further (see also #568) - else if (t >= 1.0) ip = ln1b; // ?? check further + if (t <= 0.0) ip = ln1a; + else if (t >= 1.0) ip = ln1b; else { ip.x = static_cast(ln1a.x + t * dx1); ip.y = static_cast(ln1a.y + t * dy1); - } + } return true; } +#endif inline bool SegmentsIntersect(const Point64& seg1a, const Point64& seg1b, const Point64& seg2a, const Point64& seg2b, bool inclusive = false) diff --git a/CPP/Clipper2Lib/include/clipper2/clipper.engine.h b/CPP/Clipper2Lib/include/clipper2/clipper.engine.h index f7fe0182..9dd55db7 100644 --- a/CPP/Clipper2Lib/include/clipper2/clipper.engine.h +++ b/CPP/Clipper2Lib/include/clipper2/clipper.engine.h @@ -1,6 +1,6 @@ /******************************************************************************* * Author : Angus Johnson * -* Date : 22 November 2023 * +* Date : 13 December 2023 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2023 * * Purpose : This is the main polygon clipping module * @@ -530,7 +530,7 @@ namespace Clipper2Lib { public: explicit ClipperD(int precision = 2) : ClipperBase() { - CheckPrecision(precision, error_code_); + CheckPrecisionRange(precision, error_code_); // to optimize scaling / descaling precision // set the scale to a power of double's radix (2) (#25) scale_ = std::pow(std::numeric_limits::radix, diff --git a/CPP/Clipper2Lib/include/clipper2/clipper.h b/CPP/Clipper2Lib/include/clipper2/clipper.h index 33beef6f..f9446220 100644 --- a/CPP/Clipper2Lib/include/clipper2/clipper.h +++ b/CPP/Clipper2Lib/include/clipper2/clipper.h @@ -1,6 +1,6 @@ /******************************************************************************* * Author : Angus Johnson * -* Date : 18 November 2023 * +* Date : 13 December 2023 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2023 * * Purpose : This module provides a simple interface to the Clipper Library * @@ -47,7 +47,7 @@ namespace Clipper2Lib { const PathsD& subjects, const PathsD& clips, int precision = 2) { int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); PathsD result; if (error_code) return result; ClipperD clipper(precision); @@ -63,7 +63,7 @@ namespace Clipper2Lib { { polytree.Clear(); int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (error_code) return; ClipperD clipper(precision); clipper.AddSubject(subjects); @@ -104,7 +104,7 @@ namespace Clipper2Lib { { PathsD result; int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (error_code) return result; ClipperD clipper(precision); clipper.AddSubject(subjects); @@ -149,7 +149,7 @@ namespace Clipper2Lib { int precision = 2, double arc_tolerance = 0.0) { int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (!delta) return paths; if (error_code) return PathsD(); const double scale = std::pow(10, precision); @@ -219,7 +219,7 @@ namespace Clipper2Lib { { if (rect.IsEmpty() || paths.empty()) return PathsD(); int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (error_code) return PathsD(); const double scale = std::pow(10, precision); Rect64 r = ScaleRect(rect, scale); @@ -251,7 +251,7 @@ namespace Clipper2Lib { { if (rect.IsEmpty() || lines.empty()) return PathsD(); int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (error_code) return PathsD(); const double scale = std::pow(10, precision); Rect64 r = ScaleRect(rect, scale); @@ -545,7 +545,7 @@ namespace Clipper2Lib { inline PathD TrimCollinear(const PathD& path, int precision, bool is_open_path = false) { int error_code = 0; - CheckPrecision(precision, error_code); + CheckPrecisionRange(precision, error_code); if (error_code) return PathD(); const double scale = std::pow(10, precision); Path64 p = ScalePath(path, scale, error_code);