Approximating Data with a Multi Segment Cubic Bezier Curve and a Distance as Well as a Curvature Contraint

Approximating data with a multi segment cubic bezier curve and a distance as well as a curvature contraint

I found the solution that fulfills my criterea. The solution is to first find a B-Spline that approximates the points in the least square sense and then convert that spline into a multi segment bezier curve. B-Splines do have the advantage that in contrast to bezier curves they will not pass through the control points as well as providing a way to specify a desired "smoothness" of the approximation curve. The needed functionality to generate such a spline is implemented in the FITPACK library to which scipy offers a python binding. Lets suppose I read my data into the lists x and y, then I can do:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
tck,u = interpolate.splprep([x,y],s=3)
unew = np.arange(0,1.01,0.01)
out = interpolate.splev(unew,tck)
plt.figure()
plt.plot(x,y,out[0],out[1])
plt.show()

The result then looks like this:

If I want the curve more smooth, then I can increase the s parameter to splprep. If I want the approximation closer to the data I can decrease the s parameter for less smoothness. By going through multiple s parameters programatically I can find a good parameter that fits the given requirements.

The question though is how to convert that result into a bezier curve. The answer in this email by Zachary Pincus. I will replicate his solution here to give a complete answer to my question:

def b_spline_to_bezier_series(tck, per = False):
  """Convert a parametric b-spline into a sequence of Bezier curves of the same degree.

  Inputs:
    tck : (t,c,k) tuple of b-spline knots, coefficients, and degree returned by splprep.
    per : if tck was created as a periodic spline, per *must* be true, else per *must* be false.

  Output:
    A list of Bezier curves of degree k that is equivalent to the input spline. 
    Each Bezier curve is an array of shape (k+1,d) where d is the dimension of the
    space; thus the curve includes the starting point, the k-1 internal control 
    points, and the endpoint, where each point is of d dimensions.
  """
  from fitpack import insert
  from numpy import asarray, unique, split, sum
  t,c,k = tck
  t = asarray(t)
  try:
    c[0][0]
  except:
    # I can't figure out a simple way to convert nonparametric splines to 
    # parametric splines. Oh well.
    raise TypeError("Only parametric b-splines are supported.")
  new_tck = tck
  if per:
    # ignore the leading and trailing k knots that exist to enforce periodicity 
    knots_to_consider = unique(t[k:-k])
  else:
    # the first and last k+1 knots are identical in the non-periodic case, so
    # no need to consider them when increasing the knot multiplicities below
    knots_to_consider = unique(t[k+1:-k-1])
  # For each unique knot, bring it's multiplicity up to the next multiple of k+1
  # This removes all continuity constraints between each of the original knots, 
  # creating a set of independent Bezier curves.
  desired_multiplicity = k+1
  for x in knots_to_consider:
    current_multiplicity = sum(t == x)
    remainder = current_multiplicity%desired_multiplicity
    if remainder != 0:
      # add enough knots to bring the current multiplicity up to the desired multiplicity
      number_to_insert = desired_multiplicity - remainder
      new_tck = insert(x, new_tck, number_to_insert, per)
  tt,cc,kk = new_tck
  # strip off the last k+1 knots, as they are redundant after knot insertion
  bezier_points = numpy.transpose(cc)[:-desired_multiplicity]
  if per:
    # again, ignore the leading and trailing k knots
    bezier_points = bezier_points[k:-k]
  # group the points into the desired bezier curves
  return split(bezier_points, len(bezier_points) / desired_multiplicity, axis = 0)

So B-Splines, FITPACK, numpy and scipy saved my day :)

Cubic bezier curves - get Y for given X - special case where X of control points is increasing

First off: this answer only works because your control point constraint means we're always dealing with a parametric equivalent of a normal function. Which is obviously what you want in this case, but anyone who finds this answer in the future should be aware that this answer is based on the assumption that there is only one y value for any given x value...

This is absolutely not true for Bezier curves in general

With that said, we know that—even though we've expressed this curve as a parametric curve in two dimensions—we're dealing with a curve that for all intents and purposes must have some unknown function of the form y = f(x). We also know that as long as we know the "t" value that uniquely belongs to a specific x (which is only the case because of your strictly monotonically increasing coefficients property), we can compute y as y = By(t), so the question is: can we compute the t value that we need to plug into By(t), given some known x value?

To which the answer is: yes, we can.

First, any x value we start off with can be said to come from x = Bx(t), so given that we know x, we should be able to find the corresponding value(s) of t that leads to that x.

let's look at the function for x(t):

x(t) = a(1-t)³ + 3b(1-t)²t + 3c(1-t)t² + dt³

We can rewrite this to a plain polynomial form as:

x(t) = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + a

This is a standard cubic polynomial, with only known constants as coefficients, and we can trivially rewrite this to:

(-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + (a-x) = 0

You might be wondering "where did all the -x go for the other values a, b, c, and d?" and the answer there is that they all cancel out, so the only one we actually end up needing to subtract is the one all the way at the end. Handy!

So now we just... solve this equation: we know everything except t, we just need some mathematical insight to tell us how to do this.

...Of course "just" is not the right qualifier here, there is nothing "just" about finding the roots of a cubic function, but thankfully, Gerolano Cardano laid the ground works to determining the roots back in the 16^th century, using complex numbers. Before anyone had even invented complex numbers. Quite a feat! But this is a programming answer, not a history lesson, so let's get implementing:

Given some known value for x, and a knowledge of our coordinates a, b, c, and d, we can implement our root-finding as follows:

// Find the roots for a cubic polynomial with bernstein coefficients
// {pa, pb, pc, pd}. The function will first convert those to the
// standard polynomial coefficients, and then run through Cardano's
// formula for finding the roots of a depressed cubic curve.
double[] findRoots(double x, double pa, double pb, double pc, double pd) {
  double
    pa3 = 3 * pa,
    pb3 = 3 * pb,
    pc3 = 3 * pc,
    a = -pa  +   pb3 - pc3 + pd,     
    b =  pa3 - 2*pb3 + pc3, 
    c = -pa3 +   pb3, 
    d =  pa  -     x;

  // Fun fact: any Bezier curve may (accidentally or on purpose)
  // perfectly model any lower order curve, so we want to test 
  // for that: lower order curves are much easier to root-find.
  if (approximately(a, 0)) {
    // this is not a cubic curve.
    if (approximately(b, 0)) {
      // in fact, this is not a quadratic curve either.
      if (approximately(c, 0)) {
        // in fact in fact, there are no solutions.
        return new double[]{};
      }
      // linear solution:
      return new double[]{-d / c};
    }
    // quadratic solution:
    double
      q = sqrt(c * c - 4 * b * d), 
      b2 = 2 * b;
    return new double[]{
      (q - c) / b2, 
      (-c - q) / b2
    };
  }

  // At this point, we know we need a cubic solution,
  // and the above a/b/c/d values were technically
  // a pre-optimized set because a might be zero and
  // that would cause the following divisions to error.

  b /= a;
  c /= a;
  d /= a;

  double
    b3 = b / 3,
    p = (3 * c - b*b) / 3, 
    p3 = p / 3, 
    q = (2 * b*b*b - 9 * b * c + 27 * d) / 27, 
    q2 = q / 2, 
    discriminant = q2*q2 + p3*p3*p3, 
    u1, v1;

  // case 1: three real roots, but finding them involves complex
  // maths. Since we don't have a complex data type, we use trig
  // instead, because complex numbers have nice geometric properties.
  if (discriminant < 0) {
    double
      mp3 = -p/3,
      r = sqrt(mp3*mp3*mp3), 
      t = -q / (2 * r), 
      cosphi = t < -1 ? -1 : t > 1 ? 1 : t, 
      phi = acos(cosphi), 
      crtr = crt(r), 
      t1 = 2 * crtr;
    return new double[]{
      t1 * cos(phi / 3) - b3,
      t1 * cos((phi + TAU) / 3) - b3,
      t1 * cos((phi + 2 * TAU) / 3) - b3
    };
  }

  // case 2: three real roots, but two form a "double root",
  // and so will have the same resultant value. We only need
  // to return two values in this case.
  else if (discriminant == 0) {
    u1 = q2 < 0 ? crt(-q2) : -crt(q2);
    return new double[]{
      2 * u1 - b3,
      -u1 - b3
    };
  }

  // case 3: one real root, 2 complex roots. We don't care about
  // complex results so we just ignore those and directly compute
  // that single real root.
  else {
    double sd = sqrt(discriminant);
    u1 = crt(-q2 + sd);
    v1 = crt(q2 + sd);
    return new double[]{u1 - v1 - b3};
  }
}

Okay, that's quite the slab of code, with quite a few additionals:

• crt() is the cuberoot function. We actually don't care about complex numbers in this case so the easier way to implement this is with a def, or macro, or ternary, or whatever shorthand your language of choice offers: crt(x) = x < 0 ? -pow(-x, 1f/3f) : pow(x, 1f/3f);.
• tau is just 2π. It's useful to have around when you're doing geometry programming.
• approximately is a function that compares a value to a very small interval around the target because IEEE floating point numerals are jerks. Basically we're talking about approximately(a,b) = return abs(a-b) < 0.000001 or something.

The rest should be fairly self-explanatory, if a little java-esque (I'm using Processing for these kind of things).

With this implementation, we can write our implementation to find y, given x. Which is a little more involved than just calling the function because cubic roots are complicated things. You can get up to three roots back. But only one of them will be applicable to our "t interval" [0,1]:

double x = some value we know!
double[] roots = getTforX(x);
double t;
if (roots.length > 0) {
  for (double _t: roots) {
    if (_t<0 || _t>1) continue;
    t = _t;
    break;
  }
}

And that's it, we're done: we now have the "t" value that we can use to get the associated "y" value.

How can I fit a Bézier curve to a set of data?

I have similar problem and I have found "An algorithm for automatically fitting digitized curves" from Graphics Gems (1990) about Bezier curve fitting.
Additionally to that I have found source code for that article.

Unfortunately it is written in C which I don't know very well. Also, the algorithm is quite hard to understand (at least for me). I am trying to translate it into C# code. If I will be successful, I will try to share it.

File GGVecLib.c in the same folder as FitCurves.c contains basic vectors manipulation functions.

I have found a similar Stack Overflow question, Smoothing a hand-drawn curve. The approved answer provide C# code for a curve fitting algorithm from Graphic Gems.

Computing the Radius of Curvature of a Bezier Curve Given Control Points

The radius of curvature R(t) is equal to 1/κ(t), where κ(t) is the curvature of the curve at point t, which for a parametric planar curve is:

            x'y" - y'x"
κ(t)  = --------------------
         (x'² + y'²)^(3/2)

(Where the ^(3/2) is really ^3/2, but you can't use html for superscript formatting in code blocks, and where the (t) part has been left off the functions for x and y because that makes things unnecessarily hard to read)

As such, for a Quadratic Bezier curve with controls P₁, P₂, and P₃, the first and second derivatives use the following control points:

B(t)': P₁' = 2(P₂ - P₁), and P₂' = 2(P₃ - P₂)
B(t)": P₁" = (P'₂ - P'₁)

Evaluating these for x and y are literally just "using the x or y coordinates", so:

x' = Px₁'(t-1) + Px₂'(t)
y' = Py₁'(t-1) + Py₂'(t)
x" = Px₁"
y" = Py₁"

Noting that x" and y" are just constants. We plug these values into the function for κ(t), provided that the denominator isn't zero, of course (which indicates a line segment, which has no radius of curvature) and then we know what R(t) is because that's just the inverse value.

But we don't actually need to do any of that, because we can make a simple observation: treating the circle as a parametric curve as well, we can find its radius pretty much immediately based on the fact that the vector length of the derivative of the curve, and the derivative of the circle at the same point, are equal. Additionally, we know that the length of the derivative is the same anywhere on a circle (because it's a circle, which by definition is radially symmetric about the origin). As such, we can solve the following:

B(t)  = P1(1-t)^2 + 2*P2(1-t)t + P3t^2
B(t)' = 2*(P2-P1)(1-t) + 2*(P3-P2)t

C(s)  = { r*sin(s), r*cos(s) }independent of 't')
C(s)' = { r*cos(s), -r*sin(s) }

d = |B(t)'| = |C(s = any value, so let's pick 0)'|
d = C(0)' = | (r,0) | = r

And we're done. The radius of curvature at point B(t) is equal to the vector length of B(t)', which is much easier to implement, and much faster to run.

Finally, if what you actually want is to just approximate sections of curve with circular arcs, then https://pomax.github.io/bezierinfo/#arcapproximation should cover the "how to" for that particular use-case.

How to change a Bezier curve by moving a point on it with the mouse?

You misunderstood the answer from tfinniga's post.

From tfinniga's post, we have

P = k0*P0 + k1*P1 + k2*P2 + k3*P3 and 
P' = k0*P0' + k1*P1' + k2*P2' + k3*P3'

Since you require P0 and P3 unchanged, we have two identities for V

V = k1*(P1'-P1) + k2*(P2'-P2)

and

V = P' - P0

So, you can choose

  P1' = P1 + c/k1 * V, 
  P2' = P2 + (1-c)/k2 * V

where c is a constant between 0 and 1.

How to clear my BezierCurve vector to stop writing coordinates on top of each other?

You will want to detect when mouse click changes from down to up:

bool prevMouse;
...
// at the end of Display()
prevMouse = MouseReleased;

Then we check when mouse click changes from pressed to released and add lines to BezierCurve:

if (PrevMouse == 0 && MouseReleased)
{
    // Store curvature into Bezier Points
    BezierCurve.push_back(p2);
}

The two code paths, if (Points.size() == 2), and else if (Points.size() > 2) inside the for loop can be simplified to if (Points.size() >= 2) and for that matter the for loop is extraneous, we don't need to update the bezier curve for any of the previous points, just the curve between the two newest points, Points[Points.size() - 2] and Points[Points.size() - 1].

The final code:

#include <iostream>
#include <stdlib.h>
#include <GL/glut.h>
#include <vector>
#include <math.h>

using namespace std;

//Point class for taking the points
class Point {
public:
    float x, y;

    void setxy(float x2, float y2)
    {
        x = x2; y = y2;
    }

    //operator overloading for '=' sign
    const Point& operator=(const Point& rPoint)
    {
        x = rPoint.x;
        y = rPoint.y;
        return *this;
    }

};

int SCREEN_HEIGHT = 500;
vector<Point> Points;
Point Tangent;
Point inverseTangent;
Point cursorLocationLive;
int TangentsSize = 0;
vector<Point> Tangents(TangentsSize);
vector<Point> inverseTangents(TangentsSize);
vector<Point> BezierCurve;
bool MouseReleased = false;
bool PrevMouse = false;

void drawDot(Point p1)
{
    glBegin(GL_POINTS);
    glVertex2i(p1.x, p1.y);
    glEnd();
}

void drawLine(Point p1, Point p2)
{
    glBegin(GL_LINE_STRIP);
    glVertex2f(p1.x, p1.y);
    glVertex2f(p2.x, p2.y);
    glEnd();
}

float interpolate(float n1, float n2, float perc)
{
    float diff = n2 - n1;

    return n1 + (diff * perc);
}

void myMouse(int button, int state, int x, int y)
{
    if (button == GLUT_LEFT_BUTTON)
    {
        if (state == GLUT_DOWN)
        {
            MouseReleased = false;
            // Store points into Points vector on click
            Point point;
            point.setxy(x, SCREEN_HEIGHT - y);
            Points.push_back(point);

            // Tangents are set to the cursor position
            Tangent.setxy(x, SCREEN_HEIGHT - y);
            inverseTangent.x = (2 * Points[Points.size() - 1].x) - Tangent.x;
            inverseTangent.y = (2 * Points[Points.size() - 1].y) - Tangent.y;

            /*Add new element to Tangent & inverseTangent so when we draw the curve
            the tangents are accessed at the right index*/
            TangentsSize++;
        }
        else if (state == GLUT_UP)
        {
            MouseReleased = true;
            // Upon mouse release store tangent and inverse tangent into separate vectors
            Tangents.push_back(Tangent);
            inverseTangents.push_back(inverseTangent);
        }
    }
}

void passiveMotion(int x, int y)
{
    // Sets the location of cursor while moving with no buttons pressed
    cursorLocationLive.setxy(x, SCREEN_HEIGHT - y);
}

void motion(int x, int y)
{
    // Sets the coordinates of the tangents when mouse moves with a button held down
    Tangent.setxy(x, SCREEN_HEIGHT - y);
    inverseTangent.x = (2 * Points[Points.size() - 1].x) - Tangent.x;
    inverseTangent.y = (2 * Points[Points.size() - 1].y) - Tangent.y;
}

void myDisplay()
{
    glClear(GL_COLOR_BUFFER_BIT);

    // Draw main points in red
    glColor3f(255, 0, 0);
    for (int i = 0; i < Points.size(); i++)
    {
        drawDot(Points[i]);
    }

    // If there is a starting point draw a line to cursor from last drawn point in passive motion
    if (Points.size() > 0)
    {
        glColor3f(0, 0, 0);
        drawLine(Points[Points.size() - 1], cursorLocationLive);
    }

    // Draw live tangent dots in green
    glColor3f(0, 255, 0);
    drawDot(Tangent);
    drawDot(inverseTangent);
    // Draw live tangent lines in blue
    glColor3f(0, 0, 255);
    drawLine(Tangent, inverseTangent);

    for (int i = 0; i < Tangents.size(); i++)
    {
        // Draw stored tangent dots in green
        glColor3f(0, 255, 0);
        drawDot(Tangents[i]);
        drawDot(inverseTangents[i]);
        // Draw stored tangent lines in blue
        glColor3f(0, 0, 255);
        drawLine(Tangents[i], inverseTangents[i]);
    }

    // Loop through all points
    if (Points.size() >= 2)
    {
        // p1 is the start of the curve set to second last point
        Point p1;
        p1 = Points[Points.size() - 2];

        float i;
        // Calculate curve coordinates
        for (float j = 0; j <= 100; j++)
        {
            i = j / 100;
            // The Green Lines
            float xa = interpolate(Points[Points.size() - 2].x, Tangents[TangentsSize - 2].x, i);
            float ya = interpolate(Points[Points.size() - 2].y, Tangents[TangentsSize - 2].y, i);
            float xb = interpolate(Tangents[TangentsSize - 2].x, inverseTangent.x, i);
            float yb = interpolate(Tangents[TangentsSize - 2].y, inverseTangent.y, i);
            float xc = interpolate(inverseTangent.x, Points[Points.size() - 1].x, i);
            float yc = interpolate(inverseTangent.y, Points[Points.size() - 1].y, i);

            // The Blue Line
            float xm = interpolate(xa, xb, i);
            float ym = interpolate(ya, yb, i);
            float xn = interpolate(xb, xc, i);
            float yn = interpolate(yb, yc, i);

            // The Black Dot
            float x2 = interpolate(xm, xn, i);
            float y2 = interpolate(ym, yn, i);

            Point p2;
            p2.setxy(x2, y2);

            drawLine(p1, p2);

            p1 = p2;

            // Prevents curves generated during mouse motion from being stored
            if (PrevMouse == 0 && MouseReleased)
            {
                // Store curvature into Bezier Points
                BezierCurve.push_back(p2);
            }
        }
    }
    std::cout << BezierCurve.size() << std::endl;
    
    PrevMouse = MouseReleased;

    // Draw all bezier curvature
    for (int i = 1; i < BezierCurve.size(); i++)
    {
        drawLine(BezierCurve[i - 1], BezierCurve[i]);
    }

    glutSwapBuffers();
}

void timer(int)
{
    glutTimerFunc(1000 / 60, timer, 0);
    glutPostRedisplay();
}

int main(int argc, char* argv[]) {
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
    glutInitWindowSize(640, 500);
    glutInitWindowPosition(100, 150);
    glutCreateWindow("Bezier Curve");
    glutDisplayFunc(myDisplay);
    glutIdleFunc(myDisplay);
    glutTimerFunc(0, timer, 0);
    glutMouseFunc(myMouse);
    glutPassiveMotionFunc(passiveMotion);
    glutMotionFunc(motion);
    glClearColor(255, 255, 255, 0.0);
    glPointSize(3);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0.0, 640.0, 0.0, 500.0);
    glutMainLoop();

    return 0;
}

---

Related Topics