From 5eb265080f6986988e5c8d0c9c100b82cb98454d Mon Sep 17 00:00:00 2001 From: jean-pierre charras Date: Mon, 19 Feb 2024 15:52:37 +0100 Subject: [PATCH] try to cherry pick commit 7fd9226bec58c9e5114595863a6f51820bff746c from master: Fix issue in CalcArcCenter( VECTOR2D& aStart, VECTOR2D& aMid, VECTOR2D& aEnd ) It happens when the segment (aStart, aMid) is horizontal Probably also when the segment (aEnd, aMid) is horizontal Slopes with value 0.0 are set to double:: epsilon(), but it was a too small values generating broken calculations. Now set to 1e-10 (it seems working). Fixes #16089 https://gitlab.com/kicad/code/kicad/-/issues/16089 --- libs/kimath/src/trigo.cpp | 66 +++++++++++++++++-- .../import_gfx/graphics_importer_pcbnew.cpp | 4 +- pcbnew/pcb_painter.cpp | 2 +- 3 files changed, 64 insertions(+), 8 deletions(-) diff --git a/libs/kimath/src/trigo.cpp b/libs/kimath/src/trigo.cpp index 574c896d80..78ffc1c8d1 100644 --- a/libs/kimath/src/trigo.cpp +++ b/libs/kimath/src/trigo.cpp @@ -39,6 +39,54 @@ // Returns true if the point P is on the segment S. // faster than TestSegmentHit() because P should be exactly on S // therefore works fine only for H, V and 45 deg segm (suitable for wires in eeschema) +/* +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 +*/ + +//#define USE_ALTERNATE_CENTER_ALGO + +#ifdef USE_ALTERNATE_CENTER_ALGO +bool CircleCenterFrom3Points( 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; +} +#endif + bool IsPointOnSegment( const VECTOR2I& aSegStart, const VECTOR2I& aSegEnd, const VECTOR2I& aTestPoint ) { @@ -302,8 +350,8 @@ const VECTOR2D CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aEnd, const EDA_ANGLE& aAngle ) { EDA_ANGLE angle( aAngle ); - VECTOR2I start = aStart; - VECTOR2I end = aEnd; + VECTOR2D start = aStart; + VECTOR2D end = aEnd; if( angle < ANGLE_0 ) { @@ -337,6 +385,7 @@ const VECTOR2D CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aEnd, const VECTOR2D CalcArcCenter( 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; @@ -385,10 +434,19 @@ const VECTOR2D CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, cons bSlope -= std::numeric_limits::epsilon(); } } +#ifdef USE_ALTERNATE_CENTER_ALGO + // We can call ArcCenterFrom3Points from here because special cases are filtered. + CircleCenterFrom3Points( aStart, aMid, aEnd, ¢er ); + return center; +#endif // Prevent divide by zero error + // a small value is used. std::numeric_limits::epsilon() is too small and + // generate false results if( aSlope == 0.0 ) - aSlope = std::numeric_limits::epsilon(); + aSlope = 1e-10; + if( bSlope == 0.0 ) + bSlope = 1e-10; // What follows is the calculation of the center using the slope of the two lines as well as // the propagated error that occurs when rounding to the nearest nanometer. The error can be @@ -423,7 +481,6 @@ const VECTOR2D CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, cons + daSlopeMidEndX * daSlopeMidEndX ); double centerX = ( abSlopeStartEndY + bSlopeStartMidX - aSlopeMidEndX ) / twiceBASlopeDiff; - double dCenterX = centerX * std::sqrt( ( dCenterNumeratorX / centerNumeratorX * dCenterNumeratorX / centerNumeratorX ) + ( dtwiceBASlopeDiff / twiceBASlopeDiff * dtwiceBASlopeDiff / twiceBASlopeDiff ) ); @@ -489,4 +546,3 @@ const VECTOR2I CalcArcCenter( const VECTOR2I& aStart, const VECTOR2I& aMid, cons return iCenter; } - diff --git a/pcbnew/import_gfx/graphics_importer_pcbnew.cpp b/pcbnew/import_gfx/graphics_importer_pcbnew.cpp index 02ee66e6ef..b03d00aa08 100644 --- a/pcbnew/import_gfx/graphics_importer_pcbnew.cpp +++ b/pcbnew/import_gfx/graphics_importer_pcbnew.cpp @@ -120,8 +120,8 @@ void GRAPHICS_IMPORTER_PCBNEW::AddArc( const VECTOR2D& aCenter, const VECTOR2D& // The criteria used here is radius < MAX_INT / 2. // this is not perfect, but we do not know the exact final position of the arc, so // we cannot test the coordinate values, because the arc can be moved before being placed. - VECTOR2D center = CalcArcCenter( arc->GetStart(), arc->GetEnd(), aAngle ); - double radius = ( center - arc->GetStart() ).EuclideanNorm(); + VECTOR2D center = MapCoordinate( aCenter ); + double radius = ( center - MapCoordinate( aStart ) ).EuclideanNorm(); double rd_max_value = std::numeric_limits::max() / 2.0; if( radius >= rd_max_value ) diff --git a/pcbnew/pcb_painter.cpp b/pcbnew/pcb_painter.cpp index 7d3a95c496..f922a287d6 100644 --- a/pcbnew/pcb_painter.cpp +++ b/pcbnew/pcb_painter.cpp @@ -844,7 +844,7 @@ void PCB_PAINTER::draw( const PCB_ARC* aArc, int aLayer ) // Debug only: enable this code only to test the SHAPE_ARC::ConvertToPolyline function // and display the polyline created by it. #if 0 - SHAPE_ARC arc( aArc->GetCenter(), aArc->GetStart(), aArc->GetAngle() / 10.0, aArc->GetWidth() ); + SHAPE_ARC arc( aArc->GetCenter(), aArc->GetStart(), aArc->GetAngle(), aArc->GetWidth() ); SHAPE_LINE_CHAIN arcSpine = arc.ConvertToPolyline( m_maxError ); m_gal->SetLineWidth( m_pcbSettings.m_outlineWidth ); m_gal->SetIsFill( false );