\begin{code}
{-# OPTIONS -XUnboxedTuples -cpp -XRecordWildCards -XNamedFieldPuns -XBangPatterns -XMagicHash -XScopedTypeVariables #-}
-- | This module contains the necessary calls to the other modules of Metafont'
-- to compute the outlines of a given number of pen strokes. The normal way of
-- using it is by calling 'outlines'. One other possible way would be :
--
-- @
-- let curves=cutAll curvesList in
-- remerge $contour curves$ intersections curves
-- @

module Graphics.Typography.Geometry.Outlines (cutAll, intersections, contour, remerge, outlines) where

import Algebra.Polynomials.Bernstein
import Algebra.Polynomials.Numerical
import Graphics.Typography.Geometry.Bezier
import Graphics.Typography.Geometry
import Data.List (sort)
import qualified Data.Map as M
import qualified Data.Vector as V

import Control.Parallel

(!)::V.Vector a->Int->a
(!)=(V.!)

-- | Cuts a curve into a list of consecutive non-selfintersecting curves.
cutNoSelf::Curve->[Curve]
cutNoSelf c@(Circle{})=[c]
cutNoSelf bez@(Bezier{..})=
let ix=intervalize cx
dx=derivate ix
solutions=
sort $filter (\(s,t)->(ilow$ eval ix (Interval s s))*
(iup $eval ix (Interval t t)) <= 0)$
solve 1e-10 (V.singleton dx) (t0,t1)
roots lastU []=
if lastU>=t1 then
[]
else
[bez { t0=lastU }]
roots lastU (u:s)
| u<=lastU = roots lastU s -- on ne coupe pas au debut
| otherwise =
(bez { t0=lastU, t1=u }):
(roots u s)
in
roots t0 $map (\(s,t)->(s+t)/2) solutions cutNoSelf off@(Offset{..})= -- offset let thr=1e-2 ix=intervalize cx iy=intervalize cy x'=derivate ix y'=derivate iy (Matrix2 a b c d)=intervalize matrix (Matrix2 a_ b_ c_ d_)=inverse$ intervalize matrix

xx'=(scale a_ x')+(scale b_ y')
yy'=(scale c_ x')+(scale d_ y')

xx''=derivate xx'
yy''=derivate yy'

evalC (t::Interval)=
let norm=sqrt $(eval xx' t)*(eval xx' t)+(eval yy' t)*(eval yy' t) derx=(eval yy'' t)/norm - ((eval yy' t)*((eval xx' t)*(eval xx'' t)+ (eval yy' t)*(eval yy'' t)))/(norm*norm*norm) dery=(eval xx'' t)/norm - ((eval xx' t)*((eval xx' t)*(eval xx'' t)+ (eval yy' t)*(eval yy'' t)))/(norm*norm*norm) in ((eval x' t)+(a*derx-b*dery), (eval y' t)+(c*derx-d*dery)) zerosx= let verif t lastxx | t>=t1 = [] | otherwise = let (xx,_)=evalC (Interval t t) in if (iup$ xx*lastxx)<=0 then
t:verif (t+thr) xx
else
verif (t+thr) xx

in
verif t0 $fst$ evalC (Interval t0 t0)

roots lastU []=
if lastU>=t1 then
[]
else
[off { t0=lastU }]
roots lastU (u:s)
| u<=lastU = roots lastU s -- on ne coupe pas au debut
| otherwise =
(off { t0=lastU, t1=u }):
(roots u s)
in
roots t0 zerosx

-- | @'cutAll' curves@ is the array of all the curves, cut such that
-- each part does not intersect itself.
cutAll::[[Curve]]->V.Vector (V.Vector Curve)
cutAll l=V.fromList $map (\c->V.fromList$ concatMap cutNoSelf c) l

data Topology=Dedans | SurLaLigne | Dehors deriving (Eq, Ord, Show)

minsert::Ord a=>a->b->M.Map a [b]->M.Map a [b]
minsert x y m=M.insertWith' (++) x [y] m

munion::Ord a=>M.Map a [b]->M.Map a [b]->M.Map a [b]
munion=M.unionWith (++)

mdeleteFindMin::Ord a=>M.Map a [b]->(Maybe (a,b),M.Map a [b])
mdeleteFindMin m=
if M.null m then
(Nothing, m)
else
let ((a,b),m')=M.deleteFindMin m in
case b of
[]->mdeleteFindMin m'
(h:s)->(Just (a,h), if null s then m' else M.insert a s m')

-- | Computes the intersections between any pair of curves given
-- as input, in parallel in GHC using @+RTS -N@.
intersections::V.Vector (V.Vector Curve)->
M.Map (Int,Int,Double) [(Int,Int,Double,Double)]
intersections curves=
let interAll ci cj
| ci>=V.length curves = M.empty
| cj>=V.length curves = interAll (ci+1) (ci+1)
| otherwise =
-- traceShow (ci,i,cj,j) $let next=interAll ci (cj+1) inters | ci==cj = V.ifoldl' (\s0 i curvei-> V.ifoldl' (\s1 j curvej-> foldl (\s2 (ti,ti',tj,tj')-> minsert (ci,i,ti) (cj,j+i+1,tj,tj')$
minsert (cj,j+i+1,tj) (ci,i,ti,ti') $s2) s1$
inter curvei curvej
)
s0 $V.drop (i+1) (curves!cj) ) M.empty$ V.take (V.length (curves!ci)-1) (curves!ci)
| otherwise =
V.ifoldl'
(\s0 i curvei->
V.ifoldl'
(\s1 j curvej->
foldl (\s2 (ti,ti',tj,tj')->
minsert (ci,i,ti) (cj,j,tj,tj') $minsert (ci,i,ti') (cj,j,tj,tj')$
minsert (cj,j,tj) (ci,i,ti,ti') $minsert (cj,j,tj') (ci,i,ti,ti')$ s2) s1 $inter curvei curvej ) s0 (curves!cj) ) M.empty$ V.take (V.length (curves!ci)-1) (curves!ci)
in
(nextparinters)seq
(nextmunioninters)
in
interAll 0 0

-- | 'contour' takes the curves and the intersections computed as in 'intersections',
-- and outputs a list of all simple closed paths defined by the curves in the input.
contour::V.Vector (V.Vector Curve)->
M.Map (Int,Int,Double) [(Int,Int,Double,Double)]->
[[(Int,Int,Double,Double)]]
contour curves inters0=

let allPaths inters1 passages1=
let (first,inters2)=mdeleteFindMin inters1 in
case first of
Nothing->[]
Just ((ci0,i0,ti0),(cj0,j0,tj0a,tj0b))->
--traceShow ("new path",pi0,pj0) $let walk ci i tia tib inters passages= --traceShow ("point",ci,i,tia,tib)$ traceShow (inters) $let (a,b)=M.split (ci,i,tib) inters (next,b')=mdeleteFindMin b in case next of Nothing-> -- traceShow ("echec 1")$
([],a,passages)
Just ((ci',i',ti'),(cj,j,tja,tjb))
| ci==ci0 && i==i0 && (ci',i',ti')>=(ci,i,ti0)->
-- fin du chemin
([(ci,i,tia,ti0)],amunionb',passages)

| ci==ci' && i==i' ->
let isVisible=
let tt=(tia+ti')/2
(xi,yi)=evalCurve (curves!ci!i) (Interval tt tt)
in
V.foldl (\vis cur->
vis &&
iup (distance xi yi $(cur!0) {t0=0,t1=1})>=1) True curves in if (not isVisible) then --traceShow ("invisible",pi')$
([],amunionb',passages)
else
let (_,p1)=M.split (ci,i,ti') passages in
(not $M.null p1) && (let ((ci_,i_,_),ti'_)=M.findMin p1 in ci_==ci && i_==i && ti'_<=ti') in if alreadyPassed then --traceShow ("already passed",pi')$
([],amunionb',passages)
else
--traceShow ("trying",pi') $let (nextPath,nextInters,nextPassages)= walk cj j tja tjb (amunionb')$
M.insert (ci,i,ti') tia passages
in
if null nextPath then
walk ci i tia tib (amunionb') passages
else
((ci,i,tia,ti'):nextPath,
nextInters,
M.insert (ci,i,ti') tia nextPassages)
| otherwise -> --traceShow ("echec 2",ci',i',ti') $([],inters,passages) (path,inters3,passages1')=walk cj0 j0 tj0a tj0b inters2 passages1 in if null path then --traceShow ("abandon")$
allPaths inters3 passages1'
else
--traceShow ("reussi") $path:(allPaths inters3 passages1') in allPaths inters0 M.empty -- | 'remerge' takes the curves, the output of 'contour', and outputs -- the list of "remerged" curves, i.e. where the parts free of self-intersections -- are glued back to each other. remerge::V.Vector (V.Vector Curve)->[(Int,Int,Double,Double)]->[Curve] remerge _ []=[] remerge curves [(ci,i,ti0,ti1)]=[(curves!ci!i) { t0=ti0,t1=ti1 }] remerge curves (l@((ci,i,ti0,_):s))= let (cj,j,_,tj1)=last s in if ci==cj && j+1==i && tj1==ti0 then -- dans ce cas, le dernier est colle au premier let takeFirsts []=(# [],[] #) takeFirsts ((h@(ci',_,_,_)):ss) | ci'==ci = let (# u,v #)=takeFirsts ss in (# h:u, v #) | otherwise = (# [],h:ss #) (# uu,vv #)=takeFirsts l in remerge_$ vv++uu
else
remerge_ l

where
remerge_ []=[]
remerge_ [(cj,j,tj0,tj1)]=[(curves!cj!j) { t0=tj0,t1=tj1 }]
remerge_ ((cj,j,tj0,tj1):(cck@(ck,k,tk0,_)):ss)
| cj==ck && k==j+1 && tj1==tk0 =
let (h':s')=remerge_ $cck:ss in (h' { t0=tj0 }) : s' | otherwise = ((curves!cj!j) { t0=tj0,t1=tj1 }) : (remerge_$ cck:ss)

-- | Takes a list of curves, potentially offset, and outputs the relevants part
-- of the outlines.
outlines::[[Curve]]->[[Curve]]
outlines curves=
let curves'=cutAll curves in
map (remerge curves') $contour curves'$ intersections curves'

\end{code}