KiCad PCB EDA Suite
bezier_curves.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) 2014 Jean-Pierre Charras, jp.charras at wanadoo.fr
5  * Copyright (C) 2014-2017 KiCad Developers, see CHANGELOG.TXT for contributors.
6  *
7  * This program is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License
9  * as published by the Free Software Foundation; either version 2
10  * of the License, or (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, you may find one here:
19  * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
20  * or you may search the http://www.gnu.org website for the version 2 license,
21  * or you may write to the Free Software Foundation, Inc.,
22  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
23  */
24 
25 /************************************/
26 /* routines to handle bezier curves */
27 /************************************/
28 
29 #include <fctsys.h>
30 #include <bezier_curves.h>
31 
32 
33 static inline double calc_sq_distance( int x1, int y1, int x2, int y2 )
34 {
35  int dx = x2 - x1;
36  int dy = y2 - y1;
37 
38  return (double)dx * dx + (double)dy * dy;
39 }
40 
41 
42 static inline double sqrt_len( int dx, int dy )
43 {
44  return ((double)dx * dx) + ((double)dy * dy);
45 }
46 
47 
48 void BEZIER_POLY::GetPoly( std::vector<wxPoint>& aOutput, int aMinSegLen )
49 {
50  wxCHECK( !m_ctrlPts.empty(), /* void */ );
51  m_minSegLen = std::max( 1, aMinSegLen );
52  m_output = &aOutput;
53  m_output->clear();
54  m_output->push_back( wxPoint( m_ctrlPts.front() ) );
55 
56  // Only quadratic and cubic Bezier curves are handled
57  if( m_ctrlPts.size() == 3 )
59  m_ctrlPts[1].x, m_ctrlPts[1].y,
60  m_ctrlPts[2].x, m_ctrlPts[2].y, 0 );
61 
62  else if( m_ctrlPts.size() == 4 )
64  m_ctrlPts[1].x, m_ctrlPts[1].y,
65  m_ctrlPts[2].x, m_ctrlPts[2].y,
66  m_ctrlPts[3].x, m_ctrlPts[3].y, 0 );
67 
68  m_output->push_back( wxPoint( m_ctrlPts.back() ) );
69 }
70 
71 
72 void BEZIER_POLY::recursiveBezier( int x1, int y1, int x2, int y2,
73  int x3, int y3, unsigned int level )
74 {
75  if( level > recursion_limit )
76  return;
77 
78  // Calculate all the mid-points of the line segments
79  //----------------------
80  int x12 = (x1 + x2) / 2;
81  int y12 = (y1 + y2) / 2;
82  int x23 = (x2 + x3) / 2;
83  int y23 = (y2 + y3) / 2;
84  int x123 = (x12 + x23) / 2;
85  int y123 = (y12 + y23) / 2;
86 
87  int dx = x3 - x1;
88  int dy = y3 - y1;
89  double d = fabs( ((double) (x2 - x3) * dy) - ((double) (y2 - y3) * dx ) );
90  double da;
91 
93  {
94  // Regular case
95  //-----------------
96  if( d * d <= distance_tolerance_square * (dx * dx + dy * dy) )
97  {
98  // If the curvature doesn't exceed the distance_tolerance value
99  // we tend to finish subdivisions.
100  //----------------------
102  {
103  addSegment( wxPoint( x123, y123 ) );
104  return;
105  }
106 
107  // Angle & Cusp Condition
108  //----------------------
109  da = fabs( atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) -
110  atan2( (double) ( y2 - y1 ), (double) ( x2 - x1 ) ) );
111  if( da >=M_PI )
112  da = 2 * M_PI - da;
113 
114  if( da < angle_tolerance )
115  {
116  // Finally we can stop the recursion
117  //----------------------
118  addSegment( wxPoint( x123, y123 ) );
119  return;
120  }
121  }
122  }
123  else
124  {
125  // Collinear case
126  //------------------
127  da = sqrt_len(dx, dy);
128  if( da == 0 )
129  {
130  d = calc_sq_distance( x1, y1, x2, y2 );
131  }
132  else
133  {
134  d = ( (double)(x2 - x1) * dx + (double)(y2 - y1) * dy ) / da;
135  if( d > 0 && d < 1 )
136  {
137  // Simple collinear case, 1---2---3
138  // We can leave just two endpoints
139  return;
140  }
141  if( d <= 0 )
142  d = calc_sq_distance( x2, y2, x1, y1 );
143  else if( d >= 1 )
144  d = calc_sq_distance( x2, y2, x3, y3 );
145  else
146  d = calc_sq_distance( x2, y2, x1 + (int) d * dx,
147  y1 + (int) d * dy );
148  }
149  if( d < distance_tolerance_square )
150  {
151  addSegment( wxPoint( x2, y2 ) );
152  return;
153  }
154  }
155 
156  // Continue subdivision
157  //----------------------
158  recursiveBezier( x1, y1, x12, y12, x123, y123, level + 1 );
159  recursiveBezier( x123, y123, x23, y23, x3, y3, level + 1 );
160 }
161 
162 
163 void BEZIER_POLY::recursiveBezier( int x1, int y1, int x2, int y2,
164  int x3, int y3, int x4, int y4, unsigned int level )
165 {
166  if( level > recursion_limit )
167  return;
168 
169  // Calculate all the mid-points of the line segments
170  //----------------------
171  int x12 = (x1 + x2) / 2;
172  int y12 = (y1 + y2) / 2;
173  int x23 = (x2 + x3) / 2;
174  int y23 = (y2 + y3) / 2;
175  int x34 = (x3 + x4) / 2;
176  int y34 = (y3 + y4) / 2;
177  int x123 = (x12 + x23) / 2;
178  int y123 = (y12 + y23) / 2;
179  int x234 = (x23 + x34) / 2;
180  int y234 = (y23 + y34) / 2;
181  int x1234 = (x123 + x234) / 2;
182  int y1234 = (y123 + y234) / 2;
183 
184 
185  // Try to approximate the full cubic curve by a single straight line
186  //------------------
187  int dx = x4 - x1;
188  int dy = y4 - y1;
189 
190  double d2 = fabs( (double) ( (x2 - x4) * dy - (y2 - y4) * dx ) );
191  double d3 = fabs( (double) ( (x3 - x4) * dy - (y3 - y4) * dx ) );
192  double da1, da2, k;
193 
194  switch( (int(d2 > curve_collinearity_epsilon) << 1) +
195  int(d3 > curve_collinearity_epsilon) )
196  {
197  case 0:
198 
199  // All collinear OR p1==p4
200  //----------------------
201  k = dx * dx + dy * dy;
202  if( k == 0 )
203  {
204  d2 = calc_sq_distance( x1, y1, x2, y2 );
205  d3 = calc_sq_distance( x4, y4, x3, y3 );
206  }
207  else
208  {
209  k = 1 / k;
210  da1 = x2 - x1;
211  da2 = y2 - y1;
212  d2 = k * (da1 * dx + da2 * dy);
213  da1 = x3 - x1;
214  da2 = y3 - y1;
215  d3 = k * (da1 * dx + da2 * dy);
216  if( d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1 )
217  {
218  // Simple collinear case, 1---2---3---4
219  // We can leave just two endpoints
220  return;
221  }
222  if( d2 <= 0 )
223  d2 = calc_sq_distance( x2, y2, x1, y1 );
224  else if( d2 >= 1 )
225  d2 = calc_sq_distance( x2, y2, x4, y4 );
226  else
227  d2 = calc_sq_distance( x2, y2, x1 + (int) d2 * dx,
228  y1 + (int) d2 * dy );
229 
230  if( d3 <= 0 )
231  d3 = calc_sq_distance( x3, y3, x1, y1 );
232  else if( d3 >= 1 )
233  d3 = calc_sq_distance( x3, y3, x4, y4 );
234  else
235  d3 = calc_sq_distance( x3, y3, x1 + (int) d3 * dx,
236  y1 + (int) d3 * dy );
237  }
238  if( d2 > d3 )
239  {
240  if( d2 < distance_tolerance_square )
241  {
242  addSegment( wxPoint( x2, y2 ) );
243  return;
244  }
245  }
246  else
247  {
248  if( d3 < distance_tolerance_square )
249  {
250  addSegment( wxPoint( x3, y3 ) );
251  return;
252  }
253  }
254  break;
255 
256  case 1:
257 
258  // p1,p2,p4 are collinear, p3 is significant
259  //----------------------
260  if( d3 * d3 <= distance_tolerance_square * sqrt_len(dx, dy) )
261  {
263  {
264  addSegment( wxPoint( x23, y23 ) );
265  return;
266  }
267 
268  // Angle Condition
269  //----------------------
270  da1 = fabs( atan2( (double) ( y4 - y3 ), (double) ( x4 - x3 ) ) -
271  atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) );
272  if( da1 >= M_PI )
273  da1 = 2 * M_PI - da1;
274 
275  if( da1 < angle_tolerance )
276  {
277  addSegment( wxPoint( x2, y2 ) );
278  addSegment( wxPoint( x3, y3 ) );
279  return;
280  }
281 
282  if( cusp_limit != 0.0 )
283  {
284  if( da1 > cusp_limit )
285  {
286  addSegment( wxPoint( x3, y3 ) );
287  return;
288  }
289  }
290  }
291  break;
292 
293  case 2:
294 
295  // p1,p3,p4 are collinear, p2 is significant
296  //----------------------
297  if( d2 * d2 <= distance_tolerance_square * sqrt_len(dx, dy) )
298  {
300  {
301  addSegment( wxPoint( x23, y23 ) );
302  return;
303  }
304 
305  // Angle Condition
306  //----------------------
307  da1 = fabs( atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) -
308  atan2( (double) ( y2 - y1 ), (double) ( x2 - x1 ) ) );
309  if( da1 >= M_PI )
310  da1 = 2 * M_PI - da1;
311 
312  if( da1 < angle_tolerance )
313  {
314  addSegment( wxPoint( x2, y2 ) );
315  addSegment( wxPoint( x3, y3 ) );
316  return;
317  }
318 
319  if( cusp_limit != 0.0 )
320  {
321  if( da1 > cusp_limit )
322  {
323  addSegment( wxPoint( x2, y2 ) );
324  return;
325  }
326  }
327  }
328  break;
329 
330  case 3:
331 
332  // Regular case
333  //-----------------
334  if( (d2 + d3) * (d2 + d3) <= distance_tolerance_square * sqrt_len(dx, dy) )
335  {
336  // If the curvature doesn't exceed the distance_tolerance value
337  // we tend to finish subdivisions.
338  //----------------------
340  {
341  addSegment( wxPoint( x23, y23 ) );
342  return;
343  }
344 
345  // Angle & Cusp Condition
346  //----------------------
347  k = atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) );
348  da1 = fabs( k - atan2( (double) ( y2 - y1 ),
349  (double) ( x2 - x1 ) ) );
350  da2 = fabs( atan2( (double) ( y4 - y3 ),
351  (double) ( x4 - x3 ) ) - k );
352  if( da1 >= M_PI )
353  da1 = 2 * M_PI - da1;
354  if( da2 >= M_PI )
355  da2 = 2 * M_PI - da2;
356 
357  if( da1 + da2 < angle_tolerance )
358  {
359  // Finally we can stop the recursion
360  //----------------------
361  addSegment( wxPoint( x23, y23 ) );
362  return;
363  }
364 
365  if( cusp_limit != 0.0 )
366  {
367  if( da1 > cusp_limit )
368  {
369  addSegment( wxPoint( x2, y2 ) );
370  return;
371  }
372 
373  if( da2 > cusp_limit )
374  {
375  addSegment( wxPoint( x3, y3 ) );
376  return;
377  }
378  }
379  }
380  break;
381  }
382 
383  // Continue subdivision
384  //----------------------
385  recursiveBezier( x1, y1, x12, y12, x123, y123, x1234, y1234, level + 1 );
386  recursiveBezier( x1234, y1234, x234, y234, x34, y34, x4, y4, level + 1 );
387 }
static constexpr double distance_tolerance_square
static constexpr double angle_tolerance
void GetPoly(std::vector< wxPoint > &aOutput, int aMinSegLen=0)
Converts Bezier curve to a polygon.
void recursiveBezier(int x1, int y1, int x2, int y2, int x3, int y3, unsigned int level)
std::vector< wxPoint > m_ctrlPts
Control points
Definition: bezier_curves.h:79
void addSegment(const wxPoint &aSegment)
Definition: bezier_curves.h:84
static double sqrt_len(int dx, int dy)
static constexpr double curve_collinearity_epsilon
static constexpr double cusp_limit
std::vector< wxPoint > * m_output
Pointer to the output vector
Definition: bezier_curves.h:82
static constexpr double curve_angle_tolerance_epsilon
#define max(a, b)
Definition: auxiliary.h:86
static constexpr unsigned int recursion_limit
static double calc_sq_distance(int x1, int y1, int x2, int y2)