Providence Salumu 2D Contouring

About

The CAD systems that I design use functional representations to represent solid models. These f-reps need to be converted into meshes before they can be used in other programs (e.g. to be manufactured on a 3D printer).

There are many meshing algorithms in the literature, but we want one that is:

This write-up explores 2D implementations of a few common algorithms. We start with Marching Squares, then upgrade to Dual Contouring.

Comparison

Also, we'll be using Haskell. Hold on to your hats!

Functional representations

We'll start by defining a type for 2D points and a type for our shape functions.

type Point = (Double, Double)
type Shape = Point -> Double

A shape function takes an (x,y) tuple and and returns a value. If this value is less than or equal to zero, the point (x,y) is inside the shape.

Our first shape will be a circle:

circle :: Point -> Double -> Shape
circle (x0, y0) r (x, y) = sqrt ((x0 - x)**2 + (y0 - y)**2) - r

This function parameterizes a circle with a center and radius, then returns the distance from any point in the plane to the circle's edge. Thanks to partial application, you can call it with just (x,y) r arguments to get a Shape object out.

λ> let c = circle (0,0) 1
λ> c (0,0)      -- points inside the circle are negative
-1.0
λ> c (1,0)      -- points on the contour are zero
0.0
λ> c (2,0)      -- points outside the circle are positive
1.0

We'll also define a set of half-planes, e.g. left 1 means that all points (x,y) with x <= 1 are part of the shape:

left :: Double -> Shape
left x0 (x, _) = x - x0

right :: Double -> Shape
right x0 (x, _) = x0 - x

lower :: Double -> Shape
lower y0 (_, y) = y - y0

upper :: Double -> Shape
upper y0 (_, y) = y0 - y

Next, let's create operators for union, intersection, and negation:

(∪) :: Shape -> Shape -> Shape
a ∪ b = \p -> min (a p) (b p)

(∩) :: Shape -> Shape -> Shape
a ∩ b  = \p -> max (a p) (b p)

inv :: Shape -> Shape
inv a = \p -> -(a p)

The negation operator lets us make cutouts: subtracting b from a is given by a ∩ (inv b).

We can use these CSG operators and the half-planes defined above to create a rectangle, which is the intersection of four half-planes:

type Min = Point
type Max = Point

rectangle :: Min -> Max -> Shape
rectangle (xmin, ymin) (xmax, ymax) =
    right xmin ∩ left xmax ∩ upper ymin ∩ lower ymax

Here's a simple test shape that we'll be using for the rest of this writeup:

hi :: Shape
hi = h ∪ i where
    h = (rectangle (0.1, 0.1) (0.25, 0.9) ∪
         rectangle (0.1, 0.1) (0.6, 0.35) ∪
         circle (0.35, 0.35) 0.25) ∩ inv
         (circle (0.35, 0.35) 0.1 ∪
          rectangle (0.25, 0.1) (0.45, 0.35))
    i = rectangle (0.75, 0.1) (0.9, 0.55) ∪
        circle (0.825, 0.75) 0.1

Rendering this to a bitmap is a matter of evaluating the function on many (x,y) coordinates, then coloring them depending on the function's sign. This is a very naïve rendering strategy, but it's good enough to put pixels on the screen.

Here's our testing shape, evaluated as a bitmap:

Hi!

Quadtrees

Quadtrees are a basic data structure in 2D graphics. Quadtreees are hierarchical: a 2D region is recursively split into four quadrants; each quadrant is either a leaf cell or subdivided further.

They are interesting because they represent a shape with fewer samples than a bitmap: areas that are completely filled or empty can be collapsed into a single cell.

We'll use one type of hierarchical cell and three terminal cells in the tree:

As motivation, here's a quadtree decomposition of our test shape. Green, grey, and black squares represent Leaf, Empty, and Full cells respectively.

Quadtree of 'hi'

Now that we know what we're looking for, let's start by creating our data type. Subtrees and corners will always be numbered as follows:

Quadtree numbering

For the purposes of this exercise, we'll make a generic Tree_ type that stores some kind of data at each leaf, then specialize it to storing the leaf's corner positions.

data Tree_ a = Root (Tree_ a) (Tree_ a) (Tree_ a) (Tree_ a) |
               Empty | Full | Leaf a
     deriving (Show, Eq)

type Cell = (Min, Max)
type Tree = Tree_ Cell

Now, let's build a complete tree structure:

buildTree :: Min -> Max -> Int -> Tree
buildTree min max 0 = Leaf (min, max)
buildTree (xmin, ymin) (xmax, ymax) i =
    Root (buildTree (xmin, ymin) (xmid, ymid) (i - 1))
         (buildTree (xmid, ymin) (xmax, ymid) (i - 1))
         (buildTree (xmin, ymid) (xmid, ymax) (i - 1))
         (buildTree (xmid, ymid) (xmax, ymax) (i - 1))
    where xmid = (xmin + xmax) / 2
          ymid = (ymin + ymax) / 2

This tree isn't specialized to any particular shape; it contains a uniform set of leaf cells filling a space. For example:

λ> buildTree (0,0) (1,1) 2

Full quadtree

It may look like we're building a full tree here, but that's not quite the case. Because Haskell is lazy, this tree will actually be constructed on an as-needed basis, so we need not feel guilty about building it naïvely.

Now, let's think about collapsing cells. Given a Shape function, we can mark cells that are entirely filled or empty, then recursively collapse them:

collapse :: Shape -> Tree -> Tree

collapse shape leaf@(Leaf ((xmin, ymin), (xmax, ymax))) =
    if      all (<  0) values then  Full
    else if all (>= 0) values then  Empty
    else                            leaf
    where values = [shape (x,y) | x <- [xmin, xmax],
                                  y <- [ymin, ymax]]

collapse shape (Root a b c d) =
    collapse' $ map (collapse shape) [a, b, c, d]
    where collapse' [Empty, Empty, Empty, Empty] = Empty
          collapse' [Full, Full, Full, Full] = Full
          collapse' [q, r, s, t] = Root q r s t

collapse _ t = t

Here's what the first five quadtree levels look like on our example shape:

Quadtree

Note that the contour-containing cells are all Leaf cells (and vice versa):

Marching squares, which will be introduced shortly, operates on a single cell at a time. To map it across all of the Leaf cells in the shape, let's make the Tree_ type an instance of the Foldable typeclass:

instance Foldable Tree_ where
    foldMap f (Leaf a) = f a
    foldMap f (Root a b c d) = mconcat $
                               map (foldMap f) [a, b, c, d]
    foldMap _ _ = mempty

This means we can easily map functions across every Leaf in the tree (or, more strictly, across every Cell stored in a leaf). For example, consider finding the centers of each leaf:

λ> let center ((x,y), (x',y')) = ((x+x')/2, (y+y')/2)
λ> foldMap (\cell -> [center cell]) $ buildTree (0,0) (1,1) 2
[(0.125,0.125),(0.375,0.125),(0.125,0.375),(0.375,0.375)
 (0.625,0.125),(0.875,0.125),(0.625,0.375),(0.875,0.375),
 (0.125,0.625),(0.375,0.625),(0.125,0.875),(0.375,0.875),
 (0.625,0.625),(0.875,0.625),(0.625,0.875),(0.875,0.875)]

Leaf positions

Marching Squares

Marching Squares is a fundamental algorithm for extracting isocontours from 2D samples.

Here's the idea: for every cell, we examine the corners, which match one of sixteen cases:

Marching squares table

Each case creates between zero and two edges. Run this procedure on every contour-containing cell and you'll get your shape's contour.

To get started, we'll encode the 16 cases shown above in a lookup table named lut.

data Side = Upper | Lower | Left | Right deriving Show

-- This lookup table takes a bitmask abcd and
-- returns a list of edges between which we
-- should draw contours (to outline the shape)
lut :: [[(Side, Side)]]
lut = [[],                          -- 0000
       [(Upper,Right)],             -- 000d
       [(Left,Upper)],              -- 00c0
       [(Left,Right)],              -- 00cd
       [(Right,Lower)],             -- 0b00
       [(Upper,Lower)],             -- 0b0d
       [(Right,Lower),(Left,Upper)],-- 0bc0
       [(Left,Lower)],              -- 0bcd
       [(Lower,Left)],              -- a000
       [(Lower,Left),(Upper,Right)],-- a00d
       [(Lower,Upper)],             -- a0c0
       [(Lower,Right)],             -- a0cd
       [(Right,Left)],              -- ab00
       [(Upper,Left)],              -- ab0d
       [(Right,Upper)],             -- abc0
       []]                          -- abcd

Next, we'll write helper functions that go from a Cell to an index in the lookup table and then to a list of Side pairs.

index :: Shape -> Cell -> Int
index shape ((xmin, ymin), (xmax, ymax)) =
    sum [if shape pt < 0 then 2^(3 - i) else 0 |
         (pt, i) <- zip pts [0..]]
    where pts = [(x,y) | y <- [ymin, ymax], x <- [xmin, xmax]]

edges :: Shape -> Cell -> [(Side, Side)]
edges shape c = lut !! index shape c

The edges function gives us pairs of cell sides between which we should draw edges; however, we need to pick the starting and ending point for each edge:

Finding points

In the image above, how do we find the green points? Well, the green points should be points where the shape function is zero. We'll assume that the shape function is monotonic along the edge and use binary search:

pt :: Shape -> Cell -> Side -> Point
pt shape ((xmin, ymin), (xmax, ymax)) side =
    case side of Left  -> zero shape (xmin, ymin) (xmin, ymax)
                 Right -> zero shape (xmax, ymin) (xmax, ymax)
                 Lower -> zero shape (xmin, ymin) (xmax, ymin)
                 Upper -> zero shape (xmin, ymax) (xmax, ymax)

zero :: Shape -> Point -> Point -> Point
zero s a@(ax, ay) b@(bx, by)
    | s a >= 0 = zero s b a
    | otherwise = zero' 0.5 0.25 10
    where pos f = (ax * (1-f) + bx * f, ay * (1-f) + by * f)
          zero' f step i =
            if i == 0 then pos f
            else if s (pos f) < 0 then
                    zero' (f + step) (step / 2) (i - 1)
            else    zero' (f - step) (step / 2) (i - 1)

Finally, we can combine these functions into one that converts a cell into a set of contours. Notice that this function can be applied to a tree with foldMap!

type Edge = (Point, Point)

contours :: Shape -> Cell -> [Edge]
contours shape cell = [(pt' a, pt' b) |
                       (a, b) <- edges shape cell]
    where pt' = pt shape cell

Applying this to a tree extracts the shape's contours:

λ> foldMap (contours hi) $ collapse hi $ buildTree (0,0) (1,1) 5

(factoring out hi with a monad is left as an exercise for the reader)

Marching squares

Sampled distance fields

There's an interesting nuance to our quadtrees: the collapse function treats sampled data as a binary value (based on whether it's less than zero). Can we do better by using the scalar value?

Consider the following quadtree (with one level of subdivision):

Sampled

The internal samples have different boolean values (some are inside the shape, some are outside), but all of the internal sampled values can be reconstructed from the corner values using bilinear interpolation.

This means that we can merge the cells without losing any information!
(This idea appears in the literature as Adaptively Sampled Distance Fields).

This may seem like an uncommon occurance, but it turns out that this kind of merging is possible whenever a cell's distance field is entirely determined by a single line (or face in 3D). When working with flat surfaces, this happens a lot!

Let's write some code to do this kind of reduction. First, we'll create functions that measure the error that would be introduced by this kind of merging:

-- interpolate samples a Shape at the four corners of a Cell,
-- then uses those values to estimate the function's result at
-- an arbitrary (x,y) position
interpolate :: Shape -> Cell -> Point -> Double
interpolate shape ((xmin, ymin), (xmax, ymax)) (x,y) =
    let dx = (x - xmin) / (xmax - xmin)
        dy = (y - ymin) / (ymax - ymin)
        ab = (shape (xmin, ymin)) * (1 - dx) +
             (shape (xmax, ymin)) * dx
        cd = (shape (xmin, ymax)) * (1 - dx) +
             (shape (xmax, ymax)) * dx
    in ab * (1 - dy) + cd * dy

-- score returns the difference between an interpolated
-- estimate and the function value at a given point
score :: Shape -> Cell -> Point -> Double
score shape cell pt = abs $
                      (interpolate shape cell pt) - shape pt

Now, we can use the score function and a user-defined cutoff to decide whether to merge a set of leaf cells. Note the names for nodes and cells:

Named

merge :: Shape -> Tree -> Tree
merge shape (Root a b c d) =
    merge' a' b' c' d'
    where [a', b', c', d'] = map (merge shape) [a, b, c, d]
          merge' (Leaf (min, i)) (Leaf (q, r))
                 (Leaf (s, t)) (Leaf (_, max)) =
            let scores = map (score shape (min, max))
                         [i, q, r, s, t]
            in if all (< 0.001) scores
               then Leaf (min, max)
               else Root a' b' c' d'
          merge' _ _ _ _ = Root a' b' c' d'
merge _ t = t

We apply this function before calling collapse to produce a sparser quadtree:

λ> collapse hi $ merge hi $ buildTree (0,0) (1,1) res

SDF

As expected, we now see larger leaf cells containing a single flat edge.

We can use the same contouring strategy on this sparser quadtree, producing a contour with fewer samples along flat edges:

Contoured SDF

This seems too easy — what's the catch?

Our meshing strategy assumed that all leaf cells were the same size, so a pair of leaves that shared an edge performed the same zero-crossing search. Now, with leaves of different sizes, the zero-crossing search will not be numerically identical, so the edges will not quite touch:

Cracks

This becomes a more significant issue in 3D: naïvely meshing a hierarchical structure will leave cracks in between cells of different sizes.

Dual contouring

Dual Contouring addresses two issues:

It does so by considering the dual graph of the quadtree. Instead of each leaf containing zero or more edges, each leaf stores a vertex; we connect these vertices to create the contour.

In this image, the original quadtree is drawn in light grey and the dual grid is shown in black.

Dual grid

There are three important questions to answer:

To answer the first question, we want to draw an edge between any two cells that share a side with a sign change.

Sign change

The answer to the second question is more subtle: we want to position vertices of the dual graph on features of the original shape. For example, in a cell containing a corner, the vertex should be on that corner. This strategy allows the algorithm to preserve sharp features.

We'll find features by looking at the normals of the shape function, which are given by its partial derivatives. In the past, I've written my own automatic differentiator; here, we'll do it numerically:

deriv :: Shape -> Point -> Point
deriv shape (x,y) =
    let epsilon = 0.001
        dx = shape (x + epsilon, y) - shape (x - epsilon, y)
        dy = shape (x, y + epsilon) - shape (x, y - epsilon)
        len = sqrt $ dx**2 + dy**2
    in (dx / len, dy / len)

A circle's partial derivatives are always pointing outwards from the center:

λ> deriv (circle (0,0) 1) (0,1)
(0.0,1.0)
λ> deriv (circle (0,0) 1) (1,1)
(0.7071067811865476,0.7071067811865476)
λ> deriv (circle (0,0) 1) (1,0)
(1.0,0.0)

Circle normals

To find a feature in a cell, we find the position and normal of every point along the cell's sides that intersects the isosurface, then perform a least-squares fit to estimate the feature position.

This is the same technique I use in Antimony, though we'll skip the sanity-checking for simplicity. The details are described in Feature Sensitive Surface Extraction from Volume Data.

The code below uses the HMatrix library for matrix math.

feature :: Shape -> Cell -> Maybe Point
feature shape cell =
    if length pts_ >= 2 then
        let pts = map fromTuple pts_
            nms = map (fromTuple . deriv shape) pts_
            center = sum pts / (fromIntegral $ length pts)

            a = fromRows pts
            b = col $ zipWith (\pt nm -> (pt - center) <·> nm)
                      pts nms

            p = center + (toColumns $ linearSolveSVD a b) !! 0
        in Just $ (\[x,y] -> (x,y)) $ (toList p)
    else Nothing
    where pts_ = concatMap (\(a,b) -> [a,b]) $ contours shape cell
          fromTuple = \(x,y) -> fromList [x,y]

We can run this function across every leaf cell to see where vertices are placed:

λ> let tree = collapse hi $ merge hi $ buildTree (0,0) (1,1) 5
λ> catMaybes $ foldMap (\cell -> [feature hi cell]) tree

DC point placement

This looks promising! There's one remaining step: to link these points together into a shape-outlining contour.

We use a set of three recursive functions to isolate edges of the dual graph. They're tricky, so I'll present them three times: in a diagram, in plain English, and in code.

First, the diagram:

Dual construction

faceProc is called on one cell. If that cell is a leaf, then it returns nothing. If that cell is a root, it recurses by calling faceProc on each subtree, edgeProcH on each horizontal pair of cells, and edgeProcV on each vertical pair of cells.

edgeProcH is called on a pair of cells. If both are leafs, it creates a contour edge between the two cells. Otherwise, it recursively calls itself on two horizontal pairs of cells.

edgeProcV is identical to edgeProcH, but applies to vertical pairs of cells.

Now, in code:

dc :: Shape -> Tree -> Edge
dc = faceProc

faceProc :: Shape -> Tree -> [Edge]
faceProc shape (Root a b c d) =
    concatMap (faceProc shape) [a,b,c,d] ++
    edgeProcH shape a b ++ edgeProcH shape c d ++
    edgeProcV shape a c ++ edgeProcV shape b d
faceProc _ _ = []

edgeProcH :: Shape -> Tree -> Tree -> [Edge]
edgeProcH shape (Leaf a) (Leaf b) =
    [(fromJust $ feature shape a, fromJust $ feature shape b)]
edgeProcH shape leaf@(Leaf _) (Root a _ c _) =
    edgeProcH shape leaf a ++ edgeProcH shape leaf c
edgeProcH shape (Root _ b _ d) leaf@(Leaf _) =
    edgeProcH shape b leaf ++ edgeProcH shape d leaf
edgeProcH shape (Root _ b _ d) (Root a _ c _) =
    edgeProcH shape b a ++ edgeProcH shape d c
edgeProcH _ _ _ = []

edgeProcV :: Shape -> Tree -> Tree -> [Edge]
edgeProcV shape (Leaf a) (Leaf b) = [(fromJust $ feature shape a,
                                      fromJust $ feature shape b)]
edgeProcV shape (Root _ _ c d) leaf@(Leaf _) =
    edgeProcV shape c leaf ++ edgeProcV shape d leaf
edgeProcV shape leaf@(Leaf _) (Root a b _ _) =
    edgeProcV shape leaf a ++ edgeProcV shape leaf b
edgeProcV shape (Root _ _ c d) (Root a b _ _) =
    edgeProcV shape c a ++ edgeProcV shape d b
edgeProcV _ _ _ = []

This completes our implementation of Dual Contouring:

λ> dc hi $ collapse hi $ merge hi $ buildTree (0,0) (1,1) res

Dual contouring in action

At this point, we'll declare victory: the contour is a faithful representation of the original shape and point density along the contour scales with local complexity / curvature.

There's plenty more to dive into in the literature, with references below.

Thanks to Richard Bowen and Michael Gilik for their feedback on an early draft.

Download the source here (including code used to generate diagrams).

References

Papers

Other sources

Providence Salumu