# Draw Me A Diagram

Peter Henderson’s paper, Functional Geometry (see also the original 1982 version), composes a few disarmingly succinct functions to produce striking pictures.

The diagrams Haskell package builds on these ideas, and I had been using it to generate SVG images, but its many dependencies can be unpleasant when switching GHC versions. As I only need a tiny subset of its features, why not roll my own version?

Compiling demo...

## Complex Numbers

We focus on 2D diagrams, representing points with complex numbers.

We define the dot product of two complex numbers to be:

$(a + bi) \cdot (c + di) = ac + bd$

The resulting real is the product of their magnitudes and the cosine of the angle between them.

import Map

infixl 8 |>
(|>) = flip ($) infix 6 :+ data Com = Double :+ Double deriving (Eq, Show) dot (a :+ b) (c :+ d) = a*c + b*d dilate r (a :+ b) = r*a :+ r*b realPart (a :+ _) = a conjugate (a :+ b) = a :+ -b norm x = dot x x magnitude = sqrt . norm i = 0 :+ 1 instance Ring Com where (a :+ b) + (c :+ d) = (a + c) :+ (b + d) (a :+ b) - (c :+ d) = (a - c) :+ (b - d) (a :+ b) * (c :+ d) = (a*c - b*d) :+ (a*d + b*c) fromInteger a = fromInteger a :+ 0 instance Field Com where recip x = dilate (recip$ norm x) $conjugate x We also roll our own trigonometric functions for our homebrew Haskell compiler. We approximate $$\arctan(x)$$ for $$x\in [0..\frac{1}{2}]$$ with its Taylor series expansion: $\arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} …​$ and hardcode $$\arctan(1) = \pi/4$$. We ensure the smallest items are added first with foldr to fight floating-point rounding error. pi = 3.141592654 tau = 2*pi atanTaylor 1 = pi * 0.25 atanTaylor x = foldr (+) 0$ reverse $take 25$
zipWith (/) (iterate (*(-x*x)) x) (iterate (+2) 1)

We build our phase function on top, which is also known as atan2.

phase (x :+ y)
| x < 0 = if y < 0
then phase' -x -y - pi
else pi - phase' -x y
| y < 0 = -(phase' x -y)
| otherwise = phase' x y
phase' x y
| y > x = 0.5*pi - phase'' x y
| otherwise = phase'' y x
phase'' y x
| x == 0 = if y == 0 then 0 else pi * 0.5
| 2*y > x = atanTaylor 0.5 + atanTaylor ((r - 0.5) / (1 + r*0.5))
| otherwise = atanTaylor r
where r = y / x

We use CORDIC for computing sines and cosines of a given small angle $$\theta$$.

This algorithm reminds me of binary search. In brief, we start with:

\begin{aligned} \alpha &= 0 \\ \sin\alpha &= 0 \\ \cos\alpha\ &= 1 \end{aligned}

then add to or subtract from $$\alpha$$ successively the values:

$\arctan(2^0), \arctan(2^{-1}), \arctan(2^{-2}), …​$

until $$\alpha\approx\theta$$. All the while, we update $$\sin\alpha$$ with corresponding additions or subtractions of $$2^{-k} \cos\alpha$$, with a similar update for $$\cos\alpha$$.

There is a wrinkle. At each step, a scaling factor creeps in, which we could normalize away immediately. However, it’s better to defer this so that we need only one final multiplication by a precomputed constant.

More concretely, the algorithm follows the identities:

\begin{aligned} \sqrt{1+2^{-2k}} \sin(\alpha \pm \arctan(2^{-k})) &= \sin\alpha \pm 2^{-k}\cos\alpha \\ \sqrt{1+2^{-2k}} \cos(\alpha \pm \arctan(2^{-k})) &= \cos\alpha \mp 2^{-k}\sin\alpha \end{aligned}

for $$k\in[0..n]$$, where $$n$$ is number of steps. The final normalization multiplies by:

$\prod_{k=0}^n \frac{1}{\sqrt{1+2^{-2k}}}$

We work with $$\beta = \theta - \alpha$$ instead of $$\alpha$$ directly so we can compare against zero.

cossinSmall theta = cordic lim tab theta 1 0 where
lim = 25
cordic n ((p2, a):rest) beta x y
| n == 0 = (kn * x, kn * y)
| otherwise = cordic (n - 1) rest (beta - sig a) x' y'
where
sig = if beta < 0 then negate else id
x' = x - sig p2 * y
y' = sig p2 * x + y
kn = kvalues!!lim
kvalues = scanl1 (*) $(\k -> 1/sqrt(1+0.5^(2*k))) <$> [0..]
tab = zip pows $atanTaylor <$> pows where pows = iterate (*0.5) 1

cossin theta

## Unit Circle

The unit circle is a good introductory example. We define its local origin to be the center of the circle. We compute its trace by solving a quadratic to find the points of intersection between a line and a unit circle.

unitCircle = Shape
{ _envelope = (1/) . magnitude
, _trace = ptCirc
, _svg = \zoom -> ("<circle fill='none' stroke='black' stroke-width='"++)
. shows (lineWidth*zoom) . ("' r=1 />"++)
, _named = mempty
}
where
ptCirc v dv
| disc < 0 = []
| otherwise = [(-b - sd) * aInv, (-b + sd) * aInv]
where
a = dot dv dv
b = dot v dv
c = dot v v - 1
disc = b^2 - a*c
aInv = 1 / a
sd = sqrt disc

## Regular Polygons

The $$n$$th roots of unity lie on the unit circle, and we can join them with edges to form a regular $$n$$-gon.

We compute its envelope by finding the maximum normalized projection of each vertex on to the given direction. For large $$n$$ it would be faster to test only the endpoints of the edge facing the given direction.

We are similarly wasteful when computing the trace. We compute ray-segment intersections for every edge, and sort any results.

To find the intersection of two lines, we solve equations of the following form for $$\lambda$$ and $$\mu$$:

$\P + \lambda \v = \Q + \mu \w$

As $$i \v \cdot \v = 0$$, we eliminate the $$\v$$ term by dotting both sides with $$i\v$$ to find:

$\mu = \frac{i \v \cdot (\P - \Q)}{i \v \cdot \w}$

Similarly, dotting with $$i\w$$ yields:

$\lambda = \frac{i\w \cdot (\Q - \P)}{i\w \cdot \v}$

These solutions fail when $$i\v\cdot\w = 0$$, that is, when the lines are parallel.

(Our code liberally uses the identity $$i\v\cdot\w = -i\w\cdot\v$$.)

sort [] = []
sort (x:xt) = sort (filter (<= x) xt) ++ [x] ++ sort (filter (> x) xt)

cyclogon n = Shape
{ _envelope = \dir -> foldr1 max $(\d -> dot d dir / dot dir dir) <$> vs
, _trace = \pt dir -> sort $raySegment (pt, dir) =<< zip vs (tail vs ++ vs) , _svg = \zoom -> ("<polygon fill='none' stroke='black' stroke-width='"++) . shows (lineWidth*zoom) . ("' points='"++) . foldr (.) id (((' ':) .) . screenShow <$> vs)
. ("' />"++)
, _named = mempty
}
where
vs = take n $iterate (cis(tau/fromIntegral n) *) 1 screenShow (x :+ y) = (shows x) . (' ':) . (yshows y) raySegment (p, v) (w1, w2) | d == 0 || b < 0 || b > 1 = [] | otherwise = [a] where d = dot (i*w) v x = w1 - p w = w1 - w2 a = dot (i*w) x / d b = dot (i*x) v / d ## Struts We define a horizontal strut to be an invisible 1D circle with no trace. A vertical strut is the analogous shape on the imaginary axis. hstrut = Shape { _envelope = \d@(dx :+ dy) -> abs dx/norm d , _trace = \_ _ -> [] , _svg = \zoom -> id , _named = mempty } vstrut = Shape { _envelope = \d@(dx :+ dy) -> abs dy/norm d , _trace = \_ _ -> [] , _svg = \zoom -> id , _named = mempty } ## Transforming Shapes We can easily handle some well-known transformations. Scaling: simply scale the envelope and trace by the same factor. Translation: for the trace, we undo the translation on $$\P$$ before computing the original trace; for the envelope, we compute the original envelope, then subtract the normalized projection of the translation vector on the given direction. Rotation: for both the trace and envelope, undo the rotation on the given direction before computing the original function. SVG has primitives for all these transformations. scale :: Double -> Shape -> Shape scale n prim = Shape { _envelope = \dir -> n * _envelope prim dir , _trace = \pt dir -> (n *) <$> _trace prim pt dir
, _svg = \zoom -> ("<g transform='scale("++) . shows n . (")'>"++) . _svg prim (zoom / n) . ("</g>"++)
, _named = first (((n:+0)*) .) <$> _named prim } translate :: Com -> Shape -> Shape translate d@(dx :+ dy) prim = Shape { _envelope = \dir -> _envelope prim dir + dot d dir / dot dir dir , _trace = \pt dir -> _trace prim (pt - d) dir , _svg = \zoom -> ("<g transform='translate("++) . shows dx . (' ':) . yshows dy . (")'>"++) . _svg prim zoom . ("</g>"++) , _named = first ((d+) .) <$> _named prim
}

translateX x = translate $x :+ 0 translateY y = translate$ 0 :+ y

rotateBy :: Double -> Shape -> Shape
rotateBy theta p = Shape
{ _envelope = \dir -> _envelope p (dir * conjugate z)
, _trace = \pt dir -> _trace p pt (dir * conjugate z)
, _svg = \zoom -> ("<g transform='rotate("++) . shows (-theta / pi * 180) . (")'>"++) . _svg p zoom . ("</g>"++)
, _named = first ((z*) .) <$> _named p } where z = cis theta We use a transformation to provide a handy function that returns a circle of any given radius. Hard-coding a dedicated Shape might perform better, but there’s no need to optimize yet. We define strutX and strutY similarly. It might be more consistent to have circle take a diameter parameter rather than a radius, but this breaks tradition. circle n = scale n$ unitCircle
strutX x = scale (x/2) hstrut
strutY y = scale (y/2) vstrut

We could generalize the scaling and rotation cases. If $$T$$ is an invertible linear transformation for a shape $$D$$, then to compute envelope of $$T D$$ on a vector $$\v$$ we compute the envelope of $$D$$ on $$T^{-1} \v$$, and similarly for the trace. (The scaling case then simplifies considerably due to linearity.)

Some care would be needed with SVG generation since we desire things like line width to be scale-invariant. Dividing the scaling parameter by the determinant of the matrix representing $$T$$ ought to do the trick.

## Composing Shapes

The atop function places one diagram atop another by lining up their local origins. The envelope of the combined diagrams is the maximum of their envelopes, while its trace is the union of their traces. As we represent sets with sorted lists, we combine the traces with merge sort.

This associative operation is a good choice for turning Shape into a semigroup.

mergeSort xs ys = case xs of
[] -> ys
x:xt -> case ys of
[] -> xs
y:yt | x <= y -> x:mergeSort xt ys
| True -> y:mergeSort xs yt

munion m1 m2 = foldr (uncurry insert) m1 $toAscList m2 atop :: Shape -> Shape -> Shape atop p q = Shape { _envelope = \dir -> _envelope p dir max _envelope q dir , _trace = \pt dir -> _trace p pt dir mergeSort _trace q pt dir , _svg = \zoom -> _svg p zoom . _svg q zoom , _named = _named p munion _named q } instance Semigroup Shape where (<>) = atop The pieces are in place for beside, which places one Shape next to another in a given direction so that their envelopes touch. We specialize a couple of directions so we can succinctly describe horizontal and vertical layouts. beside :: Com -> Shape -> Shape -> Shape beside dir x y = x <> translate (dilate (_envelope x dir + _envelope y (-dir)) dir) y (|||) = beside 1 (===) = beside -i hcat = foldr1 (|||) vcat = foldr1 (===) For shapes like arrows, arrowheads, and labels, we have no need for the envelope and trace. We introduce the ghost function to help define Shape values that are thin wrappers around various SVG drawings. ghost f = Shape { _envelope = const 0 , _trace = \_ _ -> [] , _svg = f , _named = mempty } text :: String -> Shape text s = ghost \zoom -> ("<text fill='black'>"++) . (s++) . ("</text>"++) svgFilledPolygon pts = ("<polygon fill='black' points='"++) . foldr (.) id (map (\(x :+ y) -> (" "++) . shows x . (" "++) . yshows y) pts) . ("'/>"++) dart = ghost \zoom -> ("<g transform='scale("++) . shows (6*lineWidth/zoom) . (")'>"++) . svgFilledPolygon [0, t1, t2, conjugate t1] . ("</g>"++) where t1 = cis (2/5 * tau) - (1 :+ 0) t2 = (realPart t1 + 1/2):+0 dubDart = ghost \zoom -> ("<g transform='scale("++) . shows (6*lineWidth/zoom) . (")'>"++) . svgFilledPolygon [0, t1, t2, conjugate t1] . svgFilledPolygon [t2, t3, t4, conjugate t3] . ("</g>"++) where t1 = cis (2/5 * tau) - (1 :+ 0) t2 = (realPart t1 + 1/2):+0 t3 = t1 + t2 t4 = t2 + t2 lineWith :: String -> Com -> Com -> Shape lineWith attrs (x1 :+ y1) (x2 :+ y2) = ghost \zoom -> ("<line stroke-width='"++) . shows (lineWidth*zoom) . ("' x1="++) . shows x1 . (" y1="++) . yshows y1 . (" x2="++) . shows x2 . (" y2="++) . yshows y2 . (" stroke='black' "++) . (attrs++) . (" />"++) dashedAttrs = "stroke-dasharray=0.1" -- Assumes rad lies in [-tau..tau]. arcline :: Com -> Com -> Double -> Shape arcline a@(x1 :+ y1) b@(x2 :+ y2) rad = ghost \zoom -> ("<path stroke-width='"++) . shows (lineWidth*zoom) . ("' fill='none' stroke='black' d='M "++) . shows x1 . (" "++) . yshows y1 . (" A "++) . shows r . (" "++) . shows r . (" 0 "++) . shows (fromEnum$ abs rad >= pi)  -- Large arc flag.
. (' ':)
. shows (fromEnum $rad < 0) -- Sweep flag . (' ':) . shows x2 . (" "++) . yshows y2 . ("' />"++) where r = magnitude (b - a) / (2 * abs (sin (rad / 2))) Next are utilities for drawing arrows between named diagrams. Here, we see the importance of changing coordinate systems: by the time we wish to draw arrows, the underlying objects may have undergone several transformations, so it would make no sense to use the original local coordinates of each endpoint. straightArrowWith tip lineAttrs aName bName p = maybe p (p <>)$ do
(fa, a) <- mlookup aName $_named p (fb, b) <- mlookup bName$ _named p
let d = fb 0 - fa 0
pa <- fa <$> maxTraceP a (0:+0) d pb <- fb <$> maxTraceP b (0:+0) (negate d)
let hd = translate pb $rotateBy (phase$ pb - pa) tip
pure $lineWith lineAttrs pa pb <> hd straightArrow = straightArrowWith dart "" existsArrow = straightArrowWith dart dashedAttrs curvedArrow rad aName bName aRad bRad p = maybe p (p <>)$ do
(fa, a) <- mlookup aName $_named p (fb, b) <- mlookup bName$ _named p
pa <- fa <$> maxTraceP a (0:+0) (cis aRad) pb <- fb <$> maxTraceP b (0:+0) (cis bRad)
let hd = translate pb $rotateBy (rad / 2)$ rotateBy (phase $pb - pa) dart pure$ arcline pa pb rad <> hd

Lastly, we have a wrapper that inserts an SVG into this webpage. The demo variable refers to a <div> element at the top of this page.

Each unit takes 20 pixels, which works well with demos with small numbers.

draw p = do
jsEval \$ "demo.insertAdjacentHTML('beforeend'," ++ svg 20 p ++ ");"
pure ()

Ben Lynn blynn@cs.stanford.edu 💡