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 )
49 {
50  wxCHECK( !m_ctrlPts.empty(), /* void */ );
51  m_output = &aOutput;
52  m_output->clear();
53  m_output->push_back( wxPoint( m_ctrlPts.front() ) );
54 
55  // Only quadratic and cubic Bezier curves are handled
56  if( m_ctrlPts.size() == 3 )
58  m_ctrlPts[1].x, m_ctrlPts[1].y,
59  m_ctrlPts[2].x, m_ctrlPts[2].y, 0 );
60 
61  else if( m_ctrlPts.size() == 4 )
63  m_ctrlPts[1].x, m_ctrlPts[1].y,
64  m_ctrlPts[2].x, m_ctrlPts[2].y,
65  m_ctrlPts[3].x, m_ctrlPts[3].y, 0 );
66 
67  m_output->push_back( wxPoint( m_ctrlPts.back() ) );
68 }
69 
70 
71 void BEZIER_POLY::recursiveBezier( int x1, int y1, int x2, int y2, int x3, int y3, unsigned int level )
72 {
73  if( level > recursion_limit )
74  return;
75 
76  // Calculate all the mid-points of the line segments
77  //----------------------
78  int x12 = (x1 + x2) / 2;
79  int y12 = (y1 + y2) / 2;
80  int x23 = (x2 + x3) / 2;
81  int y23 = (y2 + y3) / 2;
82  int x123 = (x12 + x23) / 2;
83  int y123 = (y12 + y23) / 2;
84 
85  int dx = x3 - x1;
86  int dy = y3 - y1;
87  double d = fabs( ((double) (x2 - x3) * dy) - ((double) (y2 - y3) * dx ) );
88  double da;
89 
91  {
92  // Regular case
93  //-----------------
94  if( d * d <= distance_tolerance_square * (dx * dx + dy * dy) )
95  {
96  // If the curvature doesn't exceed the distance_tolerance value
97  // we tend to finish subdivisions.
98  //----------------------
100  {
101  addSegment( wxPoint( x123, y123 ) );
102  return;
103  }
104 
105  // Angle & Cusp Condition
106  //----------------------
107  da = fabs( atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) -
108  atan2( (double) ( y2 - y1 ), (double) ( x2 - x1 ) ) );
109  if( da >=M_PI )
110  da = 2 * M_PI - da;
111 
112  if( da < angle_tolerance )
113  {
114  // Finally we can stop the recursion
115  //----------------------
116  addSegment( wxPoint( x123, y123 ) );
117  return;
118  }
119  }
120  }
121  else
122  {
123  // Collinear case
124  //------------------
125  da = sqrt_len(dx, dy);
126  if( da == 0 )
127  {
128  d = calc_sq_distance( x1, y1, x2, y2 );
129  }
130  else
131  {
132  d = ( (double)(x2 - x1) * dx + (double)(y2 - y1) * dy ) / da;
133  if( d > 0 && d < 1 )
134  {
135  // Simple collinear case, 1---2---3
136  // We can leave just two endpoints
137  return;
138  }
139  if( d <= 0 )
140  d = calc_sq_distance( x2, y2, x1, y1 );
141  else if( d >= 1 )
142  d = calc_sq_distance( x2, y2, x3, y3 );
143  else
144  d = calc_sq_distance( x2, y2, x1 + (int) d * dx,
145  y1 + (int) d * dy );
146  }
147  if( d < distance_tolerance_square )
148  {
149  addSegment( wxPoint( x2, y2 ) );
150  return;
151  }
152  }
153 
154  // Continue subdivision
155  //----------------------
156  recursiveBezier( x1, y1, x12, y12, x123, y123, level + 1 );
157  recursiveBezier( x123, y123, x23, y23, x3, y3, -(level + 1) );
158 }
159 
160 
161 void BEZIER_POLY::recursiveBezier( int x1, int y1, int x2, int y2, int x3, int y3, int x4, int y4, unsigned int level )
162 {
163  if( level > recursion_limit )
164  return;
165 
166  // Calculate all the mid-points of the line segments
167  //----------------------
168  int x12 = (x1 + x2) / 2;
169  int y12 = (y1 + y2) / 2;
170  int x23 = (x2 + x3) / 2;
171  int y23 = (y2 + y3) / 2;
172  int x34 = (x3 + x4) / 2;
173  int y34 = (y3 + y4) / 2;
174  int x123 = (x12 + x23) / 2;
175  int y123 = (y12 + y23) / 2;
176  int x234 = (x23 + x34) / 2;
177  int y234 = (y23 + y34) / 2;
178  int x1234 = (x123 + x234) / 2;
179  int y1234 = (y123 + y234) / 2;
180 
181 
182  // Try to approximate the full cubic curve by a single straight line
183  //------------------
184  int dx = x4 - x1;
185  int dy = y4 - y1;
186 
187  double d2 = fabs( (double) ( (x2 - x4) * dy - (y2 - y4) * dx ) );
188  double d3 = fabs( (double) ( (x3 - x4) * dy - (y3 - y4) * dx ) );
189  double da1, da2, k;
190 
191  switch( (int(d2 > curve_collinearity_epsilon) << 1) +
192  int(d3 > curve_collinearity_epsilon) )
193  {
194  case 0:
195 
196  // All collinear OR p1==p4
197  //----------------------
198  k = dx * dx + dy * dy;
199  if( k == 0 )
200  {
201  d2 = calc_sq_distance( x1, y1, x2, y2 );
202  d3 = calc_sq_distance( x4, y4, x3, y3 );
203  }
204  else
205  {
206  k = 1 / k;
207  da1 = x2 - x1;
208  da2 = y2 - y1;
209  d2 = k * (da1 * dx + da2 * dy);
210  da1 = x3 - x1;
211  da2 = y3 - y1;
212  d3 = k * (da1 * dx + da2 * dy);
213  if( d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1 )
214  {
215  // Simple collinear case, 1---2---3---4
216  // We can leave just two endpoints
217  return;
218  }
219  if( d2 <= 0 )
220  d2 = calc_sq_distance( x2, y2, x1, y1 );
221  else if( d2 >= 1 )
222  d2 = calc_sq_distance( x2, y2, x4, y4 );
223  else
224  d2 = calc_sq_distance( x2, y2, x1 + (int) d2 * dx,
225  y1 + (int) d2 * dy );
226 
227  if( d3 <= 0 )
228  d3 = calc_sq_distance( x3, y3, x1, y1 );
229  else if( d3 >= 1 )
230  d3 = calc_sq_distance( x3, y3, x4, y4 );
231  else
232  d3 = calc_sq_distance( x3, y3, x1 + (int) d3 * dx,
233  y1 + (int) d3 * dy );
234  }
235  if( d2 > d3 )
236  {
237  if( d2 < distance_tolerance_square )
238  {
239  addSegment( wxPoint( x2, y2 ) );
240  return;
241  }
242  }
243  else
244  {
245  if( d3 < distance_tolerance_square )
246  {
247  addSegment( wxPoint( x3, y3 ) );
248  return;
249  }
250  }
251  break;
252 
253  case 1:
254 
255  // p1,p2,p4 are collinear, p3 is significant
256  //----------------------
257  if( d3 * d3 <= distance_tolerance_square * sqrt_len(dx, dy) )
258  {
260  {
261  addSegment( wxPoint( x23, y23 ) );
262  return;
263  }
264 
265  // Angle Condition
266  //----------------------
267  da1 = fabs( atan2( (double) ( y4 - y3 ), (double) ( x4 - x3 ) ) -
268  atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) );
269  if( da1 >= M_PI )
270  da1 = 2 * M_PI - da1;
271 
272  if( da1 < angle_tolerance )
273  {
274  addSegment( wxPoint( x2, y2 ) );
275  addSegment( wxPoint( x3, y3 ) );
276  return;
277  }
278 
279  if( cusp_limit != 0.0 )
280  {
281  if( da1 > cusp_limit )
282  {
283  addSegment( wxPoint( x3, y3 ) );
284  return;
285  }
286  }
287  }
288  break;
289 
290  case 2:
291 
292  // p1,p3,p4 are collinear, p2 is significant
293  //----------------------
294  if( d2 * d2 <= distance_tolerance_square * sqrt_len(dx, dy) )
295  {
297  {
298  addSegment( wxPoint( x23, y23 ) );
299  return;
300  }
301 
302  // Angle Condition
303  //----------------------
304  da1 = fabs( atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) ) -
305  atan2( (double) ( y2 - y1 ), (double) ( x2 - x1 ) ) );
306  if( da1 >= M_PI )
307  da1 = 2 * M_PI - da1;
308 
309  if( da1 < angle_tolerance )
310  {
311  addSegment( wxPoint( x2, y2 ) );
312  addSegment( wxPoint( x3, y3 ) );
313  return;
314  }
315 
316  if( cusp_limit != 0.0 )
317  {
318  if( da1 > cusp_limit )
319  {
320  addSegment( wxPoint( x2, y2 ) );
321  return;
322  }
323  }
324  }
325  break;
326 
327  case 3:
328 
329  // Regular case
330  //-----------------
331  if( (d2 + d3) * (d2 + d3) <= distance_tolerance_square * sqrt_len(dx, dy) )
332  {
333  // If the curvature doesn't exceed the distance_tolerance value
334  // we tend to finish subdivisions.
335  //----------------------
337  {
338  addSegment( wxPoint( x23, y23 ) );
339  return;
340  }
341 
342  // Angle & Cusp Condition
343  //----------------------
344  k = atan2( (double) ( y3 - y2 ), (double) ( x3 - x2 ) );
345  da1 = fabs( k - atan2( (double) ( y2 - y1 ),
346  (double) ( x2 - x1 ) ) );
347  da2 = fabs( atan2( (double) ( y4 - y3 ),
348  (double) ( x4 - x3 ) ) - k );
349  if( da1 >= M_PI )
350  da1 = 2 * M_PI - da1;
351  if( da2 >= M_PI )
352  da2 = 2 * M_PI - da2;
353 
354  if( da1 + da2 < angle_tolerance )
355  {
356  // Finally we can stop the recursion
357  //----------------------
358  addSegment( wxPoint( x23, y23 ) );
359  return;
360  }
361 
362  if( cusp_limit != 0.0 )
363  {
364  if( da1 > cusp_limit )
365  {
366  addSegment( wxPoint( x2, y2 ) );
367  return;
368  }
369 
370  if( da2 > cusp_limit )
371  {
372  addSegment( wxPoint( x3, y3 ) );
373  return;
374  }
375  }
376  }
377  break;
378  }
379 
380  // Continue subdivision
381  //----------------------
382  recursiveBezier( x1, y1, x12, y12, x123, y123, x1234, y1234, level + 1 );
383  recursiveBezier( x1234, y1234, x234, y234, x34, y34, x4, y4, level + 1 );
384 }
static constexpr double distance_tolerance_square
Definition: bezier_curves.h:89
static constexpr double angle_tolerance
Definition: bezier_curves.h:85
void GetPoly(std::vector< wxPoint > &aOutput)
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:68
void addSegment(const wxPoint &aSegment)
Definition: bezier_curves.h:73
static double sqrt_len(int dx, int dy)
static constexpr double curve_collinearity_epsilon
Definition: bezier_curves.h:91
static constexpr double cusp_limit
Definition: bezier_curves.h:86
std::vector< wxPoint > * m_output
Pointer to the output vector
Definition: bezier_curves.h:71
static constexpr double curve_angle_tolerance_epsilon
Definition: bezier_curves.h:92
static constexpr int recursion_limit
Definition: bezier_curves.h:87
static double calc_sq_distance(int x1, int y1, int x2, int y2)