cubicbezier-0.2.0: Geom2D/CubicBezier/MetaPath.hs
{-# LANGUAGE BangPatterns #-}
-- | This module implements an extension to paths as used in
-- D.E.Knuth's /Metafont/. Metafont gives a more intuitive method to
-- specify paths than bezier curves. I'll give a brief overview of
-- the metafont curves. For a more in depth explanation look at
-- /The MetafontBook/.
--
-- Each spline has a tension parameter, which is a relative measure of
-- the length of the curve. You can specify the tension for the left
-- side and the right side of the spline separately. By default
-- metafont gives a tension of 1, which gives a good looking curve.
-- Tensions shouldn't be less than 3/4, but this implementation
-- doesn't check for it. If you want to avoid points of inflection on
-- the spline, you can use @TensionAtLeast@ instead of @Tension@,
-- which will adjust the length of the control points so they fall
-- into the /bounding triangle/, if such a triangle exist.
--
-- You can either give directions for each node, or let metafont find
-- them. Metafont will solve a set of equations to find the
-- directions. You can also let metafont find directions at corner
-- points by setting the /curl/, which is how much the point /curls/
-- at that point. At endpoints a curl of 1 is implied when it is not
-- given.
--
-- Metafont will then find the control points from the path for you.
-- You can also specify the control points explicitly.
--
-- Here is an example path from the metafont program text:
--
-- @
-- z0..z1..tension atleast 1..{curl 2}z2..z3{-1,-2}..tension 3 and 4..z4..controls z45 and z54..z5
-- @
--
-- This path is equivalent to:
--
-- @
-- z0{curl 1}..tension atleast 1 and atleast 1..{curl 2}z2{curl 2}..tension 1 and 1..
-- {-1,-2}z3{-1,-2}..tension 3 and 4..z4..controls z45 and z54..z5
-- @
--
-- This path can be used with the following datatype:
--
-- @
-- OpenMetaPath [ (z0, MetaJoin Open (Tension 1) (Tension 1) Open)
-- , (z1, MetaJoin Open (TensionAtLeast 1) (TensionAtLeast 1) (Curl 2))
-- , (z2, MetaJoin Open (Tension 1) (Tension 1) Open)
-- , (z3, MetaJoin (Direction (Point (-1) (-2))) (Tension 3) (Tension 4) Open)
-- , (z4, Controls z45 z54)
-- ] z5
-- @
--
-- Cyclic paths are similar, but use the @CyclicMetaPath@ contructor.
-- There is no ending point, since the ending point will be the same
-- as the first point.
module Geom2D.CubicBezier.MetaPath
(unmeta, MetaPath (..), MetaJoin (..), MetaNodeType (..), Tension (..))
where
import Geom2D
import Geom2D.CubicBezier.Basic
import Data.List
import Text.Printf
data MetaPath = OpenMetaPath [(Point, MetaJoin)] Point
| CyclicMetaPath [(Point, MetaJoin)]
data MetaJoin = MetaJoin { metaTypeL :: MetaNodeType
, tensionL :: Tension
, tensionR :: Tension
, metaTypeR :: MetaNodeType
}
| Controls Point Point
deriving Show
data MetaNodeType = Open
| Curl {curlgamma :: Double}
| Direction {nodedir :: Point}
deriving Show
data Tension = Tension {tensionValue :: Double}
| TensionAtLeast {tensionValue :: Double}
deriving (Eq, Show)
instance Show MetaPath where
show (CyclicMetaPath nodes) =
showPath nodes ++ "cycle"
show (OpenMetaPath nodes lastpoint) =
showPath nodes ++ showPoint lastpoint
showPath :: [(Point, MetaJoin)] -> [Char]
showPath = concatMap showNodes
where
showNodes (p, Controls u v) =
showPoint p ++ "..controls " ++ showPoint u ++ "and " ++ showPoint v ++ ".."
showNodes (p, MetaJoin m1 t1 t2 m2) =
showPoint p ++ typename m1 ++ ".." ++ tensions ++ typename m2
where
tensions
| t1 == t2 && t1 == Tension 1 = ""
| t1 == t2 = printf "tension %s.." (showTension t1)
| otherwise = printf "tension %s and %s.."
(showTension t1) (showTension t2)
showTension (TensionAtLeast t) = printf "atleast %.3f" t :: String
showTension (Tension t) = printf "%.3f" t :: String
typename Open = ""
typename (Curl g) = printf "{curl %.3f}" g :: String
typename (Direction dir) = printf "{%s}" (showPoint dir) :: String
showPoint :: Point -> String
showPoint (Point x y) = printf "(%.3f, %.3f)" x y
-- | Create a normal path from a metapath.
unmeta :: MetaPath -> Path
unmeta (OpenMetaPath nodes endpoint) =
unmetaOpen (sanitizeOpen nodes) endpoint
unmeta (CyclicMetaPath nodes) =
case span (bothOpen . snd) nodes of
(l, []) -> unmetaCyclic l
(l, (m:n)) ->
if leftOpen $ snd m
then unmetaAsOpen (l++[m]) n
else unmetaAsOpen l (m:n)
unmetaOpen :: [(Point, MetaJoin)] -> Point -> Path
unmetaOpen nodes endpoint =
let subsegs = openSubSegments nodes endpoint
path = joinSegments $ map unmetaSubSegment subsegs
in OpenPath path endpoint
-- decompose into a list of subsegments that need to be solved.
openSubSegments :: [(Point, MetaJoin)] -> Point -> [MetaPath]
openSubSegments l p = openSubSegments' (tails l) p
openSubSegments' :: [[(Point, MetaJoin)]] -> Point -> [MetaPath]
openSubSegments' [[]] _ = []
openSubSegments' [] _ = []
openSubSegments' l lastPoint = case break breakPoint l of
(m, n:o) ->
let point = case o of
(((p,_):_):_) -> p
_ -> lastPoint
in OpenMetaPath (map head (m ++ [n])) point :
openSubSegments' o lastPoint
_ -> error "openSubSegments': unexpected end of segments"
-- join subsegments into one segment
joinSegments :: [Path] -> [(Point, PathJoin)]
joinSegments = concatMap nodes
where nodes (OpenPath n _) = n
nodes (ClosedPath n) = n
-- solve a cyclic metapath where all angles depend on the each other.
unmetaCyclic :: [(Point, MetaJoin)] -> Path
unmetaCyclic nodes =
let points = map fst nodes
chords = zipWith (^-^) points (last points : points)
tensionsA = (map (tensionL . snd) nodes)
tensionsB = (map (tensionR . snd) nodes)
turnAngles = zipWith turnAngle chords (tail $ cycle chords)
thetas = solveCyclicTriD $
eqsCycle tensionsA
points
tensionsB
turnAngles
phis = zipWith (\x y -> -(x+y)) turnAngles (tail thetas ++ [head thetas])
in ClosedPath $ zip points $
zipWith6 unmetaJoin points (tail points ++ [head points])
thetas phis tensionsA tensionsB
-- solve a cyclic metapath as an open path if possible.
-- rotate to the defined node, and rotate back after
-- solving the path.
unmetaAsOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)] -> Path
unmetaAsOpen l m = ClosedPath (l'++m')
where n = length m
OpenPath o _ = unmetaOpen (sanitizeCycle (m++l)) (fst $ head m)
(m',l') = splitAt n o
-- solve a subsegment
unmetaSubSegment :: MetaPath -> Path
-- the simple case where the control points are already given.
unmetaSubSegment (OpenMetaPath [(p, Controls u v)] q) =
OpenPath [(p, JoinCurve u v)] q
-- otherwise solve the angles, and find the control points
unmetaSubSegment (OpenMetaPath nodes lastpoint) =
let points = map fst nodes ++ [lastpoint]
joins = map snd nodes
chords = zipWith (^-^) (tail points) points
tensionsA = map tensionL joins
tensionsB = map tensionR joins
turnAngles = zipWith turnAngle chords (tail chords) ++ [0]
thetas = solveTriDiagonal $
eqsOpen points joins chords turnAngles
(map tensionValue tensionsA)
(map tensionValue tensionsB)
phis = zipWith (\x y -> -x-y) turnAngles (tail thetas)
pathjoins = zipWith6 unmetaJoin points (tail points) thetas phis tensionsA tensionsB
in OpenPath (zip points pathjoins) lastpoint
unmetaSubSegment _ = error "unmetaSubSegment: subsegment should not be cyclic"
bothOpen :: MetaJoin -> Bool
bothOpen (MetaJoin Open _ _ Open) = True
bothOpen _ = False
leftOpen :: MetaJoin -> Bool
leftOpen (MetaJoin Open _ _ _) = True
leftOpen _ = False
replaceLast :: [a] -> a -> [a]
replaceLast [] _ = []
replaceLast [_] n = [n]
replaceLast (l:ls) n = l:replaceLast ls n
sanitizeCycle :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
sanitizeCycle l = replaceLast ls l'
where
(l':ls) = sanitizeRest (last l: l)
-- replace open nodetypes with more defined nodetypes if possible
sanitizeOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
sanitizeOpen [] = []
-- starting open => curl
sanitizeOpen ((p, MetaJoin Open t1 t2 m):rest) =
sanitizeRest ((p, MetaJoin (Curl 1) t1 t2 m):rest)
sanitizeOpen l = sanitizeRest l
sanitizeRest :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
sanitizeRest [] = []
-- ending open => curl
sanitizeRest [(p, MetaJoin m t1 t2 Open)] =
[(p, MetaJoin m t1 t2 (Curl 1))]
sanitizeRest (node1@(p, MetaJoin m1 tl tr m2): node2@(q, MetaJoin n1 sl sr n2): rest) =
case (m2, n1) of
(Curl g, Open) -> -- curl, open => curl, curl
node1 : sanitizeRest ((q, MetaJoin (Curl g) sl sr n2):rest)
(Open, Curl g) -> -- open, curl => curl, curl
(p, MetaJoin m1 tl tr (Curl g)) : sanitizeRest (node2:rest)
(Direction dir, Open) -> -- given, open => given, given
node1 : sanitizeRest ((q, (MetaJoin (Direction dir) sl sr n2)) : rest)
(Open, Direction dir) -> -- open, given => given, given
(p, MetaJoin m1 tl tr (Direction dir)) : sanitizeRest (node2:rest)
_ -> node1 : sanitizeRest (node2:rest)
sanitizeRest ((p, m): (q, n): rest) =
case (m, n) of
(Controls _u v, MetaJoin Open t1 t2 mt2) -> -- explicit, open => explicit, given
(p, m) : sanitizeRest ((q, MetaJoin (Direction (q^-^v)) t1 t2 mt2): rest)
(MetaJoin mt1 tl tr Open, Controls u _v) -> -- open, explicit => given, explicit
(p, MetaJoin mt1 tl tr (Direction (u^-^p))) : sanitizeRest ((q, n): rest)
_ -> (p, m) : sanitizeRest ((q, n) : rest)
sanitizeRest (n:l) = n:sanitizeRest l
-- break the subsegment if the angle to the left or the right is defined or a curl.
breakPoint :: [(Point, MetaJoin)] -> Bool
breakPoint ((_, MetaJoin _ _ _ Open):(_, MetaJoin Open _ _ _):_) = False
breakPoint _ = True
-- solve the tridiagonal system for t[i]:
-- a[n] t[i-1] + b[i] t[i] + c[b] t[i+1] = d[i]
-- where a[0] = c[n] = 0
-- by first rewriting it into
-- the system t[i] + u[i] t[i+1] = v[i]
-- where u[n] = 0
-- then solving for t[n]
-- see metafont the program: ¶ 283
solveTriDiagonal :: [(Double, Double, Double, Double)] -> [Double]
solveTriDiagonal [] = error "solveTriDiagonal: not enough equations"
solveTriDiagonal ((_, b0, c0, d0): rows) = solutions
where
((_, vn): twovars) =
reverse $ scanl nextrow (c0/b0, d0/b0) rows
nextrow (u, v) (ai, bi, ci, di) =
(ci/(bi - u*ai), (di - v*ai)/(bi - u*ai))
solutions = reverse $ scanl nextsol vn twovars
nextsol ti (u, v) = v - u*ti
-- test = ((80.0,58.0,51.0),[(-432.0,78.0,102.0,503.0),(71.0,-82.0,20.0,2130.0),(52.39,-10.43,4.0,56.0),(34.0,38.0,0.0,257.0)])
-- solve the cyclic tridiagonal system.
-- see metafont the program: ¶ 286
solveCyclicTriD :: [(Double, Double, Double, Double)] -> [Double]
solveCyclicTriD rows = solutions
where
(!un, !vn, !wn): threevars =
reverse $ tail $ scanl nextrow (0, 0, 1) rows
nextrow (!u, !v, !w) (!ai, !bi, !ci, !di) =
(ci/(bi - ai*u), (di - ai*v)/(bi - ai*u), -ai*w/(bi - ai*u))
(totvn, totwn) = foldl (\(v', w') (u, v, w) ->
(v - u*v', w - u*w'))
(0, 1) threevars
t0 = (vn - un*totvn) / (1 - (wn - un*totwn))
solutions = scanl nextsol t0
((un, vn, wn) : reverse (tail threevars))
nextsol t (!u, !v, !w) = (v + w*t0 - t)/u
turnAngle :: Point -> Point -> Double
turnAngle (Point x y) q = vectorAngle $ rotateVec p $* q
where p = Point x (-y)
zipPrev :: [a] -> [(a, a)]
zipPrev [] = []
zipPrev l = zip (last l : l) l
-- find the equations for a cycle containing only open points
eqsCycle :: [Tension] -> [Point] -> [Tension]
-> [Double] -> [(Double, Double, Double, Double)]
eqsCycle tensionsA points tensionsB turnAngles =
zipWith4 eqTension
(zipPrev (map tensionValue tensionsA))
(zipPrev dists)
(zipPrev turnAngles)
(zipPrev (map tensionValue tensionsB))
where
dists = zipWith vectorDistance points (tail $ cycle points)
-- find the equations for an path with open points.
-- The first and last node should be a curl or a given angle
eqsOpen :: [Point] -> [MetaJoin] -> [Point] -> [Double]
-> [Double] -> [Double] -> [(Double, Double, Double, Double)]
eqsOpen _ [join] [delta] _ _ _ =
case join of
MetaJoin (Curl _) _ _ (Curl _) ->
[(0, 1, 0, 0), (0, 1, 0, 0)]
MetaJoin (Curl g) t1 t2 (Direction dir) ->
[eqCurl0 g (tensionValue t1) (tensionValue t2) 0,
(0, 1, 0, turnAngle delta dir)]
MetaJoin (Direction dir) t1 t2 (Curl g) ->
[(0, 1, 0, turnAngle delta dir),
eqCurlN g (tensionValue t1) (tensionValue t2)]
MetaJoin (Direction dir) _ _ (Direction dir2) ->
[(0, 1, 0, turnAngle delta dir),
(0, 1, 0, turnAngle delta dir2)]
_ -> error "eqsOpen: illegal nodetype in subsegment"
eqsOpen points joins chords turnAngles tensionsA tensionsB =
eq0 : restEquations joins tensionsA dists turnAngles tensionsB
where
dists = zipWith vectorDistance points (tail points)
eq0 = case head joins of
(MetaJoin (Curl g) _ _ _) -> eqCurl0 g (head tensionsA) (head tensionsB) (head turnAngles)
(MetaJoin (Direction dir) _ _ _) -> (0, 1, 0, turnAngle (head chords) dir)
_ -> error "eqsOpen: illegal subsegment first nodetype"
restEquations [lastnode] (tensionA:_) _ _ (tensionB:_) =
case lastnode of
MetaJoin _ _ _ (Curl g) -> [eqCurlN g tensionA tensionB]
MetaJoin _ _ _ (Direction dir) -> [(0, 1, 0, turnAngle (last chords) dir)]
_ -> error "eqsOpen: illegal subsegment last nodetype"
restEquations (_:othernodes) (tensionA:restTA) (d:restD) (turn:restTurn) (tensionB:restTB) =
eqTension (tensionA, head restTA) (d, head restD) (turn, head restTurn) (tensionB, head restTB) :
restEquations othernodes restTA restD restTurn restTB
restEquations _ _ _ _ _ = error "eqsOpen: illegal rest equations"
-- the equation for an open node
eqTension :: (Double, Double) -> (Double, Double)
-> (Double, Double) -> (Double, Double)
-> (Double, Double, Double, Double)
eqTension (tensionA', tensionA) (dist', dist) (psi', psi) (tensionB', tensionB) =
(a, b+c, d, -b*psi' - d*psi)
where
a = (tensionB' * tensionB' / (tensionA' * dist'))
b = (3 - 1/tensionA') * tensionB' * tensionB' / dist'
c = (3 - 1/tensionB) * tensionA * tensionA / dist
d = tensionA * tensionA / (tensionB * dist)
-- the equation for a starting curl
eqCurl0 :: Double -> Double -> Double -> Double -> (Double, Double, Double, Double)
eqCurl0 gamma tensionA tensionB psi = (0, c, d, r)
where
c = chi/tensionA + 3 - 1/tensionB
d = (3 - 1/tensionA)*chi + 1/tensionB
chi = gamma*tensionB*tensionB / (tensionA*tensionA)
r = -d*psi
-- the equation for an ending curl
eqCurlN :: Double -> Double -> Double -> (Double, Double, Double, Double)
eqCurlN gamma tensionA tensionB = (a, b, 0, 0)
where
a = (3 - 1/tensionB)*chi + 1/tensionA
b = chi/tensionB + 3 - 1/tensionA
chi = gamma*tensionA*tensionA / (tensionB*tensionB)
-- magic formula for getting the control points by John Hobby
unmetaJoin :: Point -> Point -> Double -> Double -> Tension -> Tension -> PathJoin
unmetaJoin !z0 !z1 !theta !phi !alpha !beta
| abs phi < 1e-4 && abs theta < 1e-4 = JoinLine
| otherwise = JoinCurve u v
where Point dx dy = z1^-^z0
bounded = (sf <= 0 && st <= 0 && sf <= 0) ||
(sf >= 0 && st >= 0 && sf >= 0)
rr' = velocity st sf ct cf alpha
ss' = velocity sf st cf ct beta
stf = st*cf + sf*ct -- sin (theta + phi)
st = sin theta
sf = sin phi
ct = cos theta
cf = cos phi
rr = case alpha of
TensionAtLeast _ | bounded ->
min rr' (sf/stf)
_ -> rr'
ss = case beta of
TensionAtLeast _ | bounded ->
min ss' (st/stf)
_ -> ss'
u = z0 ^+^ rr *^ Point (dx*ct - dy*st) (dy*ct + dx*st) -- z0 + rr * (rotate theta chord)
v = z1 ^-^ ss *^ Point (dx*cf + dy*sf) (dy*cf - dx*sf) -- z1 - ss * (rotate (-phi) chord)
constant1, constant2, sqrt2 :: Double
constant1 = 3/2*(sqrt 5 - 1)
constant2 = 3/2*(3 - sqrt 5)
sqrt2 = sqrt 2
-- another magic formula by John Hobby.
velocity :: Double -> Double -> Double
-> Double -> Tension -> Double
velocity st sf ct cf t =
(2 + sqrt2 * (st - sf/16)*(sf - st/16)*(ct - cf)) /
((3 + constant1*ct + constant2*cf) * tensionValue t)