Converting straight lines segments into a bezier curve of same shape

Hello, everyone!
I am making a typeface, styles of which are generated by algorithm from one root style.
the main alogritm is raster-based (reaction-diffusion). Then I make a super detailed PVector curve by finding the edge pixel by pixel. At final I simplify this curve with an RDP algorithm. All this gives me a quite good result, but the curve is still too rough because it consists of straight segments.
Is there any way to convert a bunch of straight segments into nice curves? I know the way must exist.

On the picture you can see what I get now and what I want to get

You can use illustrator or Inkscape for this. They both have a raster to vector functionality.

Thanks but it’s not an option. The point is to process hundreds of glyphs automatically, in one script. Also, illustrator’s tracing functionality works quite bad to me

I found An Algorithm for Automatically Fitting Digitized Curves by Philip J. Schneider, written on C, Python and JS. Hopefully, I will manage to port it into Processing

With illustrator you can process files in batch. So if you have all your input files in the same folder, you can ask illustrator to vectorize and save each one in pretty much one click.

Also, the vectorizing tool of illustrator should be the state of the art vectorizing tool so I would be surprise if you could get better results.

Have you also tried to look at some already existing library? I found this one for example: GitHub - jankovicsandras/imagetracerjava: Simple raster image tracer and vectorizer written in Java for desktop

1 Like

This should be of use… I’ve had this sitting for 13 months in an unfinished state but your post has encouraged me to get it over the line. After fitting, the output is a sequence of cubic curves.


@micycle Thanks for sharing! I am terribly sorry for such a long silence. This task fell completely out of my focus for a long time for a serious reason. But I am getting back to it.

I’ve tried your algorithms. They work, it’s amazing! But how can I adjust the accuracy level? The result I get is quite rough. What I am trying to get must be much more precise.

I somehow ported the algorithm by Philip J. Schneider from JS and Python to Processing. It works but not perfectly. Sometimes it makes too big inaccuracies and next to them — unnecessary short straight segments. I guess it is my fault, but I can not find the cause :frowning: I would be happy if someone can help me to find mistakes in the code. I want to publish it as a library in the end.

Here is the ported code. It uses PVector_D — it is a clone of PVector class but with doubles instead of floats. With floats it was glitching sometimes.

ArrayList <PVector_D[]> fitCurve(ArrayList <PVector_D> points, double maxError) {
  //Remove duplicate points
  for (int i = points.size()-1; i>0; i--) {
    PVector_D p = points.get(i);
    PVector_D prevP = points.get(i-1);
    if (p.x == prevP.x && p.y == prevP.y) points.remove(i);

  if (points.size() < 2) {
    return null;

  int len = points.size();
  PVector_D leftTangent = createTangent(points.get(1), points.get(0));
  PVector_D rightTangent = createTangent(points.get(len - 2), points.get(len - 1));
  return fitCubic(points, leftTangent, rightTangent, maxError);

 * Fit a Bezier curve to a (sub)set of digitized points.
 * Your code should not call this function directly. Use {@link fitCurve} instead.
 * @param {Array<Array<Number>>} points - Array of digitized points, e.g. [[5,5],[5,50],[110,140],[210,160],[320,110]]
 * @param {Array<Number>} leftTangent - Unit tangent vector at start point
 * @param {Array<Number>} rightTangent - Unit tangent vector at end point
 * @param {Number} error - Tolerance, squared error between points and fitted curve
 * @returns {Array<Array<Array<Number>>>} Array of Bezier curves, where each element is [first-point, control-point-1, control-point-2, second-point] and points are [x, y]

ArrayList <PVector_D[]> fitCubic(ArrayList <PVector_D> points, PVector_D leftTangent, PVector_D rightTangent, double error) {
  int MaxIterations = 100; //Max times to try iterating (to find an acceptable curve)
  PVector_D centerVector, toCenterTangent, fromCenterTangent; //Unit tangent vector(s) at splitPoint
  PVector_D [] bezCurve = new PVector_D [4]; //Control points of fitted Bezier curve
  double [] u = new double [points.size()]; //Parameter values for point
  double [] uPrime = new double [points.size()]; //Improved parameter values
  double maxError ; //Maximum fitting error
  double prevErr; 
  int splitPointIndex; //Point to split point set at if we need more than one curve
  int prevSplitIndex; 
  ArrayList <PVector_D[]> beziers = new ArrayList <PVector_D []> (); //Array of fitted Bezier curves if we need more than one curve
  double dist;

  //Use heuristic if region only has two points in it
  if (points.size() == 2) {
    PVector_D tempVector = PVector_D.sub(points.get(0), points.get(1));
    dist = tempVector.mag() / 3.0;
    bezCurve [0] = points.get(0);
    bezCurve [1] = PVector_D.add(points.get(0), PVector_D.mult(leftTangent, dist));
    bezCurve [2] = PVector_D.add(points.get(1), PVector_D.mult(rightTangent, dist));
    bezCurve [3] = points.get(1);
    return beziers;

  //Parameterize points, and attempt to fit curve
  u = chordLengthParameterize(points);
  bezCurve = generateBezier(points, u, leftTangent, rightTangent);
  //println("first attempt:");
  //println(bezCurve[0], bezCurve[1], bezCurve[2], bezCurve[3]);

  ArrayList temp = computeMaxError(points, bezCurve, u);
  maxError =  (double)temp.get(0);
  splitPointIndex = (int)temp.get(1);

  if (maxError == 0 || maxError < error) {
    return beziers;

  //If error not too large, try some reparameterization and iteration
  if (maxError < error * error) {

    uPrime = u;
    prevErr = maxError;
    prevSplitIndex = splitPointIndex;

    for (int i = 0; i < MaxIterations; i++) {

      uPrime = reparameterize(bezCurve, points, uPrime);
      bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent);
      //println("second attempt:");
      //println(bezCurve[0], bezCurve[1], bezCurve[2], bezCurve[3]);
      ArrayList temp2 = computeMaxError(points, bezCurve, u);
      maxError =  (double)temp2.get(0);
      splitPointIndex = (int)temp2.get(1);

      if (maxError < error) {
        return beziers;
      } else if (splitPointIndex == prevSplitIndex) {
        double errChange = maxError / prevErr;
        if (errChange > 0.9999 && errChange < 1.0001) {

      prevErr = maxError;
      prevSplitIndex = splitPointIndex;
  //Fitting failed -- split at max error point and fit recursively
  beziers = new ArrayList <PVector_D []> ();

  //To create a smooth transition from one curve segment to the next, we
   //calculate the line between the points directly before and after the
   //center, and use that as the tangent both to and from the center point.
   centerVector = PVector_D.sub(points.get(splitPointIndex - 1), points.get(splitPointIndex + 1));
   //However, this won't work if they're the same point, because the line we
   //want to use as a tangent would be 0. Instead, we calculate the line from
   //that "double-point" to the center point, and use its tangent.
   if (centerVector.x == 0 && centerVector.y == 0) {
   //[x,y] -> [-y,x]:
   centerVector = PVector_D.sub(points.get(splitPointIndex - 1), points.get(splitPointIndex));
   PVector_D temp3 = new PVector_D (-centerVector.y, centerVector.x);
   centerVector.x = temp3.x;
   centerVector.y = temp3.y;
   toCenterTangent = centerVector.get();
   //To and from need to point in opposite directions:
   fromCenterTangent = PVector_D.mult(toCenterTangent, -1);

  toCenterTangent = createTangent(points.get(splitPointIndex - 1), points.get(splitPointIndex));
  fromCenterTangent = createTangent(points.get(splitPointIndex + 1), points.get(splitPointIndex));

        Note: An alternative to this "divide and conquer" recursion could be to always
   let new curve segments start by trying to go all the way to the end,
   instead of only to the end of the current subdivided polyline.
   That might let many segments fit a few points more, reducing the number of total segments.
   However, a few tests have shown that the segment reduction is insignificant
   (240 pts, 100 err: 25 curves vs 27 curves. 140 pts, 100 err: 17 curves on both),
   and the results take twice as many steps and milliseconds to finish,
   without looking any better than what we already have.
  ArrayList <PVector_D> before = new ArrayList <PVector_D> ();
  for (int i = 0; i < splitPointIndex + 1; i++) {
    PVector_D p = points.get(i);
  ArrayList <PVector_D> after = new ArrayList <PVector_D> ();
  for (int i = splitPointIndex; i < points.size(); i++) {
    PVector_D p = points.get(i);
  beziers.addAll(fitCubic(before, leftTangent, toCenterTangent, error));
  beziers.addAll(fitCubic(after, fromCenterTangent, rightTangent, error));
  return beziers;

 * Use least-squares method to find Bezier control points for region.
 * @param {Array<Array<Number>>} points - Array of digitized points
 * @param {Array<Number>} parameters - Parameter values for region
 * @param {Array<Number>} leftTangent - Unit tangent vector at start point
 * @param {Array<Number>} rightTangent - Unit tangent vector at end point
 * @returns {Array<Array<Number>>} Approximated Bezier curve: [first-point, control-point-1, control-point-2, second-point] where points are [x, y]
PVector_D [] generateBezier(ArrayList <PVector_D> points, double [] parameters, PVector_D leftTangent, PVector_D rightTangent) {
  PVector_D [] bezCurve = new PVector_D [4]; 
  PVector_D firstPoint = points.get(0);
  PVector_D lastPoint = points.get(points.size()-1);
  bezCurve [0] = firstPoint;
  bezCurve [3] = lastPoint;

  //Precomputed rhs for eqn
  PVector_D [][] A;
  //Compute the A's
  A = makeZeroVectors(parameters.length);
  for (int i = 0; i < parameters.length; i++) {
    double u = parameters[i];
    double ux = 1 - u;
    A[i][0] = PVector_D.mult(leftTangent, (3 * u * (ux * ux)));
    A[i][1] = PVector_D.mult(rightTangent, (3 * ux * (u * u)));

  //Matrices C & X
  double [][] C = {{0, 0}, {0, 0}}; 
  double [] X = {0, 0}; 

  for (int i = 0; i < points.size(); i++) {
    double u = parameters[i];

    C[0][0] +=[i][0], A[i][0]);
    C[0][1] +=[i][0], A[i][1]);
    C[1][0] +=[i][0], A[i][1]);
    C[1][1] +=[i][1], A[i][1]);

    PVector_D [] tmp1 = {firstPoint, firstPoint, lastPoint, lastPoint};
    PVector_D tmp2 = PVector_D.sub(points.get(i), bezierQ(tmp1, u));
    X[0] +=[i][0], tmp2);
    X[1] +=[i][1], tmp2);

  //Determinants of matrices C and X
  double det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
  double det_C0_X = C[0][0] * X[1] - C[1][0] * X[0];
  double det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];

  //Alpha values, left and right
  double alpha_l, alpha_r; 
  // Result = test ? expression1 : expression2
  alpha_l = det_C0_C1 == 0 ? 0 : det_X_C1 / det_C0_C1;
  alpha_r = det_C0_C1 == 0 ? 0 : det_C0_X / det_C0_C1;

  //If alpha negative, use the Wu/Barsky heuristic (see text).
  //If alpha is 0, you get coincident control points that lead to
  //divide by zero in any subsequent NewtonRaphsonRootFind() call.
  PVector_D tempVector = PVector_D.sub(firstPoint, lastPoint);
  double segLength = tempVector.mag();
  double epsilon = 1.0e-6 * segLength;
  if (alpha_l < epsilon || alpha_r < epsilon) {
    //Fall back on standard (probably inaccurate) formula, and subdivide further if needed.
    bezCurve[1] = PVector_D.add(firstPoint, PVector_D.mult(leftTangent, segLength / 3.0));
    bezCurve[2] = PVector_D.add(lastPoint, PVector_D.mult(rightTangent, segLength / 3.0));
    //println("the first pair of points: ");
    //println(bezCurve[1], bezCurve[2]);
  } else {
    //First and last control points of the Bezier curve are
    //positioned exactly at the first and last data points
    //Control points 1 and 2 are positioned an alpha distance out
    //on the tangent vectors, left and right, respectively
    //println("parameters are: alfa_l = "+alpha_l+", alfa_l = "+alpha_r+", epsilon = "+epsilon);
    bezCurve[1] = PVector_D.add(firstPoint, PVector_D.mult(leftTangent, alpha_l));
    bezCurve[2] = PVector_D.add(lastPoint, PVector_D.mult(rightTangent, alpha_r));
    //println("the second pair of points: ");
    //println(bezCurve[1], bezCurve[2]);

  return bezCurve;

 * Given set of points and their parameterization, try to find a better parameterization.
 * @param {Array<Array<Number>>} bezier - Current fitted curve
 * @param {Array<Array<Number>>} points - Array of digitized points
 * @param {Array<Number>} parameters - Current parameter values
 * @returns {Array<Number>} New parameter values
double [] reparameterize(PVector_D [] bezier, ArrayList <PVector_D> points, double [] parameters) {
        var j, len, point, results, u;
   results = [];
   for (j = 0, len = points.length; j < len; j++) {
   point = points[j], u = parameters[j];
   results.push(newtonRaphsonRootFind(bezier, point, u));
   return results;

  for (int i = 0; i < parameters.length; i++) {
    parameters[i] = newtonRaphsonRootFind(bezier, points.get(i), parameters[i]);
  return parameters;

 * Use Newton-Raphson iteration to find better root.
 * @param {Array<Array<Number>>} bez - Current fitted curve
 * @param {Array<Number>} point - Digitized point
 * @param {Number} u - Parameter value for "P"
 * @returns {Number} New u

double newtonRaphsonRootFind(PVector_D [] bez, PVector_D point, double u) {
            Newton's root finding algorithm calculates f(x)=0 by reiterating
   x_n+1 = x_n - f(x_n)/f'(x_n)
   We are trying to find curve parameter u for some point p that minimizes
   the distance from that point to the curve. Distance point to curve is d=q(u)-p.
   At minimum distance the point is perpendicular to the curve.
   We are solving
   f = q(u)-p * q'(u) = 0
   f' = q'(u) * q'(u) + q(u)-p * q''(u)
   u_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)|
  PVector_D d = PVector_D.sub(bezierQ(bez, u), point); 
  PVector_D qprime = bezierQPrime(bez, u); 
  double numerator =, qprime); // not sure
  double denominator = summaXY(squareVector(qprime)) + 2*, bezierQPrimePrime(bez, u));

  if (denominator == 0) {
    return u;
  } else {
    return u - numerator / denominator;

 * Assign parameter values to digitized points using relative distances between points.
 * @param {Array<Array<Number>>} points - Array of digitized points
 * @returns {Array<Number>} Parameter values
double [] chordLengthParameterize(ArrayList <PVector_D> points) {
  double [] u = new double [points.size()]; 
  u[0] = 1;
  for (int i = 1; i < points.size(); i ++) {
    PVector_D p = points.get(i);
    PVector_D prevP = points.get(i-1);
    PVector_D tempVector = PVector_D.sub(p, prevP);
    u[i] = u[i-1] + tempVector.mag();
    if (u[i] == 0) println("ALARM 1");

  for (int i = 1; i < u.length; i ++) {
    u[i] = u[i] / u[i-1];

  return u;

 * Find the maximum squared distance of digitized points to fitted curve.
 * @param {Array<Array<Number>>} points - Array of digitized points
 * @param {Array<Array<Number>>} bez - Fitted curve
 * @param {Array<Number>} parameters - Parameterization of points
 * @returns {Array<Number>} Maximum error (squared) and point of max error
ArrayList computeMaxError(ArrayList <PVector_D> points, PVector_D [] bez, double [] parameters) {
  double maxDist = 0; //Maximum error
  int splitPointIndex = floor(points.size() / 2); //Point of maximum error. Let it be in the middle for the beginning

  double [] t_distMap = mapTtoRelativeDistances(bez, 10);

  for (int i = 1; i < points.size()-1; i++) {
    PVector_D point = points.get(i);
    //Find 't' for a point on the bez curve that's as close to 'point' as possible:
    double t = find_t(bez, parameters[i], t_distMap, 10);
    //double t = parameters[i];

    //println("the parameters are:");
    //println("bez: " + bez[0] + bez[1] + bez[2] + bez[3] + " & t = " + t);
    PVector_D somePoint = bezierQ(bez, t);
    //println("the points are: " + somePoint + " & " + point);
    PVector_D v = PVector_D.sub(somePoint, point); //Vector from point to curve
    double dist = v.magSq(); //squared dist
    //println("The actual error is: " + dist);
    if (dist > maxDist) {
      maxDist = dist;
      splitPointIndex = i;
  ArrayList outputs = new ArrayList();
  return outputs;

//Sample 't's and map them to relative distances along the curve:
double [] mapTtoRelativeDistances(PVector_D [] bez, int B_parts) {
  PVector_D B_t_curr;
  double [] B_t_dist = new  double [B_parts];
  PVector_D B_t_prev = bez[0];
  double sumLen = 0;

  for (int i = 0; i < B_parts; i++) {
    B_t_curr = bezierQ(bez, i / B_parts);

    PVector_D tempVector = PVector_D.sub(B_t_curr, B_t_prev);
    sumLen += tempVector.mag();
    B_t_dist[i] = sumLen;
    B_t_prev = B_t_curr;

  //Normalize B_length to the same interval as the parameter distances; 0 to 1:
  for (int i = 0; i < B_t_dist.length; i ++) {
    double f = B_t_dist[i];
    f /= sumLen;
    B_t_dist[i] = f;

  return B_t_dist;

double find_t(PVector_D [] bez, double param, double [] t_distMap, int B_parts) {
  if (param < 0) return 0;
  if (param > 1) return 1;

            'param' is a value between 0 and 1 telling us the relative position
   of a point on the source polyline (linearly from the start (0) to the end (1)).
   To see if a given curve - 'bez' - is a close approximation of the polyline,
   we compare such a poly-point to the point on the curve that's the same
   relative distance along the curve's length.
   But finding that curve-point takes a little work:
   There is a function "B(t)" to find points along a curve from the parametric parameter 't'
   (also relative from 0 to 1:,
   but 't' isn't linear by length (
   So, we sample some points along the curve using a handful of values for 't'.
   Then, we calculate the length between those samples via plain euclidean distance;
   B(t) concentrates the points around sharp turns, so this should give us a good-enough outline of the curve.
   Thus, for a given relative distance ('param'), we can now find an upper and lower value
   for the corresponding 't' by searching through those sampled distances.
   Finally, we just use linear interpolation to find a better value for the exact 't'.
   More info:
  double lenMax, lenMin, tMax, tMin;
  double t = 0;

  //Find the two t-s that the current param distance lies between,
  //and then interpolate a somewhat accurate value for the exact t:
  for (int i = 0; i < B_parts; i++) {
    if (param <= t_distMap[i]) {
      tMin = (i - 1) / B_parts;
      tMax = i / B_parts;
      lenMin = t_distMap[i - 1];
      lenMax = t_distMap[i];

      t = (param - lenMin) / (lenMax - lenMin) * (tMax - tMin) + tMin;
  return t;

//Creates a vector of length 1 which shows the direction from B to A
PVector_D createTangent(PVector_D pointA, PVector_D pointB) {
  PVector_D tangent = PVector_D.sub(pointA, pointB);
  return tangent;

PVector_D [][] makeZeroVectors(int x) {
  PVector_D [][] zs = new PVector_D[x][2];
  for (int i = 0; i < x; i++) {
    zs[i][0] = new PVector_D(0, 0);
    zs[i][1] = new PVector_D(0, 0);
  return zs;

PVector_D squareVector(PVector_D p) {
  return new PVector_D (p.x*p.x, p.y*p.y);

double summaXY(PVector_D p) {
  return p.x+p.y;

PVector_D bezierQ(PVector_D [] ctrlPoly, double t) {
  double tx = 1.0 - t;
  PVector_D pA = PVector_D.mult(ctrlPoly[0], tx * tx * tx); 
  PVector_D pB = PVector_D.mult(ctrlPoly[1], 3 * tx * tx * t);
  PVector_D pC = PVector_D.mult(ctrlPoly[2], 3 * tx * t * t); 
  PVector_D pD = PVector_D.mult(ctrlPoly[3], t * t * t);
  return PVector_D.add(PVector_D.add(pA, pB), PVector_D.add(pC, pD));

PVector_D bezierQPrime(PVector_D [] ctrlPoly, double t) {
  double tx = 1.0 - t;
  PVector_D pA = PVector_D.mult(PVector_D.sub(ctrlPoly[1], ctrlPoly[0]), 3 * tx * tx);
  PVector_D pB = PVector_D.mult(PVector_D.sub(ctrlPoly[2], ctrlPoly[1]), 6 * tx * t); 
  PVector_D pC = PVector_D.mult(PVector_D.sub(ctrlPoly[3], ctrlPoly[2]), 3 * t * t);
  return PVector_D.add(PVector_D.add(pA, pB), pC);

PVector_D bezierQPrimePrime(PVector_D [] ctrlPoly, double t) {
  return PVector_D.add(PVector_D.mult(PVector_D.add(PVector_D.sub(ctrlPoly[2], PVector_D.mult(ctrlPoly[1], 2)), ctrlPoly[0]), 6 * (1.0 - t)), PVector_D.mult(PVector_D.add(PVector_D.sub(ctrlPoly[3], PVector_D.mult(ctrlPoly[2], 2)), ctrlPoly[1]), 6 * t));

It sounds like you should simplify/de-densify the curve first, and then run it through spliner. PGS can perform these operations. You may also want to look at its gaussian smoothing offering:


And what exactly are the images here? The output from Schneider algorithm (where colour denotes a different curve)?

And are you looking for a “tighter” fit?

1 Like

Thanks! The library and its smoothing offering look promising! I will check them out on the Weekend.

The images I posted were made by your algorithm, that you shared in February. They are not precise enough for me, since I want to make a typeface, where the details matter a lot.

The images that I managed to produce by the code from my previous post look like this:

These four images have different values of passable error. but even the roughest (the last) one still has straight segments here:

That’s why I think I did something wrong.

Ultimately, I want to receive in the end a contour fitted as tight as possible, but less rough and noisy. Like below.

This last picture I made with the Adobe Illustrator built-in function “Simplify Path”.

P.S. — on all the pictures source path is thin and black, and the red one is an output