Created 2022-03-30

Introduction

We're all used to using an image format for pictures and SVGs for lines and vector images. Now, JPEG XL can do both, kinda.

On this page, I'll be discussing Cubic Bezier curves and Centripetal Catmull-Rom splines. Also, the ideas and methods I use here are quite basic and should be easy to follow. It being too simple and boring could also be a reason it hasn't been explored in depth.

There's a big difference between a Bezier curve and a Catmull-Rom spline. One is a type of parametric curve and the other is a type of cubic spline interpolation. The four control points used to represent a Bezier curve can't be used to create the four control points of a Catmull-Rom spline. Yes, I've read this and this, the top two search results on the topic. They do not work.

Maybe saying they don't work is a bit inaccurate. They do perfectly recreate the curves as shown here, but these go beyond the start and end point of the desired spline. This is when I realized that Catmull-Rom splines is just Catmull-Rom interpolation between points. This may be familiar if you've used the line tool in paint.net or the curvature pen tool in Photoshop. Anyway, let's stop thinking too hard about this since all that's needed to be done is to sample the Bezier curve to get those points.

Let's first define the Bezier curve.

\(\mathbf{B}(t) = (1 - t)^3 \mathbf{P}_0 + 3t(1 - t)^2 \mathbf{P}_1 + 3t^2(1 - t) \mathbf{P}_2 + t^3 \mathbf{P}_3, 0 \leq t \leq 1\)

def bezier1d(p0, p1, p2, p3, t):
  return p0 * (1 - t)**3 + 3 * p1 * t * (1 - t)**2 + 3 * p2 * t**2 * (1 - t) + p3 * t**3
def bezier2d(p0, p1, p2, p3, t):
  return (
    bezier1d(p0[0], p1[0], p2[0], p3[0], t),
    bezier1d(p0[1], p1[1], p2[1], p3[1], t),
  )

We also want to use equidistant points. We can achieve this by using the first derivative to get the size of each step. This is an approximation of the arc length parameterization.

\(\mathbf{B}'(t) = 3(1 - t)^2(\mathbf{P}_1 - \mathbf{P}_0) + 6t(1-t)(\mathbf{P}_2 - \mathbf{P}_1) + 3t^2(\mathbf{P}_3 - \mathbf{P}_2)\)

def bezier1d_dt(p0, p1, p2, p3, t):
  return 3 * (1 - t)**2 * (p1 - p0) + 6 * t * (1 - t) * (p2 - p1) + 3 * t**2 * (p3 - p2)
def bezier2d_dt(p0, p1, p2, p3, t):
  return (
    bezier1d_dt(p0[0], p1[0], p2[0], p3[0], t),
    bezier1d_dt(p0[1], p1[1], p2[1], p3[1], t),
  )

Now we can sample the curve.

def sample_bezier(p0, p1, p2, p3):
  points = []
  t = 0
  while t < 1:
    if t > 1: break
    points.append(bezier2d(p0, p1, p2, p3, t))
    dt = bezier2d_dt(p0, p1, p2, p3, t)
    l = length(dt[0], dt[1])
    t += 1 / l # 1 is the absolute distance per step
  if t > 1:
    points.append(bezier2d(p0, p1, p2, p3, 1))
  return points

A Bezier curve with points (62, 168), (193, 18), (8, 20), (130, 170) on a 200x200 canvas:

Sampling with straight line segments:

Creating the Catmull-Rom spline

As mentioned in the introduction, these points can be used without any further modification.

Firstly, a point:

struct Point {
  double x;
  double y;
  Point operator+(const Point &a) { return {x + a.x, y + a.y}; }
  Point operator-(const Point &a) { return {x - a.x, y - a.y}; }
  Point operator+(const double &a) { return {x + a, y + a}; }
  Point operator-(const double &a) { return {x - a, y - a}; }
  Point operator*(const double &a) { return {x * a, y * a}; }
};

The definition of the Catmull-Rom spline is quite lengthy so I will only show here the simplified code:

Point catmull_rom(Point points[4], double t) {
  double alpha = 0.5;
  Point p0 = points[0];
  Point p1 = points[1];
  Point p2 = points[2];
  Point p3 = points[3];
  Point a = p1 * 3.0 - p2 * 3.0 + p3 - p0;
  Point b = p0 * 2.0 - p1 * 5.0 + p2 * 4.0 - p3;
  Point c = p2 - p0;
  Point d = p1 * 2.0;
  return (a * pow(t, 3) + b * pow(t, 2) + c * t + d) * alpha;
}

The same points input into my Catmull-Rom visualizer:

It works! Obviously, there's way too many points. At some resolution, there will be points which are redundant.

Optimizing Catmull-Rom splines

This part actually took more brain power. Not because it was hard, but because it was confusing to think of.

The method I will use to optimize is to remove a point and calculate the error. The error function can be defined as the area between the old and new spline. To do this, I will sample the two lines and then run the shoelace algorithm to find the area. If the area is smaller than a set error, I will delete the point.

This method is a bit greedy so I will store those errors and start by deleting the one with the least error.

Catmull-Rom:

Sampling:

When sampling to calculate the area, I'll collect double the samples between the second and third point in the new line. The first half of these samples will be matched with the points between the second and third point in the new line and the second half will be matched with the points between the third and fourth point in the new line.

tl;dr: see the JS demo below

int res = 10;
int start = 1;
int center = 2;
int end = 4;
std::vector<Point> line1_spls;
std::vector<Point> line2_spls;
void add_samples(std::vector<Point> &spls, std::vector<Point> points, int i, int count) {
  spls.push_back(points[i]);
  Point line[4] = {
    i ? points[i - 1] : points[0],
    points[i],
    points[i + 1],
    i + 2 < points.size() ? points[i + 2] : points[i + 1],
  };
  for (int j = 1; j < count; j++) {
    Point p = catmull_rom(line, (double)j / count);
    spls.push_back(p);
  }
  spls.push_back(i + 2 < points.size() ? points[i + 2] : points[i + 1]);
}
for (int i = start; i < end; i++) {
  add_samples(line1_spls, line1, i, res);
}
line1_spls.push_back(line1[end]);
for (int i = start; i < end - 1; i++) {
  add_samples(line2_spls, line1, i, i == center - 1 ? res * 2 : res);
}
line2_spls.push_back(line2[end]);

Shoelace:

double shoelace(Point p0, Point p1, Point p2, Point p3) {
  return abs(((p0.x * p1.y - p0.y * p1.x) + (p1.x * p2.y - p1.y * p2.x) +
              (p2.x * p3.y - p2.y * p3.x) + (p3.x * p0.y - p3.y * p0.x)) /
              2);
}

JS Demo

The same points before, optimized with an allowed error of 8 square pixels, scaled by 2:

SVG?

In SVG, paths are combinations of lines (2 points), cubic beizers (4 control points), quadratic beziers (3 control points), and arcs.

In Python we can use svgpathtools and xml.etree.ElementTree to get the paths and basic information on the SVG.

import svgpathtools
import xml.etree.ElementTree as ET
file = "curves.svg"
paths, attributes, svg_attribs = svgpathtools.svg2paths2(file)
tree = ET.parse(file)
root = tree.getroot()
viewBox = [float(x) for x in root.attrib["viewBox"].split(" ")]
off_x = -viewBox[0]
off_y = -viewBox[1]
try:
  width = float(root.attrib["width"])
except:
  width = viewBox[2]
try:
  height = float(root.attrib["height"])
except:
  height = viewBox[3]
scale_x = width / viewBox[2]
scale_y = height / viewBox[3]

And a short example converting a cubic bezier into points and optimizing them:

import subprocess
import optimize
splines = []
for path, attribute in zip(paths, attributes):
  current_spline = []
  for segment in path:
    if type(segment) == svgpathtools.CubicBezier:
      p0 = (segment.start.real + off_x, segment.start.imag + off_y)
      p1 = (segment.control1.real + off_x, segment.control1.imag + off_y)
      p2 = (segment.control2.real + off_x, segment.control2.imag + off_y)
      p3 = (segment.end.real + off_x, segment.end.imag + off_y)
      p0 = (p0[0] * scale_x, p0[1] * scale_y)
      p1 = (p1[0] * scale_x, p1[1] * scale_y)
      p2 = (p2[0] * scale_x, p2[1] * scale_y)
      p3 = (p3[0] * scale_x, p3[1] * scale_y)
      cubic_points = sample_bezier(p0, p1, p2, p3)
      cubic_points = [(p[0], p[1]) for p in cubic_points]
      new_points = optimize.optimize(cubic_points, 4)
      splines.append(new_points)

Conclusion

The final points can be used with the jxl_from_tree tool. How.

Was this a waste of time? No. I'm now able to make images smaller than the original SVGs. Obviously JPEG XL is a raster image format and has set image dimensions, but SVG also technically does and upsampling can also be done with JPEG XL.

What did I learn? I didn't realize it before, but cubic splines are great for drawing. They are much easier to draw with than Bezier curves.

Relevant code and a complete(working) program can be found here.

Notes

The optimizing section is done in C++ since Python is too slow for this task.

Instead of only a point being removed, a point can be removed and surrounding points also shifted.

Another error function can use the tangents along the spline. This may perform better in keeping the general shape.

Segments in SVG paths can be continuous and able to be connected. By joining continuous segments, more points can be purged and less splines need to be created.