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