From b7824adfb180e8ff9a9719dc23cf5da950fd3982 Mon Sep 17 00:00:00 2001 From: Alex Shvartzkop Date: Sat, 25 Nov 2023 18:37:48 +0300 Subject: [PATCH] QA: add a test for CalcArcCenter( aStart, aMid, aEnd ) --- qa/tests/libs/kimath/test_kimath.cpp | 231 +++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) diff --git a/qa/tests/libs/kimath/test_kimath.cpp b/qa/tests/libs/kimath/test_kimath.cpp index da8757d215..0f3424e4f6 100644 --- a/qa/tests/libs/kimath/test_kimath.cpp +++ b/qa/tests/libs/kimath/test_kimath.cpp @@ -31,12 +31,243 @@ #include +/* +CircleCenterFrom3Points calculate the center of a circle defined by 3 points +It is similar to CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, const VECTOR2D& aEnd ) +but it was needed to debug CalcArcCenter, so I keep it available for other issues in CalcArcCenter + +The perpendicular bisector of the segment between two points is the +set of all points equidistant from both. So if you take the +perpendicular bisector of (x1,y1) and (x2,y2) and the perpendicular +bisector of the segment from (x2,y2) to (x3,y3) and find the +intersection of those lines, that point will be the center. + +To find the equation of the perpendicular bisector of (x1,y1) to (x2,y2), +you know that it passes through the midpoint of the segment: +((x1+x2)/2,(y1+y2)/2), and if the slope of the line +connecting (x1,y1) to (x2,y2) is m, the slope of the perpendicular +bisector is -1/m. Work out the equations for the two lines, find +their intersection, and bingo! You've got the coordinates of the center. + +An error should occur if the three points lie on a line, and you'll +need special code to check for the case where one of the slopes is zero. + +see https://web.archive.org/web/20171223103555/http://mathforum.org/library/drmath/view/54323.html +*/ +static bool Ref0CircleCenterFrom3Points( const VECTOR2D& p1, const VECTOR2D& p2, const VECTOR2D& p3, + VECTOR2D* aCenter ) +{ + // Move coordinate origin to p2, to simplify calculations + VECTOR2D b = p1 - p2; + VECTOR2D d = p3 - p2; + double bc = ( b.x * b.x + b.y * b.y ) / 2.0; + double cd = ( -d.x * d.x - d.y * d.y ) / 2.0; + double det = -b.x * d.y + d.x * b.y; + + if( fabs( det ) < 1.0e-6 ) // arbitrary limit to avoid divide by 0 + return false; + + det = 1 / det; + aCenter->x = ( -bc * d.y - cd * b.y ) * det; + aCenter->y = ( b.x * cd + d.x * bc ) * det; + *aCenter += p2; + + return true; +} + + +// Based on https://paulbourke.net/geometry/circlesphere/ +// Uses the Y'b equation too, depending on slope values +static const VECTOR2D Ref1CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, const VECTOR2D& aEnd ) +{ + VECTOR2D center; + double yDelta_21 = aMid.y - aStart.y; + double xDelta_21 = aMid.x - aStart.x; + double yDelta_32 = aEnd.y - aMid.y; + double xDelta_32 = aEnd.x - aMid.x; + + // This is a special case for aMid as the half-way point when aSlope = 0 and bSlope = inf + // or the other way around. In that case, the center lies in a straight line between + // aStart and aEnd + if( ( ( xDelta_21 == 0.0 ) && ( yDelta_32 == 0.0 ) ) + || ( ( yDelta_21 == 0.0 ) && ( xDelta_32 == 0.0 ) ) ) + { + center.x = ( aStart.x + aEnd.x ) / 2.0; + center.y = ( aStart.y + aEnd.y ) / 2.0; + return center; + } + + // Prevent div=0 errors + /*if( xDelta_21 == 0.0 ) + xDelta_21 = std::numeric_limits::epsilon(); + + if( xDelta_32 == 0.0 ) + xDelta_32 = -std::numeric_limits::epsilon();*/ + + double aSlope = yDelta_21 / xDelta_21; + double bSlope = yDelta_32 / xDelta_32; + + if( aSlope == bSlope ) + { + if( aStart == aEnd ) + { + // This is a special case for a 360 degrees arc. In this case, the center is halfway between + // the midpoint and either end point + center.x = ( aStart.x + aMid.x ) / 2.0; + center.y = ( aStart.y + aMid.y ) / 2.0; + return center; + } + else + { + // If the points are colinear, the center is at infinity, so offset + // the slope by a minimal amount + // Warning: This will induce a small error in the center location + /*aSlope += std::numeric_limits::epsilon(); + bSlope -= std::numeric_limits::epsilon();*/ + } + } + + center.x = ( aSlope * bSlope * ( aStart.y - aEnd.y ) + bSlope * ( aStart.x + aMid.x ) + - aSlope * ( aMid.x + aEnd.x ) ) + / ( 2 * ( bSlope - aSlope ) ); + + if( std::abs( aSlope ) > std::abs( bSlope ) ) + { + center.y = ( ( ( aStart.x + aMid.x ) / 2.0 - center.x ) / aSlope + + ( aStart.y + aMid.y ) / 2.0 ); + } + else + { + center.y = + ( ( ( aMid.x + aEnd.x ) / 2.0 - center.x ) / bSlope + ( aMid.y + aEnd.y ) / 2.0 ); + } + + return center; +} + + +struct TEST_CALC_ARC_CENTER_CASE +{ + VECTOR2I istart; + VECTOR2I imid; + VECTOR2I iend; +}; + + /** * Declare the test suite */ BOOST_AUTO_TEST_SUITE( KiMath ) +BOOST_AUTO_TEST_CASE( TestCalcArcCenter3Pts ) +{ + // Order: start, mid, end + std::vector calc_center_cases = { + //{ { 183000000, 89000000 }, { 185333332, 89000000 }, { 185333333, 91333333 } } // Fails currently + }; + + const double tolerance = 1.0; + + // Check random points (fails currently) + /*for( int i = 0; i < 100; i++ ) + { + TEST_CALC_ARC_CENTER_CASE cas; + + cas.istart.x = rand(); + cas.istart.y = rand(); + + cas.imid.x = rand(); + cas.imid.y = rand(); + + cas.iend.x = rand(); + cas.iend.y = rand(); + + calc_center_cases.push_back( cas ); + }*/ + + for( const TEST_CALC_ARC_CENTER_CASE& entry : calc_center_cases ) + { + wxString msg; + + VECTOR2D start( entry.istart ); + VECTOR2D mid( entry.imid ); + VECTOR2D end( entry.iend ); + + VECTOR2D calcCenter = CalcArcCenter( start, mid, end ); + + double crs = ( calcCenter - start ).EuclideanNorm(); + double crm = ( calcCenter - mid ).EuclideanNorm(); + double cre = ( calcCenter - end ).EuclideanNorm(); + + double cavg = ( crs + crm + cre ) / 3.0; + + if( std::abs( crs - cavg ) > tolerance || std::abs( crm - cavg ) > tolerance + || std::abs( cre - cavg ) > tolerance ) + { + msg << "CalcArcCenter failed."; + + msg << "\nstart: " << entry.istart.Format(); + msg << "\nmid: " << entry.imid.Format(); + msg << "\nend: " << entry.iend.Format(); + + { + msg << "\nCalculated center: " << wxString::Format( "%.15f", calcCenter.x ) << "," + << wxString::Format( "%.15f", calcCenter.y ); + + msg << "\n Avg radius: " << wxString::Format( "%.15f", cavg ); + msg << "\nStart radius: " << wxString::Format( "%.15f", crs ); + msg << "\n Mid radius: " << wxString::Format( "%.15f", crm ); + msg << "\n End radius: " << wxString::Format( "%.15f", cre ); + msg << "\n"; + } + + { + VECTOR2D ref0Center; + Ref0CircleCenterFrom3Points( start, mid, end, &ref0Center ); + + double r0_rs = ( ref0Center - start ).EuclideanNorm(); + double r0_rm = ( ref0Center - mid ).EuclideanNorm(); + double r0_rre = ( ref0Center - end ).EuclideanNorm(); + + double r0_ravg = ( r0_rs + r0_rm + r0_rre ) / 3.0; + + msg << "\nReference0 center: " << wxString::Format( "%.15f", ref0Center.x ) << "," + << wxString::Format( "%.15f", ref0Center.y ); + + msg << "\nRef0 Avg radius: " << wxString::Format( "%.15f", r0_ravg ); + msg << "\nRef0 Start radius: " << wxString::Format( "%.15f", r0_rs ); + msg << "\nRef0 Mid radius: " << wxString::Format( "%.15f", r0_rm ); + msg << "\nRef0 End radius: " << wxString::Format( "%.15f", r0_rre ); + msg << "\n"; + } + + { + VECTOR2D ref1Center = Ref1CalcArcCenter( start, mid, end ); + + double r1_rs = ( ref1Center - start ).EuclideanNorm(); + double r1_rm = ( ref1Center - mid ).EuclideanNorm(); + double r1_rre = ( ref1Center - end ).EuclideanNorm(); + + double r1_ravg = ( r1_rs + r1_rm + r1_rre ) / 3.0; + + msg << "\nReference1 center: " << wxString::Format( "%.15f", ref1Center.x ) << "," + << wxString::Format( "%.15f", ref1Center.y ); + + msg << "\nRef1 Avg radius: " << wxString::Format( "%.15f", r1_ravg ); + msg << "\nRef1 Start radius: " << wxString::Format( "%.15f", r1_rs ); + msg << "\nRef1 Mid radius: " << wxString::Format( "%.15f", r1_rm ); + msg << "\nRef1 End radius: " << wxString::Format( "%.15f", r1_rre ); + msg << "\n"; + } + + + BOOST_CHECK_MESSAGE( false, msg ); + } + } +} + + BOOST_AUTO_TEST_CASE( TestInterceptsPositiveX ) { BOOST_CHECK( !InterceptsPositiveX( 10.0, 20.0 ) );