QA: add a test for CalcArcCenter( aStart, aMid, aEnd )
This commit is contained in:
Alex Shvartzkop
parent
df83e24eb7
commit
b7824adfb1
|
|
@ -31,12 +31,243 @@
|
|||
#include <trigo.h>
|
||||
|
||||
|
||||
/*
|
||||
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<double>::epsilon();
|
||||
|
||||
if( xDelta_32 == 0.0 )
|
||||
xDelta_32 = -std::numeric_limits<double>::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<double>::epsilon();
|
||||
bSlope -= std::numeric_limits<double>::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<TEST_CALC_ARC_CENTER_CASE> 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 ) );
|
||||
|
|
|
|||
Loading…
Reference in New Issue