useless rants

Bézier curves in Haskell part 2

You can find the source for this article here.

Let be a Bézier curve over a polygon , and let be its derivative. To construct , we must first choose an origin . As a result we are now working in a vector space. We define a polygon where . Then is the Bézier curve over , up to an homothety of ratio .

> bezierDerivative1 :: (AffineSpace p, AffineSpace (Diff p),
>                       VectorSpace (Diff p), VectorSpace (Diff (Diff p))) =>
>                      Seq p -> Scalar (Diff (Diff p)) -> Diff p
> bezierDerivative1 polygon t =
>     let n = length polygon - 1
>         polygon' = fmap (uncurry (.-.)) $ pairs polygon
>         bezier1 = bezierSeq polygon' t in
>     sumV (take n (repeat bezier1))
>
> pairs :: Seq a -> Seq (a, a)
> pairs = fst . foldl toPair (empty, Nothing)
>     where toPair (ps, Nothing) p0 = (ps, Just p0)
>           toPair (ps, Just p0) p1 = ((p0, p1) <| ps, Just p1)

Now that's some serious type level fun. Here the AdditiveGroup module doesn't provide functions for multiplying by integers, so I had to roll my own. The pairs function just pairs up the elements of the sequence into a sequence of pairs. Let's try this with an example:

> examplePolygon :: Seq (Double, Double)
> examplePolygon = fromList [(1, 1), (-3, -5), (-3, 5), (1, -1)]

So we have this curve for the Bézier curve over examplePolygon:

and its first derivative:

Now how do we generalize this to the th derivative of a Bézier curve, where is a natural number? As the derivative is also a Bézier curve, we can apply the same reasoning to it. The th derivative of is the Bézier curve over the polygon defined by the differences of the control points of the th derivative, up to an homothety of ratio .

So we first define a function giving us the list containing the polygon for each derivative:

> polygon' :: (AdditiveGroup p) => Seq p -> [Seq p]
> polygon' = iterate (fmap (uncurry (^-^)) . pairs)

And we use it to define our th derivative easily:

> bezierDerivative :: (AdditiveGroup p, AffineSpace p, VectorSpace (Diff p)) 
>                  => Int -> Seq p -> Scalar (Diff p) -> p
> bezierDerivative r polygon t =
>     let n = length polygon - 1
>         polygon_r = polygon' polygon !! r
>         bezier_r = bezierSeq polygon_r t in
>     sumV (take (product [r+1..n]) (repeat bezier_r))

So we have, for over examplePolygon:

I think there are many improvement to be done here, most notably the final homothety is quite ugly and inefficient. I'll try to find a better solution.

Next time, we're going to look at B-spline curves.

Bézier curves in Haskell part 1

This article was intended to be a literate Haskell program, but WordPress decided otherwise. Oh well. You can still download the compressed source file here.

I'm currently studying Bézier curves, so I figured I would implement what I know.

I know of two algorithms for constructing Bézier curves: De Casteljau's algorithm and using Bernstein polynomials. The algorithm using Bernstein polynomials is more efficient, but propagates more errors during floating point operations. De Casteljau's algorithm also provides wonderful geometric insight into the nature of Bézier curves, so this is the one I'm going to use.

You can construct a Bézier curve over any affine space, and the AffineSpace typeclass allows us to work in complete generality, so this is what we're going to do.

bezier, bezier' :: (AffineSpace p, VectorSpace (Diff p)) => 
                   [p] -> Scalar (Diff p) -> p

So far it's only a bit of type gymnastic in order to work with affine spaces.

De Casteljau's algorithm goes as follows: let a polygon in an affine space, and let be the bézier curve over . To compute where , we first compute the points for . We then iterate this process over the points we get, until we get only one point which is .

This algorithm is very imperative in nature, so we will derive a recursive algorithm from it. If we look at the formula for the algorithm gives us, we notice that the final affine linear interpolation occurs over two other Bézier curves evaluated at . And we also notice that a Bézier curve over a single point is this same point.

So our first version simply does recursion over the polygon, exploiting the fact that a Bézier curve over can be expressed in function of the Bézier curve over and the Bézier curve over :

bezier' [p] _ = p
bezier' polygon t = 
    let poly0 = reverse . tail $ reverse polygon
        poly1 = tail polygon in
    alerp (bezier' poly0 t) (bezier' poly1 t) t

Here alerp p p' t varies from p to p' as t varies from 0 to 1. We test it with an example polygon:

examplePolygon :: [(Double, Double)]
examplePolygon = [(0, 0), (1, 1), (2, 0)]
 
plotBezier f = do 
    let i = [0,0.001 .. 1]
        points = map f i
        (xs, ys) = unzip points
    plotWindow xs ys

And then running the following in ghci:

plotBezier (bezier' examplePolygon)

gives us the following result:

Now, bezier' is not without fault, as we must remove the last element of the list, which is highly inefficient for large lists (hence the ugly double reverse, which is in anyway). In practice, Bézier curves are often used with small polygons and then attached to other low degree bezier curves to form larger curves, so this is not a big deal for this use case. But we cannot let such inefficiency pass by innocently.

Basically, we need a data structure which allows efficient removal at the beginning or the end of the sequence. This is precisely what the Seq data structures offers:

bezier polygon t =
    bezierSeq (fromList polygon) t
 
bezierSeq :: (AffineSpace p, VectorSpace (Diff p)) => 
             Seq p -> Scalar (Diff p) -> p
bezierSeq polygon t =
    let poly0 :> _ = viewr polygon
        p :< poly1 = viewl polygon in
    if null poly1
    then p
    else alerp (bezierSeq poly0 t) (bezierSeq poly1 t) t

The Seq datatypes guarantees constant-time removal of the first or last element. Otherwise the algorithm is precisely the same as in bezier'.

Next time, we will look at derivatives for Bezier curves.

Prüfer sequence and Haskell

A Prüfer sequence is a very compact data structure for a labeled tree. For a tree with n nodes, it only takes a list of n - 2 integers to describe it fully.

Trees in maths are usually represented as T = (X, A), where X is a set of nodes, and A is a set of couples (n1, n2) representing an edge between the nodes n1 and n2. There is a simple, iterative algorithm to convert a tree from the T = (X, A) form to the corresponding Prüfer sequence, which is what we are going to implement in Haskell.

X = {1,2,3,4,5,6}
A = {(1,5),(2,5),(3,5),(5,4),(4,6)}
Prüfer sequence: (5,5,5,4)

First, we define the data type for trees:

data Tree a = Tree [a] [(a, a)]
              deriving Show

The first line defines a new data type called Tree. This is kind of similar to a structure definition in C. [a] stands for a list containing elements of type a, and (a, a) stands for a pair of types a, therefore [(a, a)] stands for a list of pairs of types a. Here, 'a' stands for 'any type', so Tree is basically what java-heads would call a generic class, where 'a' is '<T>'.

As you may have guessed, [a] stands for the X set, and [(a, a)] for the A set. Don't pay too much attention to the second line, it simply allows us to display the contents of a Tree easily.

Let's take a look at the algorithm to convert a tree into a Prüfer sequence:

function prufer (tree T):

S = the empty set
C = the empty list
while T has more than 2 nodes left:

j = the leaf with the smallest number
k = the only node adjacent to j
add j to S
add k to S
remove j from T

end while
return C

This algorithm looks quite haskell-unfriendly at first glance. It has sequential operations, traditional iteration and mutable state, all of which don't exist in haskell. Thankfully, it's quite easy to convert such a simple iterative algorithm into a recursive algorithm. But before diving straight into it, it is a good idea to define some functions for manipulating trees:

size :: Tree a -> Int
size (Tree a _) = length a

Size returns the number of nodes in the tree. The first line may look quite scary at first, but it's only a type signature for the function. It means that the function size takes a Tree of any type, and returns an integer. It's very similar to the notation in maths.

The second line is fairly straightforward. It simply returns the length of the A list of the tree.

contains :: (Eq a) => a -> (a, a) -> Bool
contains x (y, z) =
    x == y || x == z

This function has a slightly more complicated type signature than the previous one. '(Eq a) =>' is a type constraint, so the function only accept types that have the '==' operator defined for them.

'a -> (a, a) -> Bool' means that the function takes one parameter of a type 'a', a pair of types 'a', and returns a boolean value. So contains takes a value and a pair, and returns true if the pair contains the value, false otherwise. The body of the function itself is self-explanatory if you have had any experience with a C-like language.

adjacent :: (Eq a) => Tree a ->a -> a
adjacent (Tree a x) n =
    let pair = head (filter (contains n) x) in
    if fst pair == n
    then snd pair
    else fst pair

The adjacent function is a bit more complicated. It takes a tree and a node n, and returns a node adjacent to n.
The expression 'let x = expr1 in (expr2 that uses x)' simply replaces x in expr2 by expr1, kind of like a macro in C/C++. It's called let-binding, and it's very useful for making the code clearer and more concise.

The filter function is very straightforward to use: filter p l returns all the elements of the list l that verify the predicate p. Here, it returns the list of pairs from the X set containing n, which in natural language means the full list of the edges adjacent to n.
Head then simply takes the first element of the list, as we only need one node.

The following if p then x else y sequence should be quite clear, you just need to know that fst and snd return the first and second value from a pair respectively.

degree :: (Eq a) => Tree a -> a -> Int
degree (Tree _ x) n =
    foldr (+) 0 (map (\y -> if contains n y then 1 else 0) x)

The degree function returns the number of edges adjacent to a node in the graph. It uses 2 other list processing functions, somewhat similar to filter.

The map function take a function f and a list l, applies f to every element of the list, then returns the new list. Here we use a lambda expression to define the function to apply. A lambda expression is a quick way to define anonymous functions. It's actually quite simple: '\y -> x' is the equivalent to the mathematical notation .

Here, the map function turns the list of edge into a list of Integers, 1 if the edge is adjacent to the node n, 0 otherwise. Afterwards, we use the foldr function: it basically replaces the commas in the list by the operator given, so 'foldr (+) 0 [1,2,3,4]' processes 0+1+2+3+4 = 10. Now it should be quite clear how this function works :).

leaves :: (Eq a) => Tree a -> [a]
leaves (Tree a x) =
    filter (\y -> degree (Tree a x) y == 1) a

The leaves function returns the list of the leaves in the tree. If you're unfamiliar with graph theory, a leaf is a node with only one other adjacent node, which also means its degree is 1. Here we use the filter function to return only the elements which have a degree of 1. Nothing new here.

rm :: (Eq a) => Tree a -> a -> Tree a
rm (Tree a x) n =
    Tree (filter (/= n) a) (filter (not . contains n) x)

The rm function removes a node n from the given tree, then returns the new tree.
So we return a Tree with his A set filtered so n is not there anymore, and it's X set filtered so there are no edges containing n anymore. The trick to the last part is to use the dot '.' operator, which is the function composition operator in haskell, exactly like in maths.
This particular function should be extremely intuitive to understand :).

Now that we wrote all the building blocks, we can finally write the function to convert from a regular tree to a Prüfer sequence:

prufer :: Tree Int -> [Int]
prufer tree =
    prufer' tree []
    where prufer' t s | size t == 2     = []
                      | otherwise       =
                          let j         = minimum (leaves t)
                              k         = adjacent t j
                              new_t     = rm t j
                              new_s     = s ++ [j] in
                          [k] ++ prufer' new_t new_s

The type signature means our function takes a tree containing integers, and return a list of integers, which corresponds to the C list. Note that we didn't specify integers in our previous functions because it wasn't necessary, here we need a tree of n nodes with its nodes numbered from 1 to n.

We start the function by calling a helper function, prufer', which takes a tree of integers t and the s set. The s set is empty at the beginning, so we give it the empty list '[]'.

We then use a where clause to define the prufer' function. It works exactly like 'let ... in', except after the expression. The prufer' function is a recursive function, so we define a base case where the function will stop recursing, and a general case where the function will keep on recursing. For that we use a '|' guard to define a condition, exactly like an 'if ... then ... else' block would.

In the base case, when the size of t is lower or equal than 2, we stop the recursion by returning the empty list. In the general case, which is every other case, we start by defining the values we are going to use thanks to let-binding:

• j, the leaf of t with the smallest value
• k, the only node adjacent to j
• new_t, the current tree t but without j
• new_s, the current s list with j added (++ is the operator for concatenating lists)

As you can see, our previously defined functions do all the job for us, no need for things to be overcomplicated :).
Then, we simply add k to the C list, and call the function again on the new tree and the new S set.

And it's done! Open the file with your favorite haskell interpreter and see for yourself:

*Main> test
Tree [1,2,3,4,5,6] [(1,5),(3,5),(2,5),(5,4),(4,6)]
*Main> prufer test
[5,5,5,4]
*Main>

Here is the source code if you want to look at it more closely, I even included the function to turn a Prüfer sequence back into a regular tree, but I won't detail its inner workings here, 1500~ words are more than enough for such a boring topic :).