bezier_curves.cpp

Go to the documentation of this file.

/*
 * This program source code file is part of KiCad, a free EDA CAD application.
 *
 * Copyright (C) 2014-2023 KiCad Developers, see AUTHORS.txt for contributors.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, you may find one here:
 * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
 * or you may search the http://www.gnu.org website for the version 2 license,
 * or you may write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
 */
 
/**********************************************************************************************/
/* routines to handle bezier curves                                                           */
/* Based on "Fast, Precise Flattening of Cubic Bezier segments offset Curves" by Hain, et. al.*/
/**********************************************************************************************/
 
// Portions of this code are based draw2d
// Copyright (c) 2010, Laurent Le Goff
// All rights reserved.
 
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
 
//  * Redistributions of source code must retain the above copyright notice,
//      this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
//          disclaimer in the documentation and/or other materials provided with the distribution.
 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
// OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
// OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
// Portions of this code are base on BezierKit
// Copyright (c) 2017 Holmes Futrell
 
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
 
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
 
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
 
// Portions of this code are based on the spline-research project
// Copyright 2018 Raph Levien
 
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.
 
#include <bezier_curves.h>
#include <geometry/ellipse.h>
#include <trigo.h>
#include <math/vector2d.h>  // for VECTOR2D, operator*, VECTOR2
#include <wx/debug.h>       // for wxASSERT
#include <wx/log.h>         // for wxLogTrace
 
#define BEZIER_DBG "bezier"
 
BEZIER_POLY::BEZIER_POLY( const VECTOR2I& aStart, const VECTOR2I& aCtrl1,
                          const VECTOR2I& aCtrl2, const VECTOR2I& aEnd )
{
    m_ctrlPts.emplace_back( VECTOR2D( aStart ) );
    m_ctrlPts.emplace_back( VECTOR2D( aCtrl1 ) );
    m_ctrlPts.emplace_back( VECTOR2D( aCtrl2 ) );
    m_ctrlPts.emplace_back( VECTOR2D( aEnd ) );
 
    m_minSegLen = 0.0;
}
 
 
BEZIER_POLY::BEZIER_POLY( const std::vector<VECTOR2I>& aControlPoints )
{
    m_ctrlPts.reserve( aControlPoints.size() );
 
    for( const VECTOR2I& pt : aControlPoints )
        m_ctrlPts.emplace_back( pt );
 
    m_minSegLen = 0.0;
}
 
 
bool BEZIER_POLY::isNaN() const
{
    for( const VECTOR2D& pt : m_ctrlPts )
    {
        if( std::isnan( pt.x ) || std::isnan( pt.y ) )
            return true;
    }
 
    return false;
}
 
 
bool BEZIER_POLY::isFlat( double aMaxError ) const
{
    if( m_ctrlPts.size() == 3 )
    {
        VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
        VECTOR2D D31 = m_ctrlPts[2] - m_ctrlPts[0];
 
        double   t = D21.Dot( D31 ) / D31.SquaredEuclideanNorm();
        double   u = std::min( std::max( t, 0.0 ), 1.0 );
        VECTOR2D p = m_ctrlPts[0] + u * D31;
        VECTOR2D delta = m_ctrlPts[1] - p;
 
        return 0.5 * delta.EuclideanNorm() <= aMaxError;
    }
    else if( m_ctrlPts.size() == 4 )
    {
        VECTOR2D delta = m_ctrlPts[3] - m_ctrlPts[0];
 
        VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
        VECTOR2D D31 = m_ctrlPts[2] - m_ctrlPts[0];
 
        double cross1 = delta.Cross( D21 );
        double cross2 = delta.Cross( D31 );
 
        double inv_delta_sq = 1.0 / delta.SquaredEuclideanNorm();
        double d1 = ( cross1 * cross1 ) * inv_delta_sq;
        double d2 = ( cross2 * cross2 ) * inv_delta_sq;
 
        double factor = ( cross1 * cross2 > 0.0 ) ? 3.0 / 4.0 : 4.0 / 9.0;
        double f2 = factor * factor;
        double tol = aMaxError * aMaxError;
 
        return d1 * f2 <= tol && d2 * f2 <= tol;
    }
 
    wxASSERT( false );
    return false;
 
}
 
 
void BEZIER_POLY::GetPoly( std::vector<VECTOR2I>& aOutput, int aMaxError )
{
    aOutput.clear();
    std::vector<VECTOR2D> buffer;
    GetPoly( buffer, aMaxError );
 
    aOutput.reserve( buffer.size() );
 
    for( const VECTOR2D& pt : buffer )
        aOutput.emplace_back( KiROUND( pt.x ), KiROUND( pt.y ) );
}
 
static double approx_int( double x )
{
    const double d = 0.6744897501960817;
    const double d4 = d * d * d * d;
    return x / ( 1.0 - d + std::pow( d4 + x * x * 0.25, 0.25 ) );
}
 
static double approx_inv_int( double x )
{
    const double p = 0.39538816;
    return x * ( 1.0 - p + std::sqrt( p * p + 0.25 * x * x ) );
}
 
 
VECTOR2D BEZIER_POLY::eval( double t )
{
    double omt = 1.0 - t;
    double omt2 = omt * omt;
 
    if( m_ctrlPts.size() == 3 )
    {
        return omt2 * m_ctrlPts[0] + 2.0 * omt * t * m_ctrlPts[1] + t * t * m_ctrlPts[2];
    }
    else if( m_ctrlPts.size() == 4 )
    {
        double omt3 = omt * omt2;
        double t2 = t * t;
        double t3 = t * t2;
        return omt3 * m_ctrlPts[0] + 3.0 * t * omt2 * m_ctrlPts[1]
               + 3.0 * t2 * omt * m_ctrlPts[2] + t3 * m_ctrlPts[3];
    }
    else
    {
        wxASSERT( false );
        return VECTOR2D( 0, 0 );
    }
}
 
void BEZIER_POLY::getQuadPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
{
    double ddx = 2 * m_ctrlPts[1].x - m_ctrlPts[0].x - m_ctrlPts[2].x;
    double ddy = 2 * m_ctrlPts[1].y - m_ctrlPts[0].y - m_ctrlPts[2].y;
    double u0 =
            ( m_ctrlPts[1].x - m_ctrlPts[0].x ) * ddx + ( m_ctrlPts[1].y - m_ctrlPts[0].y ) * ddy;
    double u2 =
            ( m_ctrlPts[2].x - m_ctrlPts[1].x ) * ddx + ( m_ctrlPts[2].y - m_ctrlPts[1].y ) * ddy;
    double cross =
            ( m_ctrlPts[2].x - m_ctrlPts[0].x ) * ddy - ( m_ctrlPts[2].y - m_ctrlPts[0].y ) * ddx;
    double x0 = u0 / cross;
    double x2 = u2 / cross;
    double scale = std::abs( cross ) / ( std::hypot( ddx, ddy ) * std::abs( x2 - x0 ) );
 
    double a0 = approx_int( x0 );
    double a2 = approx_int( x2 );
 
    int n = std::ceil( 0.5 * std::abs( a2 - a0 ) * std::sqrt( scale / aMaxError ) );
 
    double v0 = approx_inv_int( a0 );
    double v2 = approx_inv_int( a2 );
 
    aOutput.emplace_back( m_ctrlPts[0] );
 
    for( int ii = 0; ii < n; ++ii )
    {
        double u = approx_inv_int( a0 + ( a2 - a0 ) * ii / n );
        double t = ( u - v0 ) / ( v2 - v0 );
        aOutput.emplace_back( eval( t ) );
    }
 
    aOutput.emplace_back( m_ctrlPts[2] );
}
 
 
int BEZIER_POLY::numberOfInflectionPoints()
{
    VECTOR2D D21 = m_ctrlPts[1] - m_ctrlPts[0];
    VECTOR2D D32 = m_ctrlPts[2] - m_ctrlPts[1];
    VECTOR2D D43 = m_ctrlPts[3] - m_ctrlPts[2];
 
    double cross1 = D21.Cross( D32 ) * D32.Cross( D43 );
    double cross2 = D21.Cross( D32 ) * D21.Cross( D43 );
 
    if( cross1 < 0.0 )
        return 1;
    else if( cross2 > 0.0 )
        return 0;
 
    bool b1 = D21.Dot( D32 ) > 0.0;
    bool b2 = D32.Dot( D43 ) > 0.0;
 
    if( b1 ^ b2 )
        return 0;
 
    wxLogTrace( BEZIER_DBG, "numberOfInflectionPoints: rare case" );
    // These are rare cases where there are potentially 2 or 0 inflection points.
    return -1;
}
 
 
double BEZIER_POLY::thirdControlPointDeviation()
{
    VECTOR2D delta = m_ctrlPts[1] - m_ctrlPts[0];
    double len_sq = delta.SquaredEuclideanNorm();
 
    if( len_sq < 1e-6 )
        return 0.0;
 
    double len = std::sqrt( len_sq );
    double r = ( m_ctrlPts[1].y - m_ctrlPts[0].y ) / len;
    double s = ( m_ctrlPts[0].x - m_ctrlPts[1].x ) / len;
    double u = ( m_ctrlPts[1].x * m_ctrlPts[0].y - m_ctrlPts[0].x * m_ctrlPts[1].y ) / len;
 
    return std::abs( r * m_ctrlPts[2].x + s * m_ctrlPts[2].y + u );
}
 
 
void BEZIER_POLY::subdivide( double aT, BEZIER_POLY& aLeft, BEZIER_POLY& aRight )
{
    if( m_ctrlPts.size() == 3 )
    {
        aLeft.m_ctrlPts[0] = m_ctrlPts[0];
        aLeft.m_ctrlPts[1] = m_ctrlPts[0] + aT * ( m_ctrlPts[1] - m_ctrlPts[0] );
        aLeft.m_ctrlPts[2] = eval( aT );
 
        aRight.m_ctrlPts[2] = m_ctrlPts[2];
        aRight.m_ctrlPts[1] = m_ctrlPts[1] + aT * ( m_ctrlPts[2] - m_ctrlPts[1] );
        aRight.m_ctrlPts[0] = aLeft.m_ctrlPts[2];
    }
    else if( m_ctrlPts.size() == 4 )
    {
        VECTOR2D left_ctrl1 = m_ctrlPts[0] + aT * ( m_ctrlPts[1] - m_ctrlPts[0] );
        VECTOR2D tmp = m_ctrlPts[1] + aT * ( m_ctrlPts[2] - m_ctrlPts[1] );
        VECTOR2D left_ctrl2 = left_ctrl1 + aT * ( tmp - left_ctrl1 );
        VECTOR2D right_ctrl2 = m_ctrlPts[2] + aT * ( m_ctrlPts[3] - m_ctrlPts[2] );
        VECTOR2D right_ctrl1 = tmp + aT * ( right_ctrl2 - tmp );
        VECTOR2D shared = left_ctrl2 + aT * ( right_ctrl1 - left_ctrl2 );
 
        aLeft.m_ctrlPts[0] = m_ctrlPts[0];
        aLeft.m_ctrlPts[1] = left_ctrl1;
        aLeft.m_ctrlPts[2] = left_ctrl2;
        aLeft.m_ctrlPts[3] = shared;
 
        aRight.m_ctrlPts[3] = m_ctrlPts[3];
        aRight.m_ctrlPts[2] = right_ctrl2;
        aRight.m_ctrlPts[1] = right_ctrl1;
        aRight.m_ctrlPts[0] = shared;
    }
    else
    {
        wxASSERT( false );
    }
}
 
 
void BEZIER_POLY::recursiveSegmentation( std::vector<VECTOR2D>& aOutput, double aThreshhold )
{
    wxLogTrace( BEZIER_DBG, "recursiveSegmentation with threshold %f", aThreshhold );
    std::vector<BEZIER_POLY> stack;
 
    BEZIER_POLY* bezier = nullptr;
    BEZIER_POLY left( std::vector<VECTOR2D>(4) );
    BEZIER_POLY right( std::vector<VECTOR2D>(4) );
 
    stack.push_back( *this );
 
    while( !stack.empty() )
    {
        bezier = &stack.back();
 
        if( bezier->m_ctrlPts[3] == bezier->m_ctrlPts[0] )
        {
            wxLogTrace( BEZIER_DBG, "recursiveSegmentation dropping zero length segment" );
            stack.pop_back();
        }
        else if( bezier->isFlat( aThreshhold ) )
        {
            aOutput.push_back( bezier->m_ctrlPts[3] );
            stack.pop_back();
        }
        else
        {
            bezier->subdivide( 0.5, left, right );
            *bezier = right;
            stack.push_back( left );
        }
    }
 
    wxLogTrace( BEZIER_DBG, "recursiveSegmentation generated %zu points", aOutput.size() );
}
 
 
int BEZIER_POLY::findInflectionPoints( double& aT1, double& aT2 )
{
    VECTOR2D A{ ( -m_ctrlPts[0].x + 3 * m_ctrlPts[1].x - 3 * m_ctrlPts[2].x + m_ctrlPts[3].x ),
                ( -m_ctrlPts[0].y + 3 * m_ctrlPts[1].y - 3 * m_ctrlPts[2].y + m_ctrlPts[3].y ) };
    VECTOR2D B{ ( 3 * m_ctrlPts[0].x - 6 * m_ctrlPts[1].x + 3 * m_ctrlPts[2].x ),
                ( 3 * m_ctrlPts[0].y - 6 * m_ctrlPts[1].y + 3 * m_ctrlPts[2].y ) };
    VECTOR2D C{ ( -3 * m_ctrlPts[0].x + 3 * m_ctrlPts[1].x ),
                ( -3 * m_ctrlPts[0].y + 3 * m_ctrlPts[1].y ) };
 
    double a = 3 * A.Cross( B );
    double b = 3 * A.Cross( C );
    double c = B.Cross( C );
 
    // Solve the quadratic equation a*t^2 + b*t + c = 0
    double r2 = ( b * b - 4 * a * c );
 
    aT1 = 0.0;
    aT2 = 0.0;
 
    if( r2 >= 0.0 && a != 0.0 )
    {
        double r = std::sqrt( r2 );
 
        double t1 = ( ( -b + r ) / ( 2 * a ) );
        double t2 = ( ( -b - r ) / ( 2 * a ) );
 
        if( ( t1 > 0.0 && t1 < 1.0 ) && ( t2 > 0.0 && t2 < 1.0 ) )
        {
            if( t1 > t2 )
            {
                std::swap( t1, t2 );
            }
 
            aT1 = t1;
            aT2 = t2;
 
            if( t2 - t1 > 0.00001 )
            {
                wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 2 inflection points at t1 = %f, t2 = %f", t1, t2 );
                return 2;
            }
            else
            {
                wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t1 );
                return 1;
            }
        }
        else if( t1 > 0.0 && t1 < 1.0 )
        {
            aT1 = t1;
            wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t1 );
            return 1;
        }
        else if( t2 > 0.0 && t2 < 1.0 )
        {
            aT1 = t2;
            wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found 1 inflection point at t = %f", t2 );
            return 1;
        }
 
        wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found no inflection points" );
        return 0;
    }
 
    wxLogTrace( BEZIER_DBG, "BEZIER_POLY Found no inflection points" );
    return 0;
}
 
 
void BEZIER_POLY::cubicParabolicApprox( std::vector<VECTOR2D>& aOutput, double aMaxError )
{
    std::vector<BEZIER_POLY> stack;
    stack.push_back( std::vector<VECTOR2D>(4) );
    stack.push_back( std::vector<VECTOR2D>(4) );
    stack.push_back( std::vector<VECTOR2D>(4) );
    stack.push_back( std::vector<VECTOR2D>(4) );
 
    BEZIER_POLY* c = this;
    BEZIER_POLY* b1 = &stack[0];
    BEZIER_POLY* b2 = &stack[1];
 
    for( ;; )
    {
        if( c->isNaN() )
        {
            wxLogDebug( "cubicParabolicApprox: NaN detected" );
            break;
        }
 
        if( c->isFlat( aMaxError ) )
        {
            wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: General Flatness detected, adding %f %f", c->m_ctrlPts[3].x, c->m_ctrlPts[3].y );
            // If the subsegment deviation satisfies the flatness criterion, store the last point and stop
            aOutput.push_back( c->m_ctrlPts[3] );
            break;
        }
 
        // Find the third control point deviation and the t values for subdivision
        double d = c->thirdControlPointDeviation();
        double t = 2 * std::sqrt( aMaxError / ( 3.0 * d ) ); // Forumla 2 in Hain et al.
 
        wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: split point t = %f", t );
 
        if( t > 1.0 )
        {
            wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: Invalid t value detected" );
            // Case where the t value calculated is invalid, so use recursive subdivision
            c->recursiveSegmentation( aOutput, aMaxError );
            break;
        }
 
        // Valid t value to subdivide at that calculated value
        c->subdivide( t, *b1, *b2 );
 
        // First subsegment should have its deviation equal to flatness
        if( b1->isFlat( aMaxError ) )
        {
            wxLogTrace( BEZIER_DBG, "cubicParabolicApprox: Flatness detected, adding %f %f", b1->m_ctrlPts[3].x, b1->m_ctrlPts[3].y );
            aOutput.push_back( b1->m_ctrlPts[3] );
        }
        else
        {
            // if not then use segment to handle any mathematical errors
            b1->recursiveSegmentation( aOutput, aMaxError );
        }
 
        // Repeat the process for the left over subsegment
        c = b2;
 
        if( b1 == &stack.front() )
        {
            b1 = &stack[2];
            b2 = &stack[3];
        }
        else
        {
            b1 = &stack[0];
            b2 = &stack[1];
        }
    }
}
 
 
void BEZIER_POLY::getCubicPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
{
    aOutput.push_back( m_ctrlPts[0] );
 
    if( numberOfInflectionPoints() == 0 )
    {
        wxLogTrace( BEZIER_DBG, "getCubicPoly Short circuit to PA" );
        // If no inflection points then apply PA on the full Bezier segment.
        cubicParabolicApprox( aOutput, aMaxError );
        return;
    }
 
    // If one or more inflection points then we will have to subdivide the curve
    double t1, t2;
    int    numOfIfP = findInflectionPoints( t1, t2 );
 
    if( numOfIfP == 2 )
    {
        wxLogTrace( BEZIER_DBG, "getCubicPoly: 2 inflection points" );
        // Case when 2 inflection points then divide at the smallest one first
        BEZIER_POLY sub1( std::vector<VECTOR2D>( 4 ) );
        BEZIER_POLY tmp1( std::vector<VECTOR2D>( 4 ) );
        BEZIER_POLY sub2( std::vector<VECTOR2D>( 4 ) );
        BEZIER_POLY sub3( std::vector<VECTOR2D>( 4 ) );
 
        subdivide( t1, sub1, tmp1 );
 
        // Now find the second inflection point in the second curve and subdivide
        numOfIfP = tmp1.findInflectionPoints( t1, t2 );
        if( numOfIfP == 2 )
            tmp1.subdivide( t1, sub2, sub3 );
        else if( numOfIfP == 1 )
            tmp1.subdivide( t1, sub2, sub3 );
        else
        {
            wxLogTrace( BEZIER_DBG, "getCubicPoly: 2nd inflection point not found" );
            return;
        }
 
        // Use PA for first subsegment
        sub1.cubicParabolicApprox( aOutput, aMaxError );
 
        // Use Segment for the second (middle) subsegment
        sub2.recursiveSegmentation( aOutput, aMaxError );
 
        // Use PA for the third curve
        sub3.cubicParabolicApprox( aOutput, aMaxError );
    }
    else if( numOfIfP == 1 )
    {
        wxLogTrace( BEZIER_DBG, "getCubicPoly: 1 inflection point" );
        // Case where there is one inflection point, subdivide once and use PA on both subsegments
        BEZIER_POLY sub1( std::vector<VECTOR2D>( 4 ) );
        BEZIER_POLY sub2( std::vector<VECTOR2D>( 4 ) );
        subdivide( t1, sub1, sub2 );
        sub1.cubicParabolicApprox( aOutput, aMaxError );
        sub2.cubicParabolicApprox( aOutput, aMaxError );
    }
    else
    {
        wxLogTrace( BEZIER_DBG, "getCubicPoly: Unknown inflection points" );
        // Case where there is no inflection, use PA directly
        cubicParabolicApprox( aOutput, aMaxError );
    }
}
 
 
void BEZIER_POLY::GetPoly( std::vector<VECTOR2D>& aOutput, double aMaxError )
{
    if( aMaxError <= 0.0 )
        aMaxError = 10.0;
 
    if( m_ctrlPts.size() == 3 )
    {
        getQuadPoly( aOutput, aMaxError );
    }
    else if( m_ctrlPts.size() == 4 )
    {
        getCubicPoly( aOutput, aMaxError );
    }
    else
    {
        wxASSERT( false );
    }
 
    wxLogTrace( BEZIER_DBG, "GetPoly generated %zu points", aOutput.size() );
}
 
 
template<typename T>
void TransformEllipseToBeziers( const ELLIPSE<T>& aEllipse, std::vector<BEZIER<T>>& aBeziers )
{
    EDA_ANGLE arcAngle = -( aEllipse.EndAngle - aEllipse.StartAngle );
 
    if( arcAngle >= ANGLE_0 )
        arcAngle -= ANGLE_360;
 
    /*
     * KiCad does not natively support ellipses or elliptical arcs.  So, we convert them to Bezier
     * splines as these are the nearest thing we have that represents them in a way that is both
     * editable and preserves their curvature accurately (enough).
     *
     * Credit to Kliment for developing and documenting this method.
     */
    const int minBeziersPerCircle = 4;
 
    const int numBeziers = static_cast<int>(
            std::ceil( std::abs( arcAngle.AsRadians() / ( 2 * M_PI / minBeziersPerCircle ) ) ) );
 
    const double angleIncrement = arcAngle.AsRadians() / numBeziers;
 
    /*
     * Now, let's assume a circle of radius 1, centered on origin, with angle startangle
     * x-axis-aligned. We'll move, scale, and rotate it later. We're creating Bezier curves that hug
     * this circle as closely as possible, with the angles that will be used on the final ellipse
     * too.
     *
     * Thanks to the beautiful and excellent https://pomax.github.io/bezierinfo we know how to
     * define a curve that hugs a circle as closely as possible.
     *
     * We need the value k, which is the optimal distance from the endpoint to the control point to
     * make the curve match the circle for a given circle arc angle.
     *
     * k = 4/3 * tan(θ/4), where θ is the angle of the arc. In our case, θ=angleIncrement
     */
    double theta = angleIncrement;
    double k     = ( 4. / 3. ) * std::tan( theta / 4 );
 
    /*
     * Define our Bezier:
     * - Start point is on the circle at the x-axis
     * - First control point just uses k as the y-value
     * - Second control point is offset from the end point
     * - End point is defined by the angle of the arc segment
     * Note that we use double here no matter what the template param is; round at the end only.
     */
    BEZIER<double> first = { { 1, 0 },
                             { 1, k },
                             { std::cos( theta ) + k * std::sin( theta ),
                               std::sin( theta ) - k * std::cos( theta ) },
                             { std::cos( theta ), std::sin( theta ) } };
 
    /*
     * Now construct the actual segments by transforming/rotating the first one
     */
    auto transformPoint =
            [&]( VECTOR2D aPoint, const double aAngle ) -> VECTOR2D
            {
                // Bring to the actual starting angle
                RotatePoint( aPoint,
                             -EDA_ANGLE( aAngle - aEllipse.StartAngle.AsRadians(), RADIANS_T ) );
 
                // Then scale to the major and minor radiuses of the ellipse
                aPoint *= VECTOR2D( aEllipse.MajorRadius, aEllipse.MinorRadius );
 
                // Now rotate to the ellipse coordinate system
                RotatePoint( aPoint, -aEllipse.Rotation );
 
                // And finally offset to the center location of the ellipse
                aPoint += aEllipse.Center;
 
                return aPoint;
            };
 
    for( int i = 0; i < numBeziers; i++ )
    {
        aBeziers.emplace_back( BEZIER<T>( {
                transformPoint( first.Start, i * angleIncrement ),
                transformPoint( first.C1,    i * angleIncrement ),
                transformPoint( first.C2,    i * angleIncrement ),
                transformPoint( first.End,   i * angleIncrement )
        } ) );
    }
}
 
 
template void TransformEllipseToBeziers( const ELLIPSE<double>& aEllipse,
                                         std::vector<BEZIER<double>>& aBeziers );
template void TransformEllipseToBeziers( const ELLIPSE<int>& aEllipse,
                                         std::vector<BEZIER<int>>& aBeziers );