KiCad PCB EDA Suite
Loading...
Searching...
No Matches
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-2023 KiCad Developers, see AUTHORS.txt for contributors.
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License
8 * as published by the Free Software Foundation; either version 2
9 * of the License, or (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, you may find one here:
18 * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
19 * or you may search the http://www.gnu.org website for the version 2 license,
20 * or you may write to the Free Software Foundation, Inc.,
21 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
22 */
23
24/**********************************************************************************************/
25/* routines to handle bezier curves */
26/* Based on "Fast, Precise Flattening of Cubic Bezier segments offset Curves" by Hain, et. al.*/
27/**********************************************************************************************/
28
29// Portions of this code are based draw2d
30// Copyright (c) 2010, Laurent Le Goff
31// All rights reserved.
32
33// Redistribution and use in source and binary forms, with or without modification,
34// are permitted provided that the following conditions are met:
35
36// * Redistributions of source code must retain the above copyright notice,
37// this list of conditions and the following disclaimer.
38// * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
39// disclaimer in the documentation and/or other materials provided with the distribution.
40
41// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
42// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
43// IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
44// OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
45// OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
46// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
47// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
48
49// Portions of this code are base on BezierKit
50// Copyright (c) 2017 Holmes Futrell
51
52// Permission is hereby granted, free of charge, to any person obtaining a copy
53// of this software and associated documentation files (the "Software"), to deal
54// in the Software without restriction, including without limitation the rights
55// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
56// copies of the Software, and to permit persons to whom the Software is
57// furnished to do so, subject to the following conditions:
58
59// The above copyright notice and this permission notice shall be included in all
60// copies or substantial portions of the Software.
61
62// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
63// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
64// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
65// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
66// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
67// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
68// SOFTWARE.
69
70// Portions of this code are based on the spline-research project
71// Copyright 2018 Raph Levien
72
73// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
74// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
75// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
76// option. This file may not be copied, modified, or distributed
77// except according to those terms.
78
79#include <bezier_curves.h>
80#include <geometry/ellipse.h>
81#include <trigo.h>
82#include <math/vector2d.h> // for VECTOR2D, operator*, VECTOR2
83#include <wx/debug.h> // for wxASSERT
84#include <wx/log.h> // for wxLogTrace
85
86#define BEZIER_DBG "bezier"
87
88BEZIER_POLY::BEZIER_POLY( const VECTOR2I& aStart, const VECTOR2I& aCtrl1,
89 const VECTOR2I& aCtrl2, const VECTOR2I& aEnd )
90{
91 m_ctrlPts.emplace_back( VECTOR2D( aStart ) );
92 m_ctrlPts.emplace_back( VECTOR2D( aCtrl1 ) );
93 m_ctrlPts.emplace_back( VECTOR2D( aCtrl2 ) );
94 m_ctrlPts.emplace_back( VECTOR2D( aEnd ) );
95
96 m_minSegLen = 0.0;
97}
98
99
100BEZIER_POLY::BEZIER_POLY( const std::vector<VECTOR2I>& aControlPoints )
101{
102 m_ctrlPts.reserve( aControlPoints.size() );
103
104 for( const VECTOR2I& pt : aControlPoints )
105 m_ctrlPts.emplace_back( pt );
106
107 m_minSegLen = 0.0;
108}
109
110
112{
113 for( const VECTOR2D& pt : m_ctrlPts )
114 {
115 if( std::isnan( pt.x ) || std::isnan( pt.y ) )
116 return true;
117 }
118
119 return false;
120}
121
122
123bool BEZIER_POLY::isFlat( double aMaxError ) const
124{
125 if( m_ctrlPts.size() == 3 )
126 {
127 VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
128 VECTOR2D D31 = m_ctrlPts[2] - m_ctrlPts[0];
129
130 double t = D21.Dot( D31 ) / D31.SquaredEuclideanNorm();
131 double u = std::min( std::max( t, 0.0 ), 1.0 );
132 VECTOR2D p = m_ctrlPts[0] + u * D31;
133 VECTOR2D delta = m_ctrlPts[1] - p;
134
135 return 0.5 * delta.EuclideanNorm() <= aMaxError;
136 }
137 else if( m_ctrlPts.size() == 4 )
138 {
140
141 VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
142 VECTOR2D D31 = m_ctrlPts[2] - m_ctrlPts[0];
143
144 double cross1 = delta.Cross( D21 );
145 double cross2 = delta.Cross( D31 );
146
147 double inv_delta_sq = 1.0 / delta.SquaredEuclideanNorm();
148 double d1 = ( cross1 * cross1 ) * inv_delta_sq;
149 double d2 = ( cross2 * cross2 ) * inv_delta_sq;
150
151 double factor = ( cross1 * cross2 > 0.0 ) ? 3.0 / 4.0 : 4.0 / 9.0;
152 double f2 = factor * factor;
153 double tol = aMaxError * aMaxError;
154
155 return d1 * f2 <= tol && d2 * f2 <= tol;
156 }
157
158 wxASSERT( false );
159 return false;
160
161}
162
163
164void BEZIER_POLY::GetPoly( std::vector<VECTOR2I>& aOutput, int aMaxError )
165{
166 aOutput.clear();
167 std::vector<VECTOR2D> buffer;
168 GetPoly( buffer, aMaxError );
169
170 aOutput.reserve( buffer.size() );
171
172 for( const VECTOR2D& pt : buffer )
173 aOutput.emplace_back( KiROUND( pt.x ), KiROUND( pt.y ) );
174}
175
176static double approx_int( double x )
177{
178 const double d = 0.6744897501960817;
179 const double d4 = d * d * d * d;
180 return x / ( 1.0 - d + std::pow( d4 + x * x * 0.25, 0.25 ) );
181}
182
183static double approx_inv_int( double x )
184{
185 const double p = 0.39538816;
186 return x * ( 1.0 - p + std::sqrt( p * p + 0.25 * x * x ) );
187}
188
189
191{
192 double omt = 1.0 - t;
193 double omt2 = omt * omt;
194
195 if( m_ctrlPts.size() == 3 )
196 {
197 return omt2 * m_ctrlPts[0] + 2.0 * omt * t * m_ctrlPts[1] + t * t * m_ctrlPts[2];
198 }
199 else if( m_ctrlPts.size() == 4 )
200 {
201 double omt3 = omt * omt2;
202 double t2 = t * t;
203 double t3 = t * t2;
204 return omt3 * m_ctrlPts[0] + 3.0 * t * omt2 * m_ctrlPts[1]
205 + 3.0 * t2 * omt * m_ctrlPts[2] + t3 * m_ctrlPts[3];
206 }
207 else
208 {
209 wxASSERT( false );
210 return VECTOR2D( 0, 0 );
211 }
212}
213
214void BEZIER_POLY::getQuadPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
215{
216 double ddx = 2 * m_ctrlPts[1].x - m_ctrlPts[0].x - m_ctrlPts[2].x;
217 double ddy = 2 * m_ctrlPts[1].y - m_ctrlPts[0].y - m_ctrlPts[2].y;
218 double u0 =
219 ( m_ctrlPts[1].x - m_ctrlPts[0].x ) * ddx + ( m_ctrlPts[1].y - m_ctrlPts[0].y ) * ddy;
220 double u2 =
221 ( m_ctrlPts[2].x - m_ctrlPts[1].x ) * ddx + ( m_ctrlPts[2].y - m_ctrlPts[1].y ) * ddy;
222 double cross =
223 ( m_ctrlPts[2].x - m_ctrlPts[0].x ) * ddy - ( m_ctrlPts[2].y - m_ctrlPts[0].y ) * ddx;
224 double x0 = u0 / cross;
225 double x2 = u2 / cross;
226 double scale = std::abs( cross ) / ( std::hypot( ddx, ddy ) * std::abs( x2 - x0 ) );
227
228 double a0 = approx_int( x0 );
229 double a2 = approx_int( x2 );
230
231 int n = std::ceil( 0.5 * std::abs( a2 - a0 ) * std::sqrt( scale / aMaxError ) );
232
233 double v0 = approx_inv_int( a0 );
234 double v2 = approx_inv_int( a2 );
235
236 aOutput.emplace_back( m_ctrlPts[0] );
237
238 for( int ii = 0; ii < n; ++ii )
239 {
240 double u = approx_inv_int( a0 + ( a2 - a0 ) * ii / n );
241 double t = ( u - v0 ) / ( v2 - v0 );
242 aOutput.emplace_back( eval( t ) );
243 }
244
245 aOutput.emplace_back( m_ctrlPts[2] );
246}
247
248
250{
251 VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
252 VECTOR2D D32 = m_ctrlPts[2] - m_ctrlPts[1];
253 VECTOR2D D43 = m_ctrlPts[3] - m_ctrlPts[2];
254
255 double cross1 = D21.Cross( D32 ) * D32.Cross( D43 );
256 double cross2 = D21.Cross( D32 ) * D21.Cross( D43 );
257
258 if( cross1 < 0.0 )
259 return 1;
260 else if( cross2 > 0.0 )
261 return 0;
262
263 bool b1 = D21.Dot( D32 ) > 0.0;
264 bool b2 = D32.Dot( D43 ) > 0.0;
265
266 if( b1 ^ b2 )
267 return 0;
268
269 wxLogTrace( BEZIER_DBG, "numberOfInflectionPoints: rare case" );
270 // These are rare cases where there are potentially 2 or 0 inflection points.
271 return -1;
272}
273
274
276{
278 double len_sq = delta.SquaredEuclideanNorm();
279
280 if( len_sq < 1e-6 )
281 return 0.0;
282
283 double len = std::sqrt( len_sq );
284 double r = ( m_ctrlPts[1].y - m_ctrlPts[0].y ) / len;
285 double s = ( m_ctrlPts[0].x - m_ctrlPts[1].x ) / len;
286 double u = ( m_ctrlPts[1].x * m_ctrlPts[0].y - m_ctrlPts[0].x * m_ctrlPts[1].y ) / len;
287
288 return std::abs( r * m_ctrlPts[2].x + s * m_ctrlPts[2].y + u );
289}
290
291
292void BEZIER_POLY::subdivide( double aT, BEZIER_POLY& aLeft, BEZIER_POLY& aRight )
293{
294 if( m_ctrlPts.size() == 3 )
295 {
296 aLeft.m_ctrlPts[0] = m_ctrlPts[0];
297 aLeft.m_ctrlPts[1] = m_ctrlPts[0] + aT * ( m_ctrlPts[1] - m_ctrlPts[0] );
298 aLeft.m_ctrlPts[2] = eval( aT );
299
300 aRight.m_ctrlPts[2] = m_ctrlPts[2];
301 aRight.m_ctrlPts[1] = m_ctrlPts[1] + aT * ( m_ctrlPts[2] - m_ctrlPts[1] );
302 aRight.m_ctrlPts[0] = aLeft.m_ctrlPts[2];
303 }
304 else if( m_ctrlPts.size() == 4 )
305 {
306 VECTOR2D left_ctrl1 = m_ctrlPts[0] + aT * ( m_ctrlPts[1] - m_ctrlPts[0] );
307 VECTOR2D tmp = m_ctrlPts[1] + aT * ( m_ctrlPts[2] - m_ctrlPts[1] );
308 VECTOR2D left_ctrl2 = left_ctrl1 + aT * ( tmp - left_ctrl1 );
309 VECTOR2D right_ctrl2 = m_ctrlPts[2] + aT * ( m_ctrlPts[3] - m_ctrlPts[2] );
310 VECTOR2D right_ctrl1 = tmp + aT * ( right_ctrl2 - tmp );
311 VECTOR2D shared = left_ctrl2 + aT * ( right_ctrl1 - left_ctrl2 );
312
313 aLeft.m_ctrlPts[0] = m_ctrlPts[0];
314 aLeft.m_ctrlPts[1] = left_ctrl1;
315 aLeft.m_ctrlPts[2] = left_ctrl2;
316 aLeft.m_ctrlPts[3] = shared;
317
318 aRight.m_ctrlPts[3] = m_ctrlPts[3];
319 aRight.m_ctrlPts[2] = right_ctrl2;
320 aRight.m_ctrlPts[1] = right_ctrl1;
321 aRight.m_ctrlPts[0] = shared;
322 }
323 else
324 {
325 wxASSERT( false );
326 }
327}
328
329
330void BEZIER_POLY::recursiveSegmentation( std::vector<VECTOR2D>& aOutput, double aThreshhold )
331{
332 wxLogTrace( BEZIER_DBG, "recursiveSegmentation with threshold %f", aThreshhold );
333 std::vector<BEZIER_POLY> stack;
334
335 BEZIER_POLY* bezier = nullptr;
336 BEZIER_POLY left( std::vector<VECTOR2D>(4) );
337 BEZIER_POLY right( std::vector<VECTOR2D>(4) );
338
339 stack.push_back( *this );
340
341 while( !stack.empty() )
342 {
343 bezier = &stack.back();
344
345 if( bezier->m_ctrlPts[3] == bezier->m_ctrlPts[0] )
346 {
347 wxLogTrace( BEZIER_DBG, "recursiveSegmentation dropping zero length segment" );
348 stack.pop_back();
349 }
350 else if( bezier->isFlat( aThreshhold ) )
351 {
352 aOutput.push_back( bezier->m_ctrlPts[3] );
353 stack.pop_back();
354 }
355 else
356 {
357 bezier->subdivide( 0.5, left, right );
358 *bezier = right;
359 stack.push_back( left );
360 }
361 }
362
363 wxLogTrace( BEZIER_DBG, "recursiveSegmentation generated %zu points", aOutput.size() );
364}
365
366
367int BEZIER_POLY::findInflectionPoints( double& aT1, double& aT2 )
368{
369 VECTOR2D A{ ( -m_ctrlPts[0].x + 3 * m_ctrlPts[1].x - 3 * m_ctrlPts[2].x + m_ctrlPts[3].x ),
370 ( -m_ctrlPts[0].y + 3 * m_ctrlPts[1].y - 3 * m_ctrlPts[2].y + m_ctrlPts[3].y ) };
371 VECTOR2D B{ ( 3 * m_ctrlPts[0].x - 6 * m_ctrlPts[1].x + 3 * m_ctrlPts[2].x ),
372 ( 3 * m_ctrlPts[0].y - 6 * m_ctrlPts[1].y + 3 * m_ctrlPts[2].y ) };
373 VECTOR2D C{ ( -3 * m_ctrlPts[0].x + 3 * m_ctrlPts[1].x ),
374 ( -3 * m_ctrlPts[0].y + 3 * m_ctrlPts[1].y ) };
375
376 double a = 3 * A.Cross( B );
377 double b = 3 * A.Cross( C );
378 double c = B.Cross( C );
379
380 // Solve the quadratic equation a*t^2 + b*t + c = 0
381 double r2 = ( b * b - 4 * a * c );
382
383 aT1 = 0.0;
384 aT2 = 0.0;
385
386 if( r2 >= 0.0 && a != 0.0 )
387 {
388 double r = std::sqrt( r2 );
389
390 double t1 = ( ( -b + r ) / ( 2 * a ) );
391 double t2 = ( ( -b - r ) / ( 2 * a ) );
392
393 if( ( t1 > 0.0 && t1 < 1.0 ) && ( t2 > 0.0 && t2 < 1.0 ) )
394 {
395 if( t1 > t2 )
396 {
397 std::swap( t1, t2 );
398 }
399
400 aT1 = t1;
401 aT2 = t2;
402
403 if( t2 - t1 > 0.00001 )
404 {
405 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 2 inflection points at t1 = %f, t2 = %f", t1, t2 );
406 return 2;
407 }
408 else
409 {
410 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t1 );
411 return 1;
412 }
413 }
414 else if( t1 > 0.0 && t1 < 1.0 )
415 {
416 aT1 = t1;
417 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t1 );
418 return 1;
419 }
420 else if( t2 > 0.0 && t2 < 1.0 )
421 {
422 aT1 = t2;
423 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t2 );
424 return 1;
425 }
426
427 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found no inflection points" );
428 return 0;
429 }
430
431 wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found no inflection points" );
432 return 0;
433}
434
435
436void BEZIER_POLY::cubicParabolicApprox( std::vector<VECTOR2D>& aOutput, double aMaxError )
437{
438 std::vector<BEZIER_POLY> stack;
439 stack.push_back( std::vector<VECTOR2D>(4) );
440 stack.push_back( std::vector<VECTOR2D>(4) );
441 stack.push_back( std::vector<VECTOR2D>(4) );
442 stack.push_back( std::vector<VECTOR2D>(4) );
443
444 BEZIER_POLY* c = this;
445 BEZIER_POLY* b1 = &stack[0];
446 BEZIER_POLY* b2 = &stack[1];
447
448 for( ;; )
449 {
450 if( c->isNaN() )
451 {
452 wxLogDebug( "cubicParabolicApprox: NaN detected" );
453 break;
454 }
455
456 if( c->isFlat( aMaxError ) )
457 {
458 wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: General Flatness detected, adding %f %f", c->m_ctrlPts[3].x, c->m_ctrlPts[3].y );
459 // If the subsegment deviation satisfies the flatness criterion, store the last point and stop
460 aOutput.push_back( c->m_ctrlPts[3] );
461 break;
462 }
463
464 // Find the third control point deviation and the t values for subdivision
465 double d = c->thirdControlPointDeviation();
466 double t = 2 * std::sqrt( aMaxError / ( 3.0 * d ) ); // Forumla 2 in Hain et al.
467
468 wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: split point t = %f", t );
469
470 if( t > 1.0 )
471 {
472 wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: Invalid t value detected" );
473 // Case where the t value calculated is invalid, so use recursive subdivision
474 c->recursiveSegmentation( aOutput, aMaxError );
475 break;
476 }
477
478 // Valid t value to subdivide at that calculated value
479 c->subdivide( t, *b1, *b2 );
480
481 // First subsegment should have its deviation equal to flatness
482 if( b1->isFlat( aMaxError ) )
483 {
484 wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: Flatness detected, adding %f %f", b1->m_ctrlPts[3].x, b1->m_ctrlPts[3].y );
485 aOutput.push_back( b1->m_ctrlPts[3] );
486 }
487 else
488 {
489 // if not then use segment to handle any mathematical errors
490 b1->recursiveSegmentation( aOutput, aMaxError );
491 }
492
493 // Repeat the process for the left over subsegment
494 c = b2;
495
496 if( b1 == &stack.front() )
497 {
498 b1 = &stack[2];
499 b2 = &stack[3];
500 }
501 else
502 {
503 b1 = &stack[0];
504 b2 = &stack[1];
505 }
506 }
507}
508
509
510void BEZIER_POLY::getCubicPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
511{
512 aOutput.push_back( m_ctrlPts[0] );
513
514 if( numberOfInflectionPoints() == 0 )
515 {
516 wxLogTrace( BEZIER_DBG, "getCubicPoly Short circuit to PA" );
517 // If no inflection points then apply PA on the full Bezier segment.
518 cubicParabolicApprox( aOutput, aMaxError );
519 return;
520 }
521
522 // If one or more inflection points then we will have to subdivide the curve
523 double t1, t2;
524 int numOfIfP = findInflectionPoints( t1, t2 );
525
526 if( numOfIfP == 2 )
527 {
528 wxLogTrace( BEZIER_DBG, "getCubicPoly: 2 inflection points" );
529 // Case when 2 inflection points then divide at the smallest one first
530 BEZIER_POLY sub1( std::vector<VECTOR2D>( 4 ) );
531 BEZIER_POLY tmp1( std::vector<VECTOR2D>( 4 ) );
532 BEZIER_POLY sub2( std::vector<VECTOR2D>( 4 ) );
533 BEZIER_POLY sub3( std::vector<VECTOR2D>( 4 ) );
534
535 subdivide( t1, sub1, tmp1 );
536
537 // Now find the second inflection point in the second curve and subdivide
538 numOfIfP = tmp1.findInflectionPoints( t1, t2 );
539 if( numOfIfP == 2 )
540 tmp1.subdivide( t1, sub2, sub3 );
541 else if( numOfIfP == 1 )
542 tmp1.subdivide( t1, sub2, sub3 );
543 else
544 {
545 wxLogTrace( BEZIER_DBG, "getCubicPoly: 2nd inflection point not found" );
546 return;
547 }
548
549 // Use PA for first subsegment
550 sub1.cubicParabolicApprox( aOutput, aMaxError );
551
552 // Use Segment for the second (middle) subsegment
553 sub2.recursiveSegmentation( aOutput, aMaxError );
554
555 // Use PA for the third curve
556 sub3.cubicParabolicApprox( aOutput, aMaxError );
557 }
558 else if( numOfIfP == 1 )
559 {
560 wxLogTrace( BEZIER_DBG, "getCubicPoly: 1 inflection point" );
561 // Case where there is one inflection point, subdivide once and use PA on both subsegments
562 BEZIER_POLY sub1( std::vector<VECTOR2D>( 4 ) );
563 BEZIER_POLY sub2( std::vector<VECTOR2D>( 4 ) );
564 subdivide( t1, sub1, sub2 );
565 sub1.cubicParabolicApprox( aOutput, aMaxError );
566 sub2.cubicParabolicApprox( aOutput, aMaxError );
567 }
568 else
569 {
570 wxLogTrace( BEZIER_DBG, "getCubicPoly: Unknown inflection points" );
571 // Case where there is no inflection, use PA directly
572 cubicParabolicApprox( aOutput, aMaxError );
573 }
574}
575
576
577void BEZIER_POLY::GetPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
578{
579 if( aMaxError <= 0.0 )
580 aMaxError = 10.0;
581
582 if( m_ctrlPts.size() == 3 )
583 {
584 getQuadPoly( aOutput, aMaxError );
585 }
586 else if( m_ctrlPts.size() == 4 )
587 {
588 getCubicPoly( aOutput, aMaxError );
589 }
590 else
591 {
592 wxASSERT( false );
593 }
594
595 wxLogTrace( BEZIER_DBG, "GetPoly generated %zu points", aOutput.size() );
596}
597
598
599template<typename T>
600void TransformEllipseToBeziers( const ELLIPSE<T>& aEllipse, std::vector<BEZIER<T>>& aBeziers )
601{
602 EDA_ANGLE arcAngle = -( aEllipse.EndAngle - aEllipse.StartAngle );
603
604 if( arcAngle >= ANGLE_0 )
605 arcAngle -= ANGLE_360;
606
607 /*
608 * KiCad does not natively support ellipses or elliptical arcs. So, we convert them to Bezier
609 * splines as these are the nearest thing we have that represents them in a way that is both
610 * editable and preserves their curvature accurately (enough).
611 *
612 * Credit to Kliment for developing and documenting this method.
613 */
615 const int minBeziersPerCircle = 4;
616
618 const int numBeziers = static_cast<int>(
619 std::ceil( std::abs( arcAngle.AsRadians() / ( 2 * M_PI / minBeziersPerCircle ) ) ) );
620
622 const double angleIncrement = arcAngle.AsRadians() / numBeziers;
623
624 /*
625 * Now, let's assume a circle of radius 1, centered on origin, with angle startangle
626 * x-axis-aligned. We'll move, scale, and rotate it later. We're creating Bezier curves that hug
627 * this circle as closely as possible, with the angles that will be used on the final ellipse
628 * too.
629 *
630 * Thanks to the beautiful and excellent https://pomax.github.io/bezierinfo we know how to
631 * define a curve that hugs a circle as closely as possible.
632 *
633 * We need the value k, which is the optimal distance from the endpoint to the control point to
634 * make the curve match the circle for a given circle arc angle.
635 *
636 * k = 4/3 * tan(θ/4), where θ is the angle of the arc. In our case, θ=angleIncrement
637 */
638 double theta = angleIncrement;
639 double k = ( 4. / 3. ) * std::tan( theta / 4 );
640
641 /*
642 * Define our Bezier:
643 * - Start point is on the circle at the x-axis
644 * - First control point just uses k as the y-value
645 * - Second control point is offset from the end point
646 * - End point is defined by the angle of the arc segment
647 * Note that we use double here no matter what the template param is; round at the end only.
648 */
649 BEZIER<double> first = { { 1, 0 },
650 { 1, k },
651 { std::cos( theta ) + k * std::sin( theta ),
652 std::sin( theta ) - k * std::cos( theta ) },
653 { std::cos( theta ), std::sin( theta ) } };
654
655 /*
656 * Now construct the actual segments by transforming/rotating the first one
657 */
658 auto transformPoint =
659 [&]( VECTOR2D aPoint, const double aAngle ) -> VECTOR2D
660 {
661 // Bring to the actual starting angle
662 RotatePoint( aPoint,
663 -EDA_ANGLE( aAngle - aEllipse.StartAngle.AsRadians(), RADIANS_T ) );
664
665 // Then scale to the major and minor radiuses of the ellipse
666 aPoint *= VECTOR2D( aEllipse.MajorRadius, aEllipse.MinorRadius );
667
668 // Now rotate to the ellipse coordinate system
669 RotatePoint( aPoint, -aEllipse.Rotation );
670
671 // And finally offset to the center location of the ellipse
672 aPoint += aEllipse.Center;
673
674 return aPoint;
675 };
676
677 for( int i = 0; i < numBeziers; i++ )
678 {
679 aBeziers.emplace_back( BEZIER<T>( {
680 transformPoint( first.Start, i * angleIncrement ),
681 transformPoint( first.C1, i * angleIncrement ),
682 transformPoint( first.C2, i * angleIncrement ),
683 transformPoint( first.End, i * angleIncrement )
684 } ) );
685 }
686}
687
688
689template void TransformEllipseToBeziers( const ELLIPSE<double>& aEllipse,
690 std::vector<BEZIER<double>>& aBeziers );
691template void TransformEllipseToBeziers( const ELLIPSE<int>& aEllipse,
692 std::vector<BEZIER<int>>& aBeziers );
static double approx_int(double x)
static double approx_inv_int(double x)
#define BEZIER_DBG
void TransformEllipseToBeziers(const ELLIPSE< T > &aEllipse, std::vector< BEZIER< T > > &aBeziers)
Transforms an ellipse or elliptical arc into a set of quadratic Bezier curves that approximate it.
Bezier curves to polygon converter.
Definition: bezier_curves.h:38
void recursiveSegmentation(std::vector< VECTOR2D > &aOutput, double aMaxError)
void getQuadPoly(std::vector< VECTOR2D > &aOutput, double aMaxError)
int numberOfInflectionPoints()
double thirdControlPointDeviation()
void cubicParabolicApprox(std::vector< VECTOR2D > &aOutput, double aMaxError)
double m_minSegLen
Control points.
Definition: bezier_curves.h:81
void getCubicPoly(std::vector< VECTOR2D > &aOutput, double aMaxError)
BEZIER_POLY(const VECTOR2I &aStart, const VECTOR2I &aCtrl1, const VECTOR2I &aCtrl2, const VECTOR2I &aEnd)
void subdivide(double aT, BEZIER_POLY &aLeft, BEZIER_POLY &aRight)
bool isFlat(double aMaxError) const
void GetPoly(std::vector< VECTOR2I > &aOutput, int aMaxError=10)
Convert a Bezier curve to a polygon.
int findInflectionPoints(double &aT1, double &aT2)
bool isNaN() const
VECTOR2D eval(double t)
std::vector< VECTOR2D > m_ctrlPts
Definition: bezier_curves.h:84
Generic cubic Bezier representation.
Definition: bezier_curves.h:95
VECTOR2< NumericType > Start
VECTOR2< NumericType > C1
VECTOR2< NumericType > C2
VECTOR2< NumericType > End
double AsRadians() const
Definition: eda_angle.h:117
This class was created to handle importing ellipses from other file formats that support them nativel...
Definition: ellipse.h:34
NumericType MinorRadius
Definition: ellipse.h:69
EDA_ANGLE Rotation
Definition: ellipse.h:70
EDA_ANGLE StartAngle
Definition: ellipse.h:71
NumericType MajorRadius
Definition: ellipse.h:68
EDA_ANGLE EndAngle
Definition: ellipse.h:72
VECTOR2< NumericType > Center
Definition: ellipse.h:67
extended_type SquaredEuclideanNorm() const
Compute the squared euclidean norm of the vector, which is defined as (x ** 2 + y ** 2).
Definition: vector2d.h:302
extended_type Cross(const VECTOR2< T > &aVector) const
Compute cross product of self with aVector.
Definition: vector2d.h:535
extended_type Dot(const VECTOR2< T > &aVector) const
Compute dot product of self with aVector.
Definition: vector2d.h:543
static constexpr EDA_ANGLE ANGLE_0
Definition: eda_angle.h:401
@ RADIANS_T
Definition: eda_angle.h:32
static constexpr EDA_ANGLE ANGLE_360
Definition: eda_angle.h:407
EDA_ANGLE abs(const EDA_ANGLE &aAngle)
Definition: eda_angle.h:390
const int scale
VECTOR2I v2(1, 0)
Test suite for KiCad math code.
constexpr int delta
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
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:121
VECTOR2< double > VECTOR2D
Definition: vector2d.h:672