# Implementing Hobby Curve

I love programming. When you love something, and find something adorable, you want to hold it in your own hands. And when I find some adorable programming idea, I have to implement my own. For example, what is the fun of printing `hello world`

if not for the fact that you wrote your own? I never used *metafont* nor its cousin *metapost*. However, when I read its manual and learned its idea, I have to learn to implement my own version. In particular, I want to implement the "Hobby Curve".

Hobby curve is a spline interpolation algorithm. Given a set of points, we want to find a curve that goes through all the points and are smooth. Practically, we want to connect the dots with curve segments that can be described with cubic polynomials. Such class of problems are called cubic spline interpolation.

Each cubic segment can be defined with 4 parameters:

To solve an equation with 4 unknowns, we need 4 equations (or 4 boundary conditions). The two end points provide 2 equations. The spline has to be smooth across all points, which in math means the first order derivative (tangent at each point) has to be continuous. This first order smoothness provides the 3rd boundary condition. Three equations with 4 unknowns have infinite number of solutions. We don't like undetermined behavior, so we need one more condition to fix the solution. *Naturally*, mathematician say, let's make the curve continuous to the 2nd order derivative. The curve yielded by this second order smoothness is called *natural spline*. However, more smooth does not means more pleasant. In this case, Professor Knuth and his student Hobby found that to be pleasant, we want to minimize the overall *curvature* (bend per arc length), even at the cost of 2nd order smoothness.

Following are the difference:

I have to say the latter appears more *natural* than the natural cubic spline.

### Bézier Curve

Rather than to describe the cubic segment with *a, b, c, d* as we learned in middle school, it is geometrically more intuitive to use four points as the four parameters -- two end points and two control points:

*1-t = t = 1/2*, which leads to the famous bisection property, which leads to an efficient rasterization algorithm for such curves.

Let's introduce some terms:

θ and φ are relative curvature angles (φ is negative as shown in the figure). ρ and σ are relative velocity (normalized by the distance between the two end points). Further we'll introduce two tension parameters (τ1 and τ2) so we'll have the control points as:

Depending on θ and φ, there are optimal functions of ρ and σ that minimizes the overall curvature. Unfortunately, such minimization are neither easy to compute nor stable. What Knuth and Hobby did was to approximate the actual optimized function with an approximate one. In **metafont**, following equations were chosen:

*a, b, c*were empirically chosen as

### Solving the curve

With above definitions, we are essentially solving for θ and φ (the control points are determined from θ and φ using velocity function, which approximately minimizes the curvature). The neighboring θ and φ need satisfy the continuity (which provides 1st order smoothness). For another equation, we require the curvature to be continuous. An approximate formula for curvature (mock curvature) is used:

### Write some code

Until we have some working code to play with, the concept will remain vague. Let's write some code! To make a easier read, we'll implement in Perl -- much shorter compared to C.

```
sub solve_path($path){
...
}
```

Let's assume `$path`

is an array reference which holds the list of points, each contains attributes of "x" and "y" coordinates as well as necessary information such as "tension" or "curl". In Perl, it is represented as a list of hash reference. We'll implement the function `solve_path`

, which upon finish, will populate each points with their control points as attributes of "x-", "y-", "x+" and "y+". I'll use MyDef, a macro like meta-layer, read it as pseudo-Perl code.

```
fncode: solve_path($path)
# setup a macro for accessing attributes
$(set:p=$path->[$1]->{"$2"})
my $n=$#$path
my $N=$n
$if $(p:0,cycle)
# if cyclic, points go round, padding the list makes loops easy
push @$path, $path->[0]
push @$path, $path->[1]
$N+=2
# d is the distance between points, psi is the turning angle between segments
$call
```**calc_d_and_psi**
# solving theta angles at each points
my (@theta)
$if $(p:0,cycle)
$psi[$N]=$psi[1]
$call **solve_cycle**
$else
$psi[$N]=0
$call **solve_open**
# calculate control points from theta using velocity function
$call **set_controls**
# clean up the cyclic padding
$if $(p:0,cycle)
$(p:0,x-)=$(p:$n+1,x-)
$(p:0,y-)=$(p:$n+1,y-)
pop @$path
pop @$path

This is the skeleton of the entire function. We'll dive straight into `solve_open`

and `solve_cycle`

. `calc_d_and_psi`

and `set_controls`

are rather trivial, which we will simply append in the end.

### Solving tri-diagonal linear equations

To make it simple, we could simply just populate a matrix then call a library to solve the equations. However, tri-diagonal form of linear equations is rather easy to solve, and the process also illustrates the idea of Gaussian elimination beautifully. Not to impose for the credit, I basically copied the implementation from Knuth's metafont code. To make reference easier, I try to keep consistent variable names and equation forms. For example, rather than to use tension directly, Knuth uses α and β to refer to the *reciprocal* of left and right tension. We'll stick to the same convention and provide reference to Knuth's book where necessary.

The basic idea is to use Gaussian elimination to reduce tri-diagonal equations into di-diagonal equations, then diagonal form, which solves the equation. Recall the general form:

If the curve is open (non-cyclic), the end equation only contains two unknown, which can be used to write each equation into the form:

I am keeping the explanation to the minimum -- just to provide necessary reference to the variable names used in code. To really understand math, one has to take his own time. Once you got it, you'll find it easy (but still leave other people puzzled no matter how you try to explain).

```
# ---- open -----------------------------------------
# - t0 + u0 t1 = v0
# - t1 + u1 t2 = v1
# - ...
subcode: solve_open
my (@u, @v)
$call
```**check_ends**
$call **start_open**
$call **build_eqns**
$call **end_open**
# ----------------
subcode: check_ends
my ($dir0, $dirn)
$if defined $(p:0,dir+x)
$dir0=atan2($(p:0,dir+y), $(p:0,dir+x))
$elif defined $(p:0,x-)
$dir0=atan2($(p:0,y)-$(p:0,y-), $(p:0,x)-$(p:0,x-))
$if defined $(p:$n,dir-x)
$dirn=atan2($(p:$n,dir-y), $(p:$n,dir-x))
$elif defined $(p:$n,x+)
$dirn=atan2($(p:$n,y+)-$(p:$n,y), $(p:$n,x+)-$(p:$n,x))
subcode: start_open
$if defined $dir0
$u[0]=0
$v[0]=reduce_angle($dir0-atan2($dy[0], $dx[0]))
$else
# -: (a0c0+3-b0)t0 + ((3-a0)c0+b1)t1 = -((3-a0)c0+b1)p1
$call **get_curl,** $curl, $(p:0,curl+)
$call **get_recip,** $a, $(p:0,tension+)
$call **get_recip,** $b, $(p:1,tension-)
my $c = $a**2*$curl/$b**2
$u[0]=((3-$a)*$c+$b) / ($a*$c+3-$b)
$v[0]=-$u[0]*$psi[1]
subcode: end_open
$if defined $dirn
NOOP
$theta[$n]=reduce_angle($dirn-atan2($dy[$n-1], $dx[$n-1]))
$else
$call **get_curl,** $curl, $(p:$n,curl-)
$call **get_recip,** $a, $(p:$n-1,tension+)
$call **get_recip,** $b, $(p:$n,tension-)
my $c = $b**2*$curl/$a**2
$u[$n] = ($b*$c+3-$a)/((3-$b)*$c+$a)
# - we have: t(n-1) + u(n-1)t(n) = v(n-1)
# - now add: t(n-1) + u(n) t(n) = 0
$theta[$n] = $v[$n-1]/($u[$n-1]-$u[$n])
# Solves all the unknowns like pulling a tread, beautiful!
$for my $i=$n-1;$i>=0;$i--
$theta[$i] = $v[$i]-$u[$i]*$theta[$i+1]

For the open curve we need establish the end equations. If the end is given directions, it is rather trivial (simply set's the angle directly). Otherwise, we need establish the equation between the first two points (or the last two points) using the curl relationship (by default, set the curvature equal). It takes some effort to derive the equation, or you could simply copy down the formula (which I listed in the comments).

In `end_open`

, once we solved the last θ using the end equation, all the other angles are solved with straight substitution -- in a single loop statement -- such is the beauty of the equation form we adopted.

We left the `build_eqns`

later, as it is shared with solving cyclic curves.

```
# ---- cycle ---------------
# -: t0 + u0 t1 = v0 + w0 t0
# -: t1 + u1 t2 = v1 + w1 t0
# -: ...
subcode: solve_cycle
my (@u, @v, @w)
$call
```**start_cycle**
$(set:cycle=1)
$call **build_eqns**
$call **end_cycle**
subcode: start_cycle
# -: t0 = t0
$u[0]=0
$v[0]=0
$w[0]=1
subcode: end_cycle
# There is no end equation, but by pulling two ends together
# we could solve theta[0]:
# -: ti + ui t(i+1) = vi + wi t0
# add: t(n+1) = t0
# solve: ti = a + b t0
my $a = 0
my $b = 1
$for my $i=$n;$i>0;$i--
$a = $v[$i] - $a * $u[$i]
$b = $w[$i] - $b * $u[$i]
# -: t1 = a + b t0
# recall: t(n+1) + u(n+1)t(n+2) = v(n+1) + w(n+1) t0
# ----->: t0 + u(n+1) t1 = v(n+1) + w(n+1) t0
my $t0= ($v[$n+1]-$a*$u[$n+1])/(1-($w[$n+1]-$b*$u[$n+1]))
# -- sustitute theta[0] into all equations
$v[0]=$t0
$for $i=1:$n+1
$v[$i] += $w[$i]*$t0
# -- realize the points go round
$theta[0]=$t0
$theta[$n+1]=$t0
# -- same as before
$for my $i=$n;$i>0;$i--
$theta[$i] = $v[$i]-$u[$i]*$theta[$i+1]
$theta[$n+2]=$theta[1]

Now the code that we have left:

```
subcode: build_eqns
$for $i=1:$N
$call
```**get_recip,** $a0, $(p:$i-1,tension+)
$call **get_recip,** $a1, $(p:$i,tension+)
$call **get_recip,** $b1, $(p:$i,tension-)
$call **get_recip,** $b2, $(p:$i+1,tension-)
my $A = $a0/$b1**2/$d[$i-1]
my $B = (3-$a0)/$b1**2/$d[$i-1]
my $C = (3-$b2)/$a1**2/$d[$i]
my $D = $b2/$a1**2/$d[$i]
my $t = $B-$u[$i-1]*$A+$C
$u[$i]=$D/$t
$v[$i]=(-$B*$psi[$i]-$D*$psi[$i+1]-$A*$v[$i-1])/$t
$(if:cycle)
$w[$i]=-$A*$w[$i-1]/$t

The "A, B, C, and D" here are not exactly those we listed previously (as in the nominal equation form). They are exactly as referenced by Knuth, where each equation has the form:

You could derive them, or trust them.

The rest of the trivial but necessary parts:

```
# ----------------------------
subcode: calc_d_and_psi
my (@dx, @dy, @d, @psi)
$for $i=0:$N
$dx[$i]=$(p:$i+1,x)-$(p:$i,x)
$dy[$i]=$(p:$i+1,y)-$(p:$i,y)
$d[$i]=sqrt($dx[$i]**2+$dy[$i]**2)
$for $i=1:$N
my $t_sin = $dy[$i-1]/$d[$i-1]
my $t_cos = $dx[$i-1]/$d[$i-1]
my $x = $dx[$i]*$t_cos + $dy[$i]*$t_sin
my $y = -$dx[$i]*$t_sin + $dy[$i]*$t_cos
$psi[$i]=atan2($y, $x)
# ---------------------------
subcode: set_controls
$for $i=0:$N
my $phi = -$psi[$i+1]-$theta[$i+1]
$call
```**get_tension**
$call **get_sin_cos**
$call **get_velocity**
$(p:$i,x+)=$(p:$i,x)+($dx[$i]*$ct-$dy[$i]*$st)*$rho
$(p:$i,y+)=$(p:$i,y)+($dx[$i]*$st+$dy[$i]*$ct)*$rho
$(p:$i+1,x-)=$(p:$i+1,x)-($dx[$i]*$cf+$dy[$i]*$sf)*$sigma
$(p:$i+1,y-)=$(p:$i+1,y)-(-$dx[$i]*$sf+$dy[$i]*$cf)*$sigma
subcode: get_tension
$call **get_recip,** $a, $(p:$i,tension+)
$call **get_recip,** $b, $(p:$i+1,tension-)
subcode: get_sin_cos
my $st=sin($theta[$i])
my $ct=cos($theta[$i])
my $sf=sin($phi)
my $cf=cos($phi)
subcode: get_velocity
$(set:a=1.41421356)
$(set:c=0.38196601125)
$(set:cc=0.61803398875)
my $alpha=$(a)*($st-$sf/16)*($sf-$st/16)*($ct-$cf)
my $beta = (1+$(cc)*$ct+$(c)*$cf)
my $rho=(2+$alpha)/$beta
my $sigma=(2-$alpha)/$beta
$rho *= $a/3
$sigma *= $b/3
# ----------------------
fncode: reduce_angle($t)
$(set:Pi=3.14159265)
$(set:TwoPi=$(Pi)*2)
$if abs($t)>$(Pi)
$if $t>0
$t-=$(TwoPi)
$else
$t+=$(TwoPi)
return $t
# ----------------------
subcode: get_curl(a, var)
my $(a)=1
$if defined $(var)
$(a)=$(var)
subcode: get_recip(v, tension)
my $(v)=1
$if defined $(tension)
$(v)=1/$(tension)

For the velocity part, reference the equation earlier.

In summary, the entire algorithm is simply solving linear equations with Gaussian elimianation. In code, it is rather complex, but ideally, one should carry out the algebra on paper. Then for the code, it is merely trascribing the derived algebra.

### Demo

I used the above code as part of my own pdf plotting library. In addition, I also implemented a MyDef `output_plot`

module where it employs some metafont DSL:

```
page: test, basic_frame, t.pdf
$set_point i=0:6, a=100, origin (300,300)
j = i+1
xj = (i%3-1) * a
yj = (0.5-i/3) * a
# I need invent some DSL to make setting graphic states easier, maybe CSS?
$call
```**linewidth,** 20
$call **linecolor,** 0.5
$call **setdraw,** stroke
$draw z5..z4..z1..tension 1.2..z3..z6..cycle
MyPlot::draw_points($page_content, \@path)

With `mydef_page -mplot test.def`

, it generates `test.pl`

, running it produces `t.pdf`

:

For reference, `test.pl`

:

```
use strict;
use MyPlot;
our $page_content;
{
MyPlot::init();
$page_content=MyPlot::newpage();
MyPlot::line_width($page_content, 2);
MyPlot::line_cap($page_content, "round");
MyPlot::line_join($page_content, "round");
my ($x1, $y1) = (200, 350);
my ($x2, $y2) = (300, 350);
my ($x3, $y3) = (400, 350);
my ($x4, $y4) = (200, 250);
my ($x5, $y5) = (300, 250);
my ($x6, $y6) = (400, 250);
MyPlot::line_width($page_content, 20);
MyPlot::stroke_color($page_content, "0.5");
my @path;
push @path, {"cycle"=>1, "x"=>$x5, "y"=>$y5};
push @path, {"x"=>$x4, "y"=>$y4};
push @path, {"tension+"=>1.2, "x"=>$x1, "y"=>$y1};
push @path, {"tension-"=>1.2, "x"=>$x3, "y"=>$y3};
push @path, {"x"=>$x6, "y"=>$y6};
MyPlot::solve_path(\@path);
MyPlot::do_path($page_content, \@path, 0);
MyPlot::close_path($page_content);
push @$page_content, "S";
MyPlot::draw_points($page_content, \@path);
MyPlot::write_pdf("t.pdf");
print " --> t.pdf\n";
}
```

### Carry the idea forward

In programming, we have always been reusing the code via library or as application. A library or an application is not one idea, it is rather many ideas integrated together. When ideas are bundled together, it is complex and it is often beyond our ability to comprehend, and it will not easily be reused. Within Knuth's TeX and Metafont there are many good ideas, many of which address problems we are still troubling with today. However, as long as they are bundled, we do not see the solution goes beyond TeX or Metafont, despite the fact that both of them are written with *literate programming*.

In this post, on the surface I may simply reinvented a wheel that by all means less elegant than the original. However, my effort isolates the idea from its complex context. Now that I have understood the idea (via implementing my own), I will carry these cute Metafont DSL wherever I go, whichever language I use, on what ever problem I encounter. That is the idea.