= Draw Me A Diagram = https://eprints.soton.ac.uk/257577/1/funcgeo2.pdf[Peter Henderson's _Functional Geometry_] composes a few disarmingly succinct functions to produce striking pictures. It is another satisfying example of simple building blocks unexpectedly leading to complex creations. See also https://dl.acm.org/doi/pdf/10.1145/800068.802148[the original 1982 version of the paper]. https://diagrams.github.io/[The `diagrams` Haskell package] builds on these ideas to provide a rich powerful language for vector graphics. I had used it to generate images for my webpages, 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?
[ ]:
-- Placeholder for examples.
== 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.
[ ]:
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 mirror (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 https://en.wikipedia.org/wiki/Atan2[`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 https://en.wikipedia.org/wiki/CORDIC[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
| theta < 0 = second negate $ cossin (-theta)
| theta <= pi/4 = cossinSmall theta
| theta <= pi/2 = (\(a,b) -> (b,a)) $ cossin (pi/2 - theta)
| theta <= pi = first negate $ cossin (pi - theta)
| theta <= 2*pi = (\(a,b) -> (-a, -b)) $ cossin $ theta - pi
| otherwise = cossin $ theta - 2*pi
cos = fst . cossin
sin = snd . cossin
cis = uncurry (:+) . cossin
Outside of elementary arithmetic, the only mathematical operation our code depends on is taking the square root, which WebAssembly conveniently provides. But if I had to code it myself, I'd likely refer to http://www.numberworld.org/y-cruncher/algorithms.html[the algorithms used by y-cruncher]. == Shapes == Our core data type is `Shape`:
[ ]:
data Shape = Shape
{ _envelope :: Com -> Double
, _trace :: Com -> Com -> [Double]
, _svg :: Double -> String -> String
, _named :: [(String, (Com -> Com, Shape))]
}
Each `Shape` is implicitly equipped with a point that we call its _local
origin_.
\(
\newcommand\O{\textbf{O}}
\newcommand\P{\textbf{P}}
\newcommand\Q{\textbf{Q}}
\newcommand\v{\textbf{v}}
\newcommand\w{\textbf{w}}
\)
Let \(D\) be a `Shape` with local origin \(\O\).
The _envelope_ of \(D\) is a function that takes a direction \(\v\) and
returns the scalar \(s\) given by:
\[
s = \sup \{ (\P - \O) \cdot \v | \P \in D \}
\]
In other words, the scalar \(s\) is the smallest value for which the plane
through \(\O + s \v\) normal to \(\v\) partitions space so that one half
contains the entirety of \(D\). Roughly speaking, if you were to walk in the
direction \(v\) starting from \(O\), then \(s\) tells you how far you must
travel so you no longer see \(D\), not even in your peripheral vision.
*Example*: for the unit circle whose local origin is its center, the envelope
is `recip . magnitude`.
*Example*: for the 1D unit circle, that is, the points -1 and 1, with local
origin 0, the envelope is the normalized projection along the real axis:
------------------------------------------------------------------------
\v@(x :+ _) -> abs x/norm v
------------------------------------------------------------------------
The _trace_ of \(D\) is a function that takes a point \(\P\) and a direction
\(\v\) and returns the set:
\[
\{ s | \P + s \v \in D \}
\]
(I haven't decided what to do about intervals within this set. Perhaps I could
replace them with their endpoints, or remove them entirely.)
In other words, the `trace` function identifies all the boundary points along a
given ray. In fact, this function is so named because of ray-tracing, and has
nothing to do with other trace functions in mathematics.
While the envelope function always returns results with respect to the local
origin \(\O\) of a diagram, the trace function must be given a starting point
\(\P\).
We represent the set of scalars with a sorted list.
The examples on this page only use the largest element, that is, the outermost
boundary point:
[ ]:
maxTraceV p pt dir = case _trace p pt dir of [] -> Nothing ss -> let s = last ss in if s <= 0 then Nothing else Just s maxTraceP p pt dir = ($ dir) . dilate <$> maxTraceV p pt dir
The `_svg` function returns a snippet of SVG that draws the shape as a difference list. It takes a scaling factor as a parameter so we can generate scale-invariant SVG for line widths, arrow heads, and so on. We hardcode the line-width to a value that works well for diagrams around the same size as a unit circle. SVG uses screen coordinates, which we make a little less confusing with `yshows` rather than mysteriously negate \(y\) coordinates here and there. However, we do simply negate the angle of rotation when needed. We define a helper that exports a `Shape` to SVG given a desired number of pixels per unit length. in the diagram with 1.1 units of padding. We call the envelope function to size the SVG appropriately.
[ ]:
lineWidth = 0.04 yshows = shows . negate svg pxPerUnit p = concat [ "<svg style='font-family:MJXZERO,MJXTEX-I;'" , " width=", show wPx , " height=", show hPx , " viewBox='", unwords (map show [x,y,w,h]), "'" , "><g font-size='0.8px'>", _svg p 1.0 "</g>" , "</svg>" ] where pad = 1.1 x0 = -(_envelope p -1) y0 = -(_envelope p i) w0 = _envelope p 1 - x0 h0 = _envelope p -i - y0 x = x0 - pad y = y0 - pad w = w0 + 2*pad h = h0 + 2*pad wPx = w * pxPerUnit hPx = h / w * wPx
We may name a `Shape` with a string so we can easily, say, connect two previously declared shapes. We implement this feature with the `_named` function. For a `Shape` \(D\), it returns a `Map` where each entry's key is the name of a component `Shape` of \(D\). The corresponding value is a tuple `(f, p)` where `p` is the component `Shape`, and `f` is a function that transforms coordinates with respect to the local origin of `p` to coordinates with respect to the local origin of \(D\). The following assigns a string name to a `Shape`:
[ ]:
named s p = let p' = p { _named = (s, (id, p')) : _named p } in p'
== 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: 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.
[ ]:
onNamed f p = second (first (f .)) <$> _named p
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 = onNamed (dilate n) 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 = onNamed (d+) 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 = onNamed (z*) p
}
where z = cis theta
reflectY :: Shape -> Shape
reflectY prim = Shape
{ _envelope = \dir -> _envelope prim $ mirror dir
, _trace = \pt dir -> _trace prim (mirror pt) $ mirror dir
, _svg = \zoom -> ("<g transform='scale(-1,1)'>"++) . _svg prim zoom . ("</g>"++)
, _named = onNamed mirror prim
}
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
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 <> _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.
[ ]:
innerPoints aName bName p = do (fa, a) <- lookup aName $ _named p (fb, b) <- lookup bName $ _named p let d = fb 0 - fa 0 pa <- fa <$> maxTraceP a (0:+0) d pb <- fb <$> maxTraceP b (0:+0) (negate d) pure (pa, pb) anglePoints aName bName aRad bRad p = do (fa, a) <- lookup aName $ _named p (fb, b) <- lookup bName $ _named p pa <- fa <$> maxTraceP a (0:+0) (cis aRad) pb <- fb <$> maxTraceP b (0:+0) (cis bRad) pure (pa, pb) straightArrowWith tip lineAttrs aName bName p = maybe p (p <>) do (pa, pb) <- innerPoints aName bName p 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 (pa, pb) <- anglePoints aName bName aRad bRad p 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 $ "repl.outdiv.insertAdjacentHTML('beforeend',`" ++ svg 20 p ++ "`);"
pure ()
== Teach a functional programmer to fish... == Let's retread the path taken by Henderson to draw a version of M. C. Escher's woodcut _Square Limit_. We start with a fish tile that we view as a fancy isoceles right triangle. The envelope and trace functions treat the shape as such a triangle, though we actually modify the border and add a center line and eyes. The redrawing of the border must obey constraints so that certain transformations of the tile fit perfectly with the original. We may freely modify one of the short edges of the triangle, but then the other short edge must be a duplicate, and the long edge must be two scaled-down flipped duplicates joined end to end after rotating one of them by 180 degrees.
[ ]:
fish = Shape
{ _envelope = \dir -> foldr1 max $ (\d -> dot d dir / dot dir dir) <$> tri
, _trace = \pt dir -> sort $ raySegment (pt, dir) =<< zip tri (tail tri ++ tri)
, _svg = \zoom -> ("<polygon fill='none' stroke='black' stroke-width='"++)
. shows (lineWidth*zoom)
. ("' points='"++)
. foldr (.) id (((' ':) .) . screenShow <$> ws)
. ("' />"++)
. ("<polyline points='-2,-2 -0.8,0.8 1,1.6' fill='none' stroke='black' stroke-width='"++) . shows (lineWidth*zoom) . ("' />"++)
. ("<polygon fill='none' stroke='black' stroke-width='"++) . shows (lineWidth*zoom) . ("' points='-1.8,-1 -1.85,-0.2 -1.55,-0.6' />"++)
. ("<polygon fill='none' stroke='black' stroke-width='"++) . shows (lineWidth*zoom) . ("' points='-1.55,-1.3 -1.35,-0.8 -1.1,-0.9' />"++)
, _named = mempty
}
where
tri = [-2:+ -2,2:+ -2,-2:+2]
ws = (- (2:+2)) <$> 0 : ws1 ++ ws2 ++ ws3 ++ ws4
ws1 = [1:+0.8,2,4]
ws2 = reverse $ (2:+2 +) . (1/sqrt 2 :+ 0 *) . (cis (-tau/8) *) <$> cs
ws3 = (2:+2 +) . (1/sqrt 2 :+ 0 *) . (cis (3*tau/8) *) <$> cs
ws4 = reverse $ (0:+1 *) <$> ws1
cs = conjugate <$> ws1
screenShow (x :+ y) = (shows x) . (' ':) . (yshows y)
draw $ scale 2 $ fish <> circle 4
An isoceles right triangle is half of a square, and the other half can be bisected into two isoceles right triangles. Thus we can adjoin two scaled-down copies of the fish tile to form a square that we call `t`. A square can also be cut into 4 congruent isoceles right triangles, which leads forming a tile `u` by joining 4 copies of our fish tile. We now see the reason for the above constraints. We treat `t` and `u` as squares. though they have wonky borders when drawn.
[ ]:
ro = rotateBy $ tau/4 unro = rotateBy $ -tau/4 fish2 = reflectY . translateY 2 . scale (1/sqrt 2) . rotateBy (tau/8) $ fish t = fish <> fish2 <> unro fish2 u = foldr1 (<>) $ take 4 $ iterate ro fish2 draw $ scale 2 t draw $ scale 2 u
We build a basic side tile with two blank tiles and two `t` tiles with one rotated so that they align. We iterate to glue on smaller copies of `t` tiles, resulting in more intricate side tiles. The corner tiles are similar, except we also need side tiles at the same level recursion.
[ ]:
sidesCorners = iterate (\(s, c) -> ( translateX -1 $ scale 0.5 $ (s ||| s) === (ro t ||| t) , scale 0.5 $ (c ||| s) === (ro s ||| u)) ) (blank, blank) where blank = scale 2 $ hstrut <> vstrut mapM (draw . scale 2 . fst) $ tail $ take 4 sidesCorners mapM (draw . scale 2 . snd) $ tail $ take 4 sidesCorners
To finish off, we surround a central tile `u` with side and corner tiles at the desired level of recursion.
[ ]:
sqlims = uncurry go <$> sidesCorners where
go s c = vcat
[ c ||| s ||| unro c
, ro s ||| u ||| unro s
, ro c ||| ro (ro s) ||| ro (ro c)
]
mapM (draw . scale 2) $ take 4 sqlims
We deviate from Henderson's paper slightly: * The `cycle` and `v` definitions are unneeded. Speaking of which, the paper seems to have the wrong definition of `cycle`: the 3rd and 4th arguments should be swapped. * Instead of `quartet` and `nonet`, we take advantage of our infix operators `(|||)` and `(===)`. * We compute side and corner tiles simultaneously so their recursion depths are automatically in sync. * Escher's work really has two kinds of tiles: the centermost fish have one squarish fin, but all the other ones, the fin stretches out and curves a little so that the two fins look more similar in size and shape, increasing the aesthetic appeal. Henderson's code gets away with a single tile because certain parts of the tile actually overlap when constructing the `u` tile. The lines in the overlapping regions superimpose one another so there is no clash, but it does cause an extra line to appear on one of the fins of all fish outside the center. We, on the other hand, just put up with lopsided fins.
[ ]:
draw $ (circle 3 === circle 1) ||| (circle 4 <> (circle 1 === circle 5))
[ ]:
draw $ hcat $ scale 1.5 . rotateBy (tau/5) . cyclogon <$> [3..10]
[ ]:
label = translate (-0.4 :+ -0.25) . text object s = label s <> (unitCircle |> named s) draw $ hcat [ object "Z" , translateY (-2) (label "h") <> translateY 2 (label "g") <> strutX 4 , object "X" , translateY 0.5 (label "f") <> strutX 4 , object "Y" ] |> straightArrow "X" "Y" |> curvedArrow (-tau/6) "Z" "X" (tau/8) (tau*3/8) |> curvedArrow (tau/6) "Z" "X" (-tau/8) (-tau*3/8)
[ ]:
jsEval [r|
const cells = convo.getElementsByClassName("cell");
const c2 = cursor.previousSibling;
const c1 = c2.previousSibling;
const c0 = c1.previousSibling;
cells[1].replaceWith(c0, c1, c2);
cursor.remove();
cursor = undefined;
|]