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:

$\dpi{110}&space;y&space;=&space;ax^3+bx^2+cx+d$

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.

$\dpi{110}&space;\kappa&space;=&space;\frac{d\phi}{ds}&space;=&space;\frac{x'y''&space;-&space;y'x''}{(x'^2+y'^2)^\frac{3}{2}}$

Following are the difference:

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

Bé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:

$\dpi{110}&space;z(t)=(1-t)^3z_1&space;+&space;3(1-t)^2tz_2&space;+&space;3(1-t)t^2z_3&space;+&space;t^3z_4$
where t varies from 0 to 1. In particular, when t=1/2, 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:

$\dpi{110}&space;\vec{z_2}&space;=&space;\vec{z_1}&space;+&space;\frac{\rho}{3\tau_1}\vec{w_1}l$
$\dpi{110}&space;\vec{z_3}&space;=&space;\vec{z_4}&space;+&space;\frac{\sigma}{3\tau_2}\vec{w_2}l$
where w_1 and w_2 are unit vectors point from end points to their control points.

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:

$\dpi{110}&space;\rho&space;=&space;f(\theta,&space;\phi)&space;=&space;\frac{2+\alpha}{\beta}$
$\dpi{110}&space;\sigma=f(\phi,&space;\theta)&space;=&space;\frac{2-\alpha}{\beta}$
where
$\dpi{110}&space;\alpha&space;=&space;a(\sin\theta&space;-&space;b\sin\phi)(\sin\phi-b\sin\theta)(\cos\theta-\cos\phi)$
$\dpi{110}&space;\beta&space;=&space;1+(1-c)\cos\theta+c\cos\phi$
The coefficients a, b, c were empirically chosen as
$\dpi{110}&space;a=\sqrt&space;2,&space;b=\frac{1}{16},&space;c=\frac{3-\sqrt&space;5}{2}$
(It is interesting that mathematician choose square roots even when the number is a an empirical heuristic approximation).

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:

$\dpi{110}&space;\hat{\kappa}(\theta,&space;\phi,&space;\tau,&space;\bar{\tau})=\tau^2(\frac{2(\theta+\phi)}{\bar{\tau}}-6\theta)$
And the continuity of mock curvature gives
$\dpi{110}&space;\hat{\kappa}(\phi_i,&space;\theta_{i-1},&space;\bar{\tau}_i,&space;\tau_{i-1})/d_{i-1}&space;=&space;\hat{\kappa}(\theta_i,&space;\phi_{i+1},&space;\tau_i,&space;\bar{\tau}_{i+1})/d_i$
Remember
$\dpi{110}&space;\theta_i+\phi_i+\psi_i=0$
where ψ is the turning angle at each point (the angle change between the two neighboring line segments). So φ can be substituted with θ, and we will have the form:
$\dpi{110}&space;A_k\theta_{k-1}+B_k\theta_{k}+C_k\theta_{k+1}=D_k$
which forms a tri-diagonal set of linear equations.

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: $\dpi{110}&space;A_k\theta_{k-1}+B_k\theta_{k}+C_k\theta_{k+1}=D_k$ We'll rewrite into the form of $\dpi{110}&space;\theta_{i}+u_{i}\theta_{i+1}&space;=&space;v_{i}&space;+&space;w_{i}\theta_0$ Each such equation essentially expresses $\dpi{110}&space;\theta_i$ in the form of $\dpi{110}&space;\theta_{i+1}$ and $\dpi{110}&space;\theta_0$ , which rolls into the next equation. 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: $\dpi{110}&space;\theta_i+u_i\theta_{i+1}=v_i$ It is one step less than the cycle curve. 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
# 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: $\dpi{110}&space;A_k\theta_{k-1}+(B_k+C_k)\theta_k+D_k\theta_{k+1}=-B_k\psi_k-D_k\psi_{k+1}$ where $\dpi{110}&space;A_k=\frac{\alpha_{k-1}}{\beta_k^2d_{k-1,k}}$ $\dpi{110}&space;B_k=\frac{3-\alpha_{k-1}}{\beta_k^2d_{k-1,k}}$ $\dpi{110}&space;C_k=\frac{3-\beta_{k+1}}{\alpha_k^2d_{k,k+1}}$ $\dpi{110}&space;D_k=\frac{\beta_{k+1}}{\alpha_k^2d_{k,k+1}}$ 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.

[X]