KiCad PCB EDA Suite
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 <roberto.fer.bau@gmail.com>
5  * Copyright (C) 2021 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 GetArcMid
28 
29 
31 {
32  Center = { 0, 0 };
33  Radius = 0;
34 }
35 
36 
37 CIRCLE::CIRCLE( const VECTOR2I& aCenter, int aRadius )
38 {
39  Center = aCenter;
40  Radius = aRadius;
41 }
42 
43 
44 CIRCLE::CIRCLE( const CIRCLE& aOther )
45 {
46  Center = aOther.Center;
47  Radius = aOther.Radius;
48 }
49 
50 
51 CIRCLE& 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, "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, "Lines do not intersect but are not parallel?" );
116  intersectPoint = intersectCalc.get();
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 = GetArcMid( 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, "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, "No solutions exist!" );
158 
159  // Find circle center by perpendicular intersection with the angle bisector
160  SEG perpendicularToTanA = aLineA.PerpendicularSeg( actTanA.get() );
161  OPT_VECTOR2I actCenter = perpendicularToTanA.IntersectLines( anglebisector );
162  wxCHECK_MSG( actCenter, *this, "No solutions exist!" );
163 
164  Center = actCenter.get();
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, "No solutions exist!" );
173 
174  // Find circle center by perpendicular intersection with the angle bisector
175  SEG perpendicularToTanB = aLineB.PerpendicularSeg( actTanB.get() );
176  OPT_VECTOR2I actCenter = perpendicularToTanB.IntersectLines( anglebisector );
177  wxCHECK_MSG( actCenter, *this, "No solutions exist!" );
178 
179  Center = actCenter.get();
180  Radius = aLineB.LineDistance( Center );
181  }
182  }
183 
184  return *this;
185 }
186 
187 
188 bool 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 
209 std::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  double rotAngle = vecCtoC.Angle();
257  VECTOR2I solution1( x, y );
258  solution1 = solution1.Rotate( rotAngle );
259  solution1 += Center;
260  retval.push_back( solution1 );
261 
262  if( y != 0 )
263  {
264  VECTOR2I solution2( x, -y );
265  solution2 = solution2.Rotate( rotAngle );
266  solution2 += Center;
267  retval.push_back( solution2 );
268  }
269 
270  return retval;
271 }
272 
273 
274 std::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 
288 std::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  // O1= O2 = sqrt( Radius^2 - OM^2)
312  //
313 
314  VECTOR2I m = aLine.LineProject( Center );
315  int64_t om = ( m - Center ).EuclideanNorm();
316 
317  if( om > ( (int64_t) Radius + SHAPE::MIN_PRECISION_IU ) )
318  {
319  return retval; // does not intersect
320  }
321  else if( om <= ( (int64_t) Radius + SHAPE::MIN_PRECISION_IU )
322  && om >= ( (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 omSquared = om * om;
330 
331  int mTo1 = sqrt( radiusSquared - omSquared );
332 
333  VECTOR2I mTo1vec = ( aLine.B - aLine.A ).Resize( mTo1 );
334  VECTOR2I mTo2vec = -mTo1vec;
335 
336  retval.push_back( mTo1vec + m );
337  retval.push_back( mTo2vec + m );
338 
339  return retval;
340 }
341 
342 
343 bool CIRCLE::Contains( const VECTOR2I& aP )
344 {
345  return ( aP - Center ).EuclideanNorm() < Radius;
346 }
double EuclideanNorm(const wxPoint &vector)
Euclidean norm of a 2D vector.
Definition: trigo.h:146
SEG PerpendicularSeg(const VECTOR2I &aP) const
Compute a segment perpendicular to this one, passing through point aP.
Definition: seg.cpp:165
OPT_VECTOR2I IntersectLines(const SEG &aSeg) const
Compute the intersection point of lines passing through ends of (this) and aSeg.
Definition: seg.h:209
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:297
SEG ParallelSeg(const VECTOR2I &aP) const
Compute a segment parallel to this one, passing through point aP.
Definition: seg.cpp:174
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 LineProject(const VECTOR2I &aP) const
Compute the perpendicular projection point of aP on a line passing through ends of the segment.
Definition: seg.cpp:268
OPT< VECTOR2I > OPT_VECTOR2I
Definition: seg.h:38
const VECTOR2I GetArcMid(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:163
bool ApproxParallel(const SEG &aSeg) const
Definition: seg.h:290
static const int MIN_PRECISION_IU
This is the minimum precision for all the points in a shape.
Definition: shape.h:122
std::vector< VECTOR2I > IntersectLine(const SEG &aLine) const
Compute the intersection points between this circle and aLine.
Definition: circle.cpp:288
int Radius
Public to make access simpler.
Definition: circle.h:115
Represent basic circle geometry with utility geometry functions.
Definition: circle.h:32
static float distance(const SFVEC2UI &a, const SFVEC2UI &b)
double Angle() const
Compute the angle of the vector.
Definition: vector2d.h:307
Definition: seg.h:40
VECTOR2< T > Resize(T aNewLength) const
Return a vector of the same direction, but length specified in aNewLength.
Definition: vector2d.h:404
VECTOR2< T > Rotate(double aAngle) const
Rotate the vector by a given angle.
Definition: vector2d.h:371
std::vector< VECTOR2I > Intersect(const CIRCLE &aCircle) const
Compute the intersection points between this circle and aCircle.
Definition: circle.cpp:209
bool Contains(const VECTOR2I &aP) const
Return true if aP is on the circumference of this circle.
Definition: circle.cpp:188
VECTOR2I A
Definition: seg.h:48
VECTOR2I Center
Public to make access simpler.
Definition: circle.h:116
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:73
T EuclideanNorm() const
Compute the Euclidean norm of the vector, which is defined as sqrt(x ** 2 + y ** 2).
Definition: vector2d.h:293
CIRCLE()
Definition: circle.cpp:30
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 SEG &aSeg) const
Definition: seg.h:331
VECTOR2I B
Definition: seg.h:49