KiCad PCB EDA Suite
Loading...
Searching...
No Matches
circle.cpp
Go to the documentation of this file.
1/*
2 * This program source code file is part of KiCad, a free EDA CAD application.
3 *
4 * Copyright (C) 2021 Roberto Fernandez Bautista <[email protected]>
5 * Copyright (C) 2021-2022 KiCad Developers, see AUTHORS.txt for contributors.
6 *
7 * This program is free software: you can redistribute it and/or modify it
8 * under the terms of the GNU General Public License as published by the
9 * Free Software Foundation, either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License along
18 * with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21#include <geometry/circle.h>
22#include <geometry/seg.h>
23#include <geometry/shape.h> // for MIN_PRECISION_IU
24#include <math/util.h> // for KiROUND
25#include <math/vector2d.h> // for VECTOR2I
26#include <math.h> // for sqrt
27#include <trigo.h> // for CalcArcMid
28
29
31{
32 Center = { 0, 0 };
33 Radius = 0;
34}
35
36
37CIRCLE::CIRCLE( const VECTOR2I& aCenter, int aRadius )
38{
39 Center = aCenter;
40 Radius = aRadius;
41}
42
43
44CIRCLE::CIRCLE( const CIRCLE& aOther )
45{
46 Center = aOther.Center;
47 Radius = aOther.Radius;
48}
49
50
51CIRCLE& CIRCLE::ConstructFromTanTanPt( const SEG& aLineA, const SEG& aLineB, const VECTOR2I& aP )
52{
53 //fixme: There might be more efficient / accurate solution than using geometrical constructs
54
55 SEG anglebisector;
56 VECTOR2I intersectPoint;
57
58 auto furthestFromIntersect =
59 [&]( VECTOR2I aPt1, VECTOR2I aPt2 ) -> VECTOR2I
60 {
61 if( ( aPt1 - intersectPoint ).EuclideanNorm()
62 > ( aPt2 - intersectPoint ).EuclideanNorm() )
63 {
64 return aPt1;
65 }
66 else
67 {
68 return aPt2;
69 }
70 };
71
72 auto closestToIntersect =
73 [&]( VECTOR2I aPt1, VECTOR2I aPt2 ) -> VECTOR2I
74 {
75 if( ( aPt1 - intersectPoint ).EuclideanNorm()
76 <= ( aPt2 - intersectPoint ).EuclideanNorm() )
77 {
78 return aPt1;
79 }
80 else
81 {
82 return aPt2;
83 }
84 };
85
86 if( aLineA.ApproxParallel( aLineB ) )
87 {
88 // Special case, no intersection point between the two lines
89 // The center will be in the line equidistant between the two given lines
90 // The radius will be half the distance between the two lines
91 // The possible centers can be found by intersection
92
93 SEG perpendicularAtoB( aLineA.A, aLineB.LineProject( aLineA.A ) );
94 VECTOR2I midPt = perpendicularAtoB.Center();
95 Radius = ( midPt - aLineA.A ).EuclideanNorm();
96
97 anglebisector = aLineA.ParallelSeg( midPt );
98
99 Center = aP; // use this circle as a construction to find the actual centers
100 std::vector<VECTOR2I> possibleCenters = IntersectLine( anglebisector );
101
102 wxCHECK_MSG( possibleCenters.size() > 0, *this, wxT( "No solutions exist!" ) );
103 intersectPoint = aLineA.A; // just for the purpose of deciding which solution to return
104
105 // For the special case of the two segments being parallel, we will return the solution
106 // whose center is closest to aLineA.A
107 Center = closestToIntersect( possibleCenters.front(), possibleCenters.back() );
108 }
109 else
110 {
111 // General case, using homothety.
112 // All circles inscribed in the same angle are homothetic with center at the intersection
113 // In this code, the prefix "h" denotes "the homothetic image"
114 OPT_VECTOR2I intersectCalc = aLineA.IntersectLines( aLineB );
115 wxCHECK_MSG( intersectCalc, *this, wxT( "Lines do not intersect but are not parallel?" ) );
116 intersectPoint = *intersectCalc;
117
118 if( aP == intersectPoint )
119 {
120 //Special case: The point is at the intersection of the two lines
121 Center = aP;
122 Radius = 0;
123 return *this;
124 }
125
126 // Calculate bisector
127 VECTOR2I lineApt = furthestFromIntersect( aLineA.A, aLineA.B );
128 VECTOR2I lineBpt = furthestFromIntersect( aLineB.A, aLineB.B );
129 VECTOR2I bisectorPt = CalcArcMid( lineApt, lineBpt, intersectPoint, true );
130
131 anglebisector.A = intersectPoint;
132 anglebisector.B = bisectorPt;
133
134 // Create an arbitrary circle that is tangent to both lines
135 CIRCLE hSolution;
136 hSolution.Center = anglebisector.LineProject( aP );
137 hSolution.Radius = aLineA.LineDistance( hSolution.Center );
138
139 // Find the homothetic image of aP in the construction circle (hSolution)
140 SEG throughaP( intersectPoint, aP );
141 std::vector<VECTOR2I> hProjections = hSolution.IntersectLine( throughaP );
142 wxCHECK_MSG( hProjections.size() > 0, *this, wxT( "No solutions exist!" ) );
143
144 // We want to create a fillet, so the projection of homothetic projection of aP
145 // should be the one closest to the intersection
146 VECTOR2I hSelected = closestToIntersect( hProjections.front(), hProjections.back() );
147
148 VECTOR2I hTanLineA = aLineA.LineProject( hSolution.Center );
149 VECTOR2I hTanLineB = aLineB.LineProject( hSolution.Center );
150
151 // To minimise errors, use the furthest away tangent point from aP
152 if( ( hTanLineA - aP ).EuclideanNorm() > ( hTanLineB - aP ).EuclideanNorm() )
153 {
154 // Find the tangent at line A by homothetic inversion
155 SEG hT( hTanLineA, hSelected );
156 OPT_VECTOR2I actTanA = hT.ParallelSeg( aP ).IntersectLines( aLineA );
157 wxCHECK_MSG( actTanA, *this, wxT( "No solutions exist!" ) );
158
159 // Find circle center by perpendicular intersection with the angle bisector
160 SEG perpendicularToTanA = aLineA.PerpendicularSeg( *actTanA );
161 OPT_VECTOR2I actCenter = perpendicularToTanA.IntersectLines( anglebisector );
162 wxCHECK_MSG( actCenter, *this, wxT( "No solutions exist!" ) );
163
164 Center = *actCenter;
165 Radius = aLineA.LineDistance( Center );
166 }
167 else
168 {
169 // Find the tangent at line B by inversion
170 SEG hT( hTanLineB, hSelected );
171 OPT_VECTOR2I actTanB = hT.ParallelSeg( aP ).IntersectLines( aLineB );
172 wxCHECK_MSG( actTanB, *this, wxT( "No solutions exist!" ) );
173
174 // Find circle center by perpendicular intersection with the angle bisector
175 SEG perpendicularToTanB = aLineB.PerpendicularSeg( *actTanB );
176 OPT_VECTOR2I actCenter = perpendicularToTanB.IntersectLines( anglebisector );
177 wxCHECK_MSG( actCenter, *this, wxT( "No solutions exist!" ) );
178
179 Center = *actCenter;
180 Radius = aLineB.LineDistance( Center );
181 }
182 }
183
184 return *this;
185}
186
187
188bool CIRCLE::Contains( const VECTOR2I& aP ) const
189{
190 int distance = ( aP - Center ).EuclideanNorm();
191
192 return distance <= ( (int64_t) Radius + SHAPE::MIN_PRECISION_IU )
193 && distance >= ( (int64_t) Radius - SHAPE::MIN_PRECISION_IU );
194}
195
196
198{
199 VECTOR2I vec = aP - Center;
200
201 // Handle special case where aP is equal to this circle's center
202 if( vec.x == 0 && vec.y == 0 )
203 vec.x = 1; // Arbitrary, to ensure the return value is always on the circumference
204
205 return vec.Resize( Radius ) + Center;
206}
207
208
209std::vector<VECTOR2I> CIRCLE::Intersect( const CIRCLE& aCircle ) const
210{
211 // From https://mathworld.wolfram.com/Circle-CircleIntersection.html
212 //
213 // Simplify the problem:
214 // Let this circle be centered at (0,0), with radius r1
215 // Let aCircle be centered at (d, 0), with radius r2
216 // (i.e. d is the distance between the two circle centers)
217 //
218 // The equations of the two circles are
219 // (1) x^2 + y^2 = r1^2
220 // (2) (x - d)^2 + y^2 = r2^2
221 //
222 // Combining (1) into (2):
223 // (x - d)^2 + r1^2 - x^2 = r2^2
224 // Expanding:
225 // x^2 - 2*d*x + d^2 + r1^2 - x^2 = r2^2
226 // Rearranging for x:
227 // (3) x = (d^2 + r1^2 - r2^2) / (2 * d)
228 //
229 // Rearranging (1) gives:
230 // (4) y = sqrt(r1^2 - x^2)
231
232 std::vector<VECTOR2I> retval;
233
234 VECTOR2I vecCtoC = aCircle.Center - Center;
235 int64_t d = vecCtoC.EuclideanNorm();
236 int64_t r1 = Radius;
237 int64_t r2 = aCircle.Radius;
238
239 if( d > ( r1 + r2 ) || ( d < ( std::abs( r1 - r2 ) ) ) )
240 return retval; //circles do not intersect
241
242 if( d == 0 )
243 return retval; // circles are co-centered. Don't return intersection points
244
245 // Equation (3)
246 int64_t x = ( ( d * d ) + ( r1 * r1 ) - ( r2 * r2 ) ) / ( int64_t( 2 ) * d );
247 int64_t r1sqMinusXsq = ( r1 * r1 ) - ( x * x );
248
249 if( r1sqMinusXsq < 0 )
250 return retval; //circles do not intersect
251
252 // Equation (4)
253 int64_t y = KiROUND( sqrt( r1sqMinusXsq ) );
254
255 // Now correct back to original coordinates
256 EDA_ANGLE rotAngle( vecCtoC );
257 VECTOR2I solution1( x, y );
258 RotatePoint( solution1, -rotAngle );
259 solution1 += Center;
260 retval.push_back( solution1 );
261
262 if( y != 0 )
263 {
264 VECTOR2I solution2( x, -y );
265 RotatePoint( solution2, -rotAngle );
266 solution2 += Center;
267 retval.push_back( solution2 );
268 }
269
270 return retval;
271}
272
273
274std::vector<VECTOR2I> CIRCLE::Intersect( const SEG& aSeg ) const
275{
276 std::vector<VECTOR2I> retval;
277
278 for( VECTOR2I& intersection : IntersectLine( aSeg ) )
279 {
280 if( aSeg.Contains( intersection ) )
281 retval.push_back( intersection );
282 }
283
284 return retval;
285}
286
287
288std::vector<VECTOR2I> CIRCLE::IntersectLine( const SEG& aLine ) const
289{
290 std::vector<VECTOR2I> retval;
291
292 //
293 // . * .
294 // * *
295 // -----1-------m-------2----
296 // * *
297 // * O *
298 // * *
299 // * *
300 // * *
301 // * *
302 // ' * '
303 // Let O be the center of this circle, 1 and 2 the intersection points of the line
304 // and M be the center of the chord connecting points 1 and 2
305 //
306 // M will be O projected perpendicularly to the line since a chord is always perpendicular
307 // to the radius.
308 //
309 // The distance M1 = M2 can be computed by pythagoras since O1 = O2 = Radius
310 //
311 // M1= M2 = sqrt( Radius^2 - OM^2)
312 //
313
314 VECTOR2I m = aLine.LineProject( Center ); // O projected perpendicularly to the line
315 int64_t omDist = ( m - Center ).EuclideanNorm();
316
317 if( omDist > ( (int64_t) Radius + SHAPE::MIN_PRECISION_IU ) )
318 {
319 return retval; // does not intersect
320 }
321 else if( omDist <= ( (int64_t) Radius + SHAPE::MIN_PRECISION_IU )
322 && omDist >= ( (int64_t) Radius - SHAPE::MIN_PRECISION_IU ) )
323 {
324 retval.push_back( m );
325 return retval; //tangent
326 }
327
328 int64_t radiusSquared = (int64_t) Radius * (int64_t) Radius;
329 int64_t omDistSquared = omDist * omDist;
330
331 int mTo1dist = sqrt( radiusSquared - omDistSquared );
332
333 VECTOR2I mTo1vec = ( aLine.B - aLine.A ).Resize( mTo1dist );
334 VECTOR2I mTo2vec = -mTo1vec;
335
336 retval.push_back( mTo1vec + m );
337 retval.push_back( mTo2vec + m );
338
339 return retval;
340}
341
342
343bool CIRCLE::Contains( const VECTOR2I& aP )
344{
345 return ( aP - Center ).EuclideanNorm() < Radius;
346}
Represent basic circle geometry with utility geometry functions.
Definition: circle.h:33
VECTOR2I Center
Public to make access simpler.
Definition: circle.h:116
int Radius
Public to make access simpler.
Definition: circle.h:115
std::vector< VECTOR2I > Intersect(const CIRCLE &aCircle) const
Compute the intersection points between this circle and aCircle.
Definition: circle.cpp:209
CIRCLE()
Definition: circle.cpp:30
std::vector< VECTOR2I > IntersectLine(const SEG &aLine) const
Compute the intersection points between this circle and aLine.
Definition: circle.cpp:288
CIRCLE & ConstructFromTanTanPt(const SEG &aLineA, const SEG &aLineB, const VECTOR2I &aP)
Construct this circle such that it is tangent to the given segments and passes through the given poin...
Definition: circle.cpp:51
VECTOR2I NearestPoint(const VECTOR2I &aP) const
Compute the point on the circumference of the circle that is the closest to aP.
Definition: circle.cpp:197
bool Contains(const VECTOR2I &aP) const
Return true if aP is on the circumference of this circle.
Definition: circle.cpp:188
Definition: seg.h:42
VECTOR2I A
Definition: seg.h:49
int LineDistance(const VECTOR2I &aP, bool aDetermineSide=false) const
Return the closest Euclidean distance between point aP and the line defined by the ends of segment (t...
Definition: seg.cpp:341
VECTOR2I B
Definition: seg.h:50
VECTOR2I Center() const
Definition: seg.h:362
bool ApproxParallel(const SEG &aSeg, int aDistanceThreshold=1) const
Definition: seg.cpp:402
SEG ParallelSeg(const VECTOR2I &aP) const
Compute a segment parallel to this one, passing through point aP.
Definition: seg.cpp:216
OPT_VECTOR2I IntersectLines(const SEG &aSeg) const
Compute the intersection point of lines passing through ends of (this) and aSeg.
Definition: seg.h:210
bool Contains(const SEG &aSeg) const
Definition: seg.h:307
VECTOR2I LineProject(const VECTOR2I &aP) const
Compute the perpendicular projection point of aP on a line passing through ends of the segment.
Definition: seg.cpp:312
SEG PerpendicularSeg(const VECTOR2I &aP) const
Compute a segment perpendicular to this one, passing through point aP.
Definition: seg.cpp:207
static const int MIN_PRECISION_IU
This is the minimum precision for all the points in a shape.
Definition: shape.h:131
T EuclideanNorm() const
Compute the Euclidean norm of the vector, which is defined as sqrt(x ** 2 + y ** 2).
Definition: vector2d.h:265
VECTOR2< T > Resize(T aNewLength) const
Return a vector of the same direction, but length specified in aNewLength.
Definition: vector2d.h:350
EDA_ANGLE abs(const EDA_ANGLE &aAngle)
Definition: eda_angle.h:424
static float distance(const SFVEC2UI &a, const SFVEC2UI &b)
std::optional< VECTOR2I > OPT_VECTOR2I
Definition: seg.h:39
const VECTOR2I CalcArcMid(const VECTOR2I &aStart, const VECTOR2I &aEnd, const VECTOR2I &aCenter, bool aMinArcAngle=true)
Return the middle point of an arc, half-way between aStart and aEnd.
Definition: trigo.cpp:208
void RotatePoint(int *pX, int *pY, const EDA_ANGLE &aAngle)
Calculate the new point of coord coord pX, pY, for a rotation center 0, 0.
Definition: trigo.cpp:228
double EuclideanNorm(const VECTOR2I &vector)
Definition: trigo.h:128
constexpr ret_type KiROUND(fp_type v)
Round a floating point number to an integer using "round halfway cases away from zero".
Definition: util.h:85