Bézier Curves

Bézier Curves

Bézier curves are a type of parametric curve that is widely used in computer graphics, computer-aided design, and various other fields for creating smooth, elegant curves. They are named after the French engineer Pierre Bézier, who developed the mathematics behind these curves while working at Renault in the 1960s.

Bézier curves come in different degrees, with the most common being the quadratic (degree 2), cubic (degree 3), and higher-degree Bézier curves. Higher-degree Bézier curves can represent more complex shapes.

Linear Segment

Given two points of initial segment, we can calculate the new point on the line AB by the Linear equation:

Quadratic Bézier

Divide segments in half and connect midpoints to get the new position for the original vertex. In this process, the new points and the points moved will keep their position for remaining iterations. The endpoints never move.

Let u be the distance between the begin point p0 and the mid point p3 of the first segment.

Write a simple algorithm implementation for it, the left renderer shows original four points, and the right renderer shows the result of smoothing.

point.hpp

#ifndef POINT_HPP
#define POINT_HPP

#include <iostream>
#include <math.h>
using namespace std;

class Point
{
public:
    double point[3];
    Point() {}
    Point(double *ptr) { point[0]=ptr[0]; point[1]=ptr[1]; point[2]=ptr[2];}
    Point(double x, double y, double z) { point[0]=x; point[1]=y; point[2]=z;}
    Point & SetPoint(double p[]) { point[0]=p[0]; point[1]=p[1];  point[2]=p[2];  return *this; }
    double operator[](int i) const { return point[i]; }

    Point operator+(double v) { return Point(point[0]+v, point[1]+v, point[2]+v); }
    Point operator-(double v) { return Point(point[0]-v, point[1]-v, point[2]-v); }
    Point operator*(double v) { return Point(point[0]*v, point[1]*v, point[2]*v); }
    Point operator/(double v) { return Point(point[0]/v, point[1]/v, point[2]/v); }
    Point operator-() { return Point(-point[0], -point[1], -point[2]); }
    Point operator/=(double v) { point[0]/=v; point[1]/=v; point[2]/=v;   return *this; }
    double Dot(const Point &p) { return point[0]*p.point[0]+point[1]*p.point[1]+point[2]*p.point[2]; }
    double Length() { return sqrt(point[0]*point[0]+point[1]*point[1]+point[2]*point[2]); }
    Point Unit() { return (operator/=(Length())); }

    friend Point operator*(double v,Point &p) { return p*v; };
    friend Point operator/( Point &p, double v) { return p*(1./v); }
    friend Point operator+(const Point &p1, const Point &p2) { return Point(p1[0]+p2[0],p1[1]+p2[1],p1[2]+p2[2]); }
    friend Point operator-(const Point &p1, const Point &p2) { return Point(p1[0]-p2[0],p1[1]-p2[1],p1[2]-p2[2]); }
    friend Point operator^(const Point &p1, const Point &p2) { return Point(p1[1]*p2[2]-p1[2]*p2[1],p1[2]*p2[0]-p1[0]*p2[2],p1[0]*p2[1]-p1[1]*p2[0]); }
    friend bool operator==( Point &p1,  Point &p2)
    {
        bool notEqual = (p1 != p2);
        return !notEqual;
    }
    friend bool operator!=( Point &p1,  Point &p2)
    {
        if( fabs( p1[0] - p2[0] ) > 1e-6 ||
            fabs( p1[1] - p2[1] ) > 1e-6 ||
            fabs( p1[2] - p2[2] ) > 1e-6 )
        {
            return true;
        }
        return false;
    }
    bool operator ==( const Point &p2 ) const
    {
        if( fabs( point[0] - p2[0] ) > 1e-6 ||
            fabs( point[1] - p2[1] ) > 1e-6 ||
            fabs( point[2] - p2[2] ) > 1e-6 )
        {
            return false;
        }
        return true;
    }
    bool operator<( const Point &p2 ) const
    {
        return ( point[0] <=  p2[0] && point[1] <= p2[1] && point[2] <= p2[2] );
    }
};

#endif //POINT_HPP

main.cpp

#include <vtkActor.h>
#include <vtkCleanPolyData.h>
#include <vtkDistancePolyDataFilter.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkPolyDataMapper.h>
#include <vtkPolyDataReader.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkScalarBarActor.h>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkSTLReader.h>

#include <vector>
#include <list>
#include "point.hpp"

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

std::vector<Point> originalPts = { Point(1, 0, 0), Point(3, 3, 0), Point(7, 5, 0), Point(9, 0, 0) };

vtkSmartPointer<vtkPolyData> ShowListByLine(std::vector<Point> list)
{
    vtkSPtrNew( result, vtkPolyData );
    vtkSPtrNew( resultPts, vtkPoints );
    vtkSPtrNew( resultLines, vtkCellArray );
    for( int i = 0; i < list.size(); ++i )
    {
        resultPts->InsertNextPoint( list[i].point );
    }
    for( int i = 0; i < list.size()-1; ++i )
    {
        vtkIdType pts[2] = { i, i + 1 };
        resultLines->InsertNextCell( 2, pts );
    }
    result->SetPoints( resultPts );
    result->SetLines( resultLines );
    result->Modified();
    return result;
}

std::vector<Point> HandleNumSet(std::vector<Point> numSet)
{
    Point p3 = (numSet[0] + numSet[1])*0.5;
    Point p4 = (numSet[1] + numSet[2])*0.5;
    Point p1_ = (p3 + p4)*0.5;
    return std::vector<Point>{ p3, p1_, p4 };
}

std::list<Point> SubdivisionProcedure( std::vector<Point> inputPoints, std::vector<Point> &markedPts )
{
    std::list<Point> curList;
    for( auto value : inputPoints ){ curList.push_back( value ); }

    for( int i = 1; i < inputPoints.size() - 1; ++i )
    {
        // the points had been changed will be ignore.
        if( std::find( markedPts.begin(), markedPts.end(), inputPoints[i] ) != markedPts.end() ) continue;

        std::list<Point>::iterator lastOne = std::find(curList.begin(), curList.end(), inputPoints[i] );
        std::advance( lastOne, -1);
        std::list<Point>::iterator nextOne = std::find(curList.begin(), curList.end(), inputPoints[i] );
        std::advance( nextOne, 1);

        std::vector<Point> numSet = { *lastOne, inputPoints[i], *nextOne };
        std::vector<Point> newPts = HandleNumSet( numSet );

        std::list<Point>::iterator pos = std::find(curList.begin(), curList.end(), inputPoints[i] );
        curList.insert( pos, newPts[0] );
        curList.insert( pos, newPts[1] );
        curList.insert( pos, newPts[2] );
        curList.erase( pos );

        markedPts.push_back( newPts[1] );
    }
    return curList;
}

int main(int argc, char* argv[])
{
    vtkSmartPointer<vtkPolyData> input1 = ShowListByLine( originalPts );

    std::vector<Point> markedPts;
    std::vector<Point> samplePts = originalPts;
    std::list<Point> pts;
    int iterationTimes = 13;
    while( iterationTimes-- )
    {
        pts = SubdivisionProcedure( samplePts, markedPts );
        samplePts.clear();
        for( auto element: pts  )
        {
            samplePts.push_back( element );
        }
    }

    std::vector<Point> inputPts;
    for (const auto& element : pts) {
        inputPts.push_back( element );
    }

    vtkSmartPointer<vtkPolyData> input2 = ShowListByLine( inputPts );

    vtkNew<vtkRenderWindow> renWin;
    renWin->SetSize(1200, 500);

    vtkNew<vtkRenderWindowInteractor> renWinInteractor;
    renWinInteractor->SetRenderWindow(renWin);

    double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
    double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};

    vtkSPtrNew( mapper1, vtkPolyDataMapper );
    mapper1->SetInputData( input1 );
    vtkSPtrNew( actor1, vtkActor );
    actor1->SetMapper( mapper1 );

    vtkSPtrNew( mapper2, vtkPolyDataMapper );
    mapper2->SetInputData( input2 );
    vtkSPtrNew( actor2, vtkActor );
    actor2->SetMapper( mapper2 );

    // Setup renderers
    vtkSPtrNew( leftRenderer, vtkRenderer );
    leftRenderer->SetViewport( leftViewport );
    leftRenderer->AddActor( actor1 );
    leftRenderer->SetBackground(.6, .5, .4);
    leftRenderer->ResetCamera();

    vtkSPtrNew( rightRenderer, vtkRenderer );
    rightRenderer->SetViewport(rightViewport);
    rightRenderer->AddActor( actor2 );
    rightRenderer->SetBackground(.4, .5, .6);
    rightRenderer->SetActiveCamera( leftRenderer->GetActiveCamera() );

    renWin->AddRenderer( leftRenderer );
    renWin->AddRenderer( rightRenderer );

    renWin->Render();
    renWinInteractor->Start();


    return EXIT_SUCCESS;
}

Cubic Bézier

Similarlly, we can take first four points to get midpoints and change original vertices’ positions.
Here is a process to smooth polygon, the image shows the de Casteljau algorithm to get cubic curve.

if we always use mid point of segments, we will get the position of point AD is:

The general form of it about distance from A to AB (u) is :