\documentclass{article} %include lhs2TeX.fmt \begin{document} \begin{code}
{-# OPTIONS -XUnboxedTuples -XBangPatterns -XNamedFieldPuns -XRecordWildCards -XMagicHash -cpp #-}
-- | This module contains the basic functions for manipulating Bezier curves. It is heavily
-- based on the book by N. M. Patrikalakis and T. Maekawa, Shape Interrogation for Computer
-- Aided Design and Manufacturing.

module Graphics.Typography.Geometry.Bezier (
Curve(..),line,bezier3,
offset,
inter,
evalCurve,distance,
left,bottom,right,top)  where

import Algebra.Polynomials.Bernstein
import Algebra.Polynomials.Numerical
import Graphics.Typography.Geometry
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as UV
import Data.List (partition,sort)

-- | The type for representing all types of curves.
data Curve=
Bezier { cx::Bernsteinp Int Double,
cy::Bernsteinp Int Double,
t0::Double,
t1::Double }

| Offset { cx::Bernsteinp Int Double,
cy::Bernsteinp Int Double,
t0::Double,
t1::Double,
matrix::Matrix2 Double
}
| Circle { cx0::Double,
cy0::Double,
t0::Double,
t1::Double,
matrix::Matrix2 Double
}
deriving (Show)

-- | The basic constructor for lines : a line is a degree 1 Bezier curve
line::Double->Double->Double->Double->Curve
line px py px' py'=Bezier { cx=Bernsteinp 2 $UV.fromList [px,px'], cy=Bernsteinp 2$ UV.fromList [py,py'],
t0=0,t1=1 }

-- | A shortcut to define degree 3 Bezier curves from points. If the control
-- points are @a,b,c,d@, the function should be called with
-- @'bezier3' xa ya xb yb xc yc xd yd@.
bezier3::Double->Double->Double->Double->Double->Double->Double->Double->Curve
bezier3 px0 py0 px1 py1 px2 py2 px3 py3=
Bezier { cx=Bernsteinp 4 $UV.fromList [px0,px1,px2,px3], cy=Bernsteinp 4$ UV.fromList [py0,py1,py2,py3],
t0=0,t1=1 }

instance Geometric Curve where
translate x y cur@(Circle{cx0,cy0})=
cur { cx0=cx0+x,cy0=cy0+y }
translate x y cur=
cur { cx=(cx cur) { coefs=UV.map (+x) $coefs$ cx cur},
cy=(cy cur) { coefs=UV.map (+y) $coefs$ cy cur} }

apply m0@(Matrix2 a b c d) cir@(Circle{cx0,cy0,matrix})=
cir { cx0=a*cx0+b*cy0, cy0=c*cx0+d*cy0, matrix=m0*matrix }
apply (Matrix2 a b c d) cur=
cur { cx=(scale a $cx cur)+(scale b$ cy cur),
cy=(scale c $cx cur)+(scale d$ cy cur) }

-- | Gives the point corresponding to the given value of the parameter
evalCurve::Curve->Interval->(Interval,Interval)
evalCurve (Offset{..}) t=
let ix=intervalize cx
iy=intervalize cy
xt=eval ix t
yt=eval iy t
m@(Matrix2 a b c d)=intervalize matrix
(Matrix2 a_ b_ c_ d_)=inverse m
xt0'=eval (derivate ix) t
yt0'=eval (derivate iy) t
xt'=a_*xt0' + b_*yt0'
yt'=c_*xt0' + d_*yt0'
dd=sqrt $xt'*xt' + yt'*yt' in (xt+(a*yt'-b*xt')/dd, yt+(c*yt'-d*xt')/dd) evalCurve (Circle{..}) alpha= let xx=cos alpha yy=sin alpha (Matrix2 a b c d)=intervalize matrix in (interval cx0+a*xx+b*yy, interval cy0+c*xx+d*yy) evalCurve (Bezier{..}) t= let ix=intervalize cx iy=intervalize cy xx=eval ix t yy=eval iy t in (xx,yy) data Topo=Dehors | SurLaLigne | Dedans deriving Eq -- | @'inter' c0 c1@ is a list of all possible points of intersection -- between curves @c0@ and @c1@ : if @(u,v,w,x)@ is returned by 'inter', -- then curve @c0@ may intersect with @c1@ between parameter values @u@ -- and @v@, which corresponds to parameter values between @w@ and @x@ for -- @c1@. The implementation guarantees that all actual solutions are found, -- but possibly false solutions may also be returned. inter::Curve->Curve->[((Double,Double,Double,Double))] inter op@(Offset { cx=bxp_,cy=byp_,matrix=mp,t0=t0a,t1=t1a }) (Offset { cx=bxq_,cy=byq_,matrix=mq,t0=t0b,t1=t1b })= -- Attention : verifier si c'est la meme generatrice let thrx=1e-5 solutions _ []=[] solutions thr boxes@(_:_)= let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2,eq3])) boxes (correct,toRefine)=partition (\(u,v,_,_,_,_,_,_)-> let (xu,yu)=evalCurve op (Interval u u) (xv,yv)=evalCurve op (Interval v v) in (iup$ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
in
correct++(solutions (thr/2) toRefine)
in
map (\(u,v,w,x,_,_,_,_)->(u,v,w,x)) $solutions 1e-2$
[(t0a,t1a,t0b,t1b,0,1,0,1)::
(Double,Double,Double,Double,Double,Double,Double,Double)]

where

imp@(Matrix2 ap bp cp dp)=intervalize mp
imq@(Matrix2 aq bq cq dq)=intervalize mq
(Matrix2 ap_ bp_ cp_ dp_)=inverse imp
(Matrix2 aq_ bq_ cq_ dq_)=inverse imq

bxp=intervalize bxp_
byp=intervalize byp_
bxq=intervalize bxq_
byq=intervalize byq_

bxp4=promote 1 bxp
byp4=promote 1 byp
bxq4=promote 2 bxq
byq4=promote 2 byq

bxp'=derivate bxp
byp'=derivate byp
bxq'=derivate bxq
byq'=derivate byq

bXp'=promote 1 $(scale ap_ bxp')+(scale bp_ byp') bYp'=promote 1$ (scale cp_ bxp')+(scale dp_ byp')

bXq'=promote 2 $(scale aq_ bxq')+(scale bq_ byq') bYq'=promote 2$ (scale cq_ bxq')+(scale dq_ byq')

bomp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp') :: Bernsteinp (Int,Int,Int,Int) Interval
bomq@(Bernsteinp _ omegaq)=(bXq'*bXq')+(bYq'*bYq') :: Bernsteinp (Int,Int,Int,Int) Interval

au=
let au_=sqrt $UV.minimum$ UV.map ilow omegap in
max 0 $fpred au_ bu= let bu_=sqrt$ UV.maximum $UV.map iup omegap in fsucc bu_ av= let av_=sqrt$ UV.minimum $UV.map ilow omegaq in max 0$ fpred av_
bv=
let bv_=sqrt $UV.maximum$ UV.map iup omegaq in
fsucc bv_

alphau=Bernsteinp (1,1,2,1) $UV.fromList [Interval au au,Interval bu bu] alphav=Bernsteinp (1,1,1,2)$ UV.fromList [Interval av av,Interval bv bv]

eq0=
((bxp4*alphau*alphav) + (scale ap $bYp'*alphav) - (scale bp$ bXp'*alphav)
-(bxq4*alphau*alphav) - (scale aq $bYq'*alphau) + (scale bq$ bXq'*alphau))
eq1=
((byp4*alphau*alphav) + (scale cp $bYp'*alphav) - (scale dp$ bXp'*alphav)
-(byq4*alphau*alphav) - (scale cq $bYq'*alphau) + (scale dq$ bXq'*alphau))
eq2=bomp-(alphau*alphau)
eq3=bomq-(alphav*alphav)

inter b@(Circle{}) a@(Offset{})=
map (\(i,j,k,l)->(k,l,i,j)) $inter a b inter o@(Offset { cx=bxp, cy=byp, matrix=mp }) cir@(Circle{cx0,cy0,matrix=mq})= let ix=intervalize bxp iy=intervalize byp m@(Matrix2 a b c d)=intervalize mp (Matrix2 a_ b_ c_ d_)=inverse m x'=derivate ix y'=derivate iy xx'=(scale a_ x')+(scale b_ y') yy'=(scale c_ x')+(scale d_ y') omega@(Bernsteinp _ omegap)=xx'*xx'+yy'*yy' au= let au_=sqrt$ UV.minimum $UV.map ilow omegap in max 0$ fpred au_
bu=
let bu_=sqrt $UV.maximum$ UV.map iup omegap in
fsucc bu_

alphau=Bernsteinp (1,2) $UV.fromList [Interval au au,Interval bu bu] lambda=(promote 1 omega) - alphau*alphau -- Avant multiplication par M_C^-1 xx0=(promote 1$ ix-(intervalize $constant cx0))*alphau +(promote 1$ scale a yy'-scale b xx')
yy0=(promote 1 $iy-(intervalize$ constant cy0))*alphau
+(promote 1 $scale c yy'-scale d xx') :: Bernsteinp (Int,Int) Interval (Matrix2 ac_ bc_ cc_ dc_)=inverse$ intervalize mq
xx1=(scale ac_ xx0)+(scale bc_ yy0)
yy1=(scale cc_ xx0)+(scale dc_ yy0)

eqc=xx1*xx1+yy1*yy1-alphau*alphau

thrx=1e-5

solutions _ []=[]
solutions thr boxes@(_:_)=
let sol0=concatMap (solve thr (V.fromList [eqc,lambda])) boxes

(correct,toRefine)=partition (\(u,v,_,_)->
let (xu,yu)=evalCurve o (Interval u u)
(xv,yv)=evalCurve o (Interval v v)
in
(iup $(xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0 in correct++(solutions (thr/2) toRefine) -- Removing false positives by computing the distance to the center of -- the circle (this is quite fast). removeFalse cl0 (h@(_,v,_,_):h'@(u',_,_,_):s)= let u''=(v+u')/2 (xu,yu)=evalCurve o (Interval u'' u'') Interval dl du=distance xu yu cir cl1 | du<1 = Dedans | dl>1 = Dehors | otherwise = SurLaLigne in if cl0/=cl1 then h:(removeFalse cl1 (h':s)) else removeFalse cl1 (h':s) removeFalse _ l=l initCl= let (x0,y0)=evalCurve o (Interval (t0 o) (t0 o)) Interval dl du=distance x0 y0 cir in if dl>1 then Dehors else if du<1 then Dedans else SurLaLigne in foldl (\l (u,v,_,_)-> let (Interval xl xu,Interval yl yu)=evalCurve o (Interval u v) in case angle (Interval xl xu) (Interval yl yu) cir of Just (Interval a0 a1)-> (u,v,a0,a1):l Nothing->l ) []$ removeFalse initCl $sort$ solutions 1e-2 [(t0 o,t1 o,0::Double,1::Double)]

inter a@(Circle{cx0=x0a,cy0=y0a,matrix=ma})
b@(Circle{cx0=x0b,cy0=y0b,matrix=mb})=

if (intervalize ma)intersects(intervalize mb) && x0a==x0b && y0a==y0b then
let up ix@(Interval _ x_) tt0 tt1
| x_<tt0 =
up (ix+(2*interval pi)) tt0 tt1
| otherwise = down ix tt0 tt1
down ix@(Interval x_ x__) tt0 tt1
| x_>tt1 =
down (ix-(2*interval pi)) tt0 tt1
| x__<tt0 =
Nothing
| otherwise =
Just ix

alpha=up (interval $t0 a) (t0 b) (t1 b) beta=up (interval$ t0 b) (t0 b) (t1 b)
in

case (alpha,beta) of
(Just aa,Just ab)->
case (up aa (t0 a) (t1 a),
up ab (t0 a) (t1 a)) of

(Just ba,Just bb)
| ilow aa<=iup ab -> [(ilow ba, iup bb,
ilow aa, iup ab)]
| otherwise->
case (up (interval $t0 b) (t0 a) (t1 a), up (interval$ t1 b) (t0 a) (t1 a)) of
(Just b0,Just b1)->
[(ilow b0,iup bb,
t0 b, iup ab),
(ilow ba,iup b1,
ilow aa, t1 b)]
_->[]
_->[]
_->[]

else
let thr=1e-5
solutions=solve thr (V.fromList [eq0,eq1]) (fpred u0,fsucc v0,
fpred w0,fsucc x0)
in
foldl (\l (u,v,w,x)->
let alpha=angle (Interval u v) (Interval w x) a
beta=angle (Interval u v) (Interval w x) b
in
case alpha of
Just (Interval a0l a0u)->
case beta of
Just (Interval b0l b0u)->(a0l,a0u,b0l,b0u):l
_->l
_->l
) [] solutions
where

ima@(Matrix2 am bm cm dm)=intervalize ma

maxa=max (iup $abs am+abs bm) (iup$ abs cm+abs dm)
(u0,v0,w0,x0)=(x0a-maxa,x0a+maxa,y0a-maxa,y0a+maxa)

-- x-x0
xxa0=intervalize $Bernsteinp (2,1)$ UV.fromList [-x0a,1-x0a] :: Bernsteinp (Int,Int) Interval
yya0=intervalize $Bernsteinp (1,2)$ UV.fromList [-y0a,1-y0a] :: Bernsteinp (Int,Int) Interval
(Matrix2 aa_ ba_ ca_ da_)=inverse ima
xxa=(scale aa_ xxa0)+(scale ba_ yya0)::Bernsteinp (Int,Int) Interval
yya=(scale ca_ xxa0)+(scale da_ yya0)

xxb0=intervalize $Bernsteinp (2,1)$ UV.fromList [-x0b,1-x0b]
yyb0=intervalize $Bernsteinp (1,2)$ UV.fromList [-y0b,1-y0b]
(Matrix2 ab_ bb_ cb_ db_)=inverse $intervalize mb xxb=(scale ab_ xxb0)+(scale bb_ yyb0) yyb=(scale cb_ xxb0)+(scale db_ yyb0) c1=Bernsteinp (1,1)$ UV.singleton 1

eq0=xxa*xxa+yya*yya-c1
eq1=xxb*xxb+yyb*yyb-c1

inter op@(Bezier{cx=bxa,cy=bya,t0=t0a,t1=t1a}) (Bezier{cx=xb,cy=yb,t0=t0b,t1=t1b})=

let p0=(promote 1 $intervalize bxa)-(promote 2$ intervalize xb) :: Bernsteinp (Int,Int) Interval
p1=(promote 1 $intervalize bya)-(promote 2$ intervalize yb) :: Bernsteinp (Int,Int) Interval
thrx=1e-2
solutions _ []=[]
solutions thr boxes@(_:_)=
let sol0=concatMap (solve thr (V.fromList [p0,p1])) boxes

(correct,toRefine)=partition (\(u,v,_,_)->
let (xu,yu)=evalCurve op (Interval u u)
(xv,yv)=evalCurve op (Interval v v)
in
(iup $(xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0 in correct++(solutions (thr/2) toRefine) in solutions 1e-2 [(t0a,t1a,t0b,t1b)] inter cir@(Circle{}) bez@(Bezier{})=map (\(u,v,w,x)->(w,x,u,v))$ inter bez cir

inter bez@(Bezier{}) cir@(Circle{})=
let xx=(intervalize $cx bez)-(intervalize$ constant $cx0 cir) yy=(intervalize$ cy bez)-(intervalize $constant$ cy0 cir)
(Matrix2 a b c d)=inverse $intervalize$ matrix cir
xx0=scale a xx+scale b yy
yy0=scale c xx+scale d yy

thrx=1e-5

solutions _ []=[]
solutions thr boxes@(_:_)=
let sol0=concatMap (solve thr (V.singleton (xx0*xx0+yy0*yy0-(constant 1)))) boxes

(correct,toRefine)=partition (\(u,v)->
let (xu,yu)=evalCurve bez (Interval u u)
(xv,yv)=evalCurve bez (Interval v v)
in
(iup $(xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0 in correct++(solutions (thr/2) toRefine) in foldl (\l (u,v)-> let (Interval xl xu,Interval yl yu)=evalCurve bez (Interval u v) in case angle (Interval xl xu) (Interval yl yu) cir of Just (Interval a0 a1)-> (u,v,a0,a1):l Nothing->l ) []$!
solutions (1e-2) [(t0 bez,t1 bez)]

inter bez@(Bezier{}) off@(Offset{})=map (\(u,v,w,x)->(w,x,u,v)) $inter off bez inter off@(Offset{}) bez@(Bezier{})= let thr=1e-2 thrx=1e-5 solutions _ []=[] solutions thr boxes@(_:_)= let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2])) boxes (correct,toRefine)=partition (\(u,v,_,_,_,_)-> let (xu,yu)=evalCurve off (Interval u u) (xv,yv)=evalCurve off (Interval v v) in (iup$ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
in
correct++(solutions (thr/2) toRefine)
in
map (\(a,b,c,d,_,_)->(a,b,c,d)) $solutions 1e-2$
[(0,1,0,1,0,1)::(Double,Double,Double,Double,Double,Double)]
where

bxp=intervalize $cx off byp=intervalize$ cy off
bxq=intervalize $cx bez byq=intervalize$ cy bez

bxp'=derivate bxp
byp'=derivate byp
bxp3=promote 1 bxp
byp3=promote 1 byp
bxq3=promote 2 bxq
byq3=promote 2 byq

mp@(Matrix2 ap bp cp dp)=intervalize $matrix$ off
(Matrix2 ap_ bp_ cp_ dp_)=inverse mp
bXp'=promote 1 $(scale ap_ bxp')+(scale bp_ byp') bYp'=promote 1$ (scale cp_ bxp')+(scale dp_ byp')

omp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp')
au=
let au_=sqrt $UV.minimum$ UV.map ilow omegap in
max 0 $fpred au_ bu= let bu_=sqrt$ UV.maximum $UV.map iup omegap in fsucc bu_ alphau=Bernsteinp (1,1,2)$ UV.fromList [Interval au au,Interval bu bu]
eq0=bxp3*alphau + (scale ap bYp') - (scale bp bXp') - bxq3
eq1=byp3*alphau + (scale cp bYp') - (scale dp bXp') - byq3
eq2=alphau*alphau-omp

angle::Interval->Interval->Curve->Maybe Interval
angle x y (Circle { cx0,cy0,matrix,t0,t1 })=
let vx=x-interval cx0
vy=y-interval cy0

Matrix2 a b c d=inverse $intervalize matrix -- L'arithmetique d'intervalles fait un peu n'importe quoi -- quand le vecteur est trop long. On le raccourcit. alpha= let co@(Interval col cou)=a*vx+b*vy Interval sil siu=c*vx+d*vy co2= let (col2,cou2)=if col*col<cou*cou then (col*col,cou*cou) else (cou*cou,col*col) in Interval (fpred col2) (fsucc cou2) si2= let (sil2,siu2)=if sil*sil<siu*siu then (sil*sil,siu*siu) else (siu*siu,sil*sil) in Interval (fpred sil2) (fsucc siu2) coco=co/(sqrt (co2+si2)) ac@(Interval acl acu)=acos$ Interval (max (-1) $ilow coco) (min 1$ iup coco)
in
if siu<0 then negate ac else
if sil>=0 then ac else
Interval (negate $min (abs acl) (abs acu)) (max (abs acl) (abs acu)) up ix | iup ix<t0 = up$ ix+(2*interval pi)
| otherwise =
down ix
down ix
| ilow ix>t1 =
down $ix-(2*interval pi) | iup ix<t0 = Nothing | otherwise = Just ix in up alpha angle _ _ _=error "angle" -- | Pseudo-distance from a point to a curve. Is the result is -- smaller than 1, the point is inside the curve. If it is greater -- than 1, the point is outside. Else we don't know (as usual with -- interval arithmetic). distance::Interval->Interval->Curve->Interval distance x0 y0 (Bezier{..})= distance x0 y0 (Offset{cx,cy,t0,t1,matrix=Matrix2 1 0 0 1}) distance x0 y0 (Offset{..})= let (Matrix2 a b c d)=inverse$ intervalize matrix
vx_=intervalize cx-(constant x0)
vy_=intervalize cy-(constant y0)
vx=scale a vx_+scale b vy_
vy=scale c vx_+scale d vy_

dist=vx*vx+vy*vy
in
foldl (\di (u,v)->let di'=eval dist (Interval u v) in
if iup di<iup di' then di else di') (Interval (1/0) (1/0)) $(t0,t0):(t1,t1):(solve 1e-5 (V.singleton (derivate dist)) (t0,t1)) distance x1 y1 (Circle{..})= let (Matrix2 a b c d)=inverse$ intervalize matrix
vx_=x1-Interval cx0 cx0
vy_=y1-Interval cy0 cy0
vx=a*vx_+b*vy_
vy=c*vx_+d*vy_
in
vx*vx+vy*vy

-- | Offsets a given Bezier curve with the given pen matrix. The original
-- pen is a circle of radius one, the matrix, if inversible, is applied to it.

offset::Matrix2 Double->Curve->[Curve]
offset m (Bezier{cx=x@(Bernsteinp nx bx),cy=y@(Bernsteinp ny by)})=
if nx <=1 && ny <=1 then
[Circle { cx0=UV.head bx,cy0=UV.head by,t0=ilow 0,t1=iup $2*pi,matrix=m }] else [ c0,c1,c2,c3 ] where im=intervalize m (Matrix2 a_ b_ c_ d_)=inverse im ibx=intervalize x iby=intervalize y lastCoef (Bernsteinp n c) | n>=1 = UV.last c | otherwise = 0 firstCoef (Bernsteinp n c) | n>=1 = UV.head c | otherwise = 0 -- Premiere courbe offset c0=Offset { cx=x, cy=y, t0=0,t1=1,matrix=m } -- Demi-cercle 1 ibx'=derivate ibx iby'=derivate iby -- Calcul du vecteur tangent au bout du premier alpha0= let xx0=lastCoef ibx' yy0=lastCoef iby' xx0'=a_*xx0+b_*yy0 yy0'=c_*xx0+d_*yy0 norm0=sqrt$ xx0'*xx0'+yy0'*yy0'

xx'=xx0'/norm0
yy'=yy0'/norm0
in
if ilow xx'>=0 then
-(acos yy')
else
if iup xx'<=0 then
acos yy'
else
let Interval u v=acos yy' in
Interval (negate $max (abs u) (abs v)) (max (abs u) (abs v)) alpha0'=alpha0+interval pi c1=Circle { cx0=lastCoef x, cy0=lastCoef y, t0=ilow alpha0, t1=iup alpha0', matrix=m } -- Deuxieme courbe offset c2=Offset { cx=reorient x, cy=reorient y, t0=0,t1=1, matrix=m } -- Deuxieme demi-cercle alpha1= let xx0=firstCoef ibx' yy0=firstCoef iby' xx0'=a_*xx0+b_*yy0 yy0'=c_*xx0+d_*yy0 norm0=sqrt$ xx0'*xx0'+yy0'*yy0'

xx'=xx0'/norm0
yy'=yy0'/norm0
in
if ilow xx'>=0 then
-(acos yy')
else
if iup xx'<=0 then
acos yy'
else
let Interval u v=acos yy' in
Interval (negate $max (abs u) (abs v)) (max (abs u) (abs v)) alpha1'=alpha1-pi c3=Circle { cx0=firstCoef x, cy0=firstCoef y, t0=ilow alpha1', t1=iup alpha1, matrix=m } offset _ _=error "offset : undefined yet for other than Bezier" rnd::Interval->Double rnd (Interval a b)=(a+b)/2 derivRoots::Double->Curve->([(Double,Double)],[(Double,Double)]) derivRoots thr (Bezier{..})= (solve thr (V.singleton$ derivate $intervalize cx) (t0,t1), solve thr (V.singleton$ derivate $intervalize cy) (t0,t1)) derivRoots thr (Offset{..})= let ix=intervalize cx iy=intervalize cy x'=derivate ix y'=derivate iy m@(Matrix2 a b c d)=intervalize matrix (Matrix2 a_ b_ c_ d_)=inverse m xx'=(scale a_ x')+(scale b_ y') yy'=(scale c_ x')+(scale d_ y') xx''=derivate xx' yy''=derivate yy' omega=xx'*xx'+yy'*yy' au= let au_=sqrt$ UV.minimum $UV.map ilow$ coefs omega in
max 0 $fpred au_ bu= let bu_=sqrt$ UV.maximum $UV.map iup$ coefs omega in
fsucc bu_
alphau=Bernsteinp (1,2) $UV.fromList [Interval au au,Interval bu bu] eqx0=(promote 1 yy'')*(alphau^(2::Int))-(promote 1$ yy'*yy''+xx'*xx'')
eqy0=(promote 1 $yy'*yy''+xx'*xx'')-(promote 1 xx'')*(alphau^(2::Int)) eqx=(promote 1 x')*(alphau^(3::Int))+(scale a eqx0)+(scale b eqy0) :: Bernsteinp (Int,Int) Interval eqy=(promote 1 y')*(alphau^(3::Int))+(scale c eqx0)+(scale d eqy0) eq=(promote 1 omega)-alphau^(2::Int) :: Bernsteinp (Int,Int) Interval in (map (\(u,v,_,_)->(u,v))$ solve thr (V.fromList [eqx,eq])
((t0,t1,0,1)::(Double,Double,Double,Double)),
map (\(u,v,_,_)->(u,v)) $solve thr (V.fromList [eqy,eq]) ((t0,t1,0,1)::(Double,Double,Double,Double))) derivRoots _ (Circle{..})= let (Matrix2 a b c d)=intervalize matrix aa=sqrt$ a*a+b*b
cc=sqrt $c*c+d*d sx | ilow a>=0 = acos$ b/aa
| otherwise = negate $acos$ b/aa
sy
| ilow c>=0 = acos $d/cc | otherwise = negate$ acos $d/cc Interval ux vx=(-sx)-pi/2 Interval uy vy=(-sy)-pi/2 in ([(ux,vx)],[(uy,vy)]) -- | The leftmost point on a curve left::Curve->(Double,Double) left cur= let (x,y)= foldl (\m@(xx,_) (s,t)-> let m'@(xx',_)=evalCurve cur$ interval $(s+t)/2 in if ilow xx<ilow xx' then m else m') (1/0,1/0)$
(t0 cur,t0 cur):(t1 cur,t1 cur):(fst $derivRoots 1e-5 cur) in (rnd x,rnd y) -- | The bottommost point on a curve bottom::Curve->(Double,Double) bottom cur= let (x,y)= foldl (\m@(_,yy) (s,t)-> let m'@(_,yy')=evalCurve cur$ interval $(s+t)/2 in if ilow yy<ilow yy' then m else m') (1/0,1/0)$
(t0 cur,t0 cur):(t1 cur,t1 cur):(snd $derivRoots 1e-5 cur) in (rnd x,rnd y) -- | The rightmost point on a curve right::Curve->(Double,Double) right cur= let (x,y)= foldl (\m@(xx,_) (s,t)-> let m'@(xx',_)=evalCurve cur$ interval $(s+t)/2 in if iup xx>iup xx' then m else m') (-1/0,-1/0)$
(t0 cur,t0 cur):(t1 cur,t1 cur):(fst $derivRoots 1e-5 cur) in (rnd x,rnd y) -- | The topmost point on a curve top::Curve->(Double,Double) top cur= let (x,y)= foldl (\m@(_,yy) (s,t)-> let m'@(_,yy')=evalCurve cur$ interval $(s+t)/2 in if iup yy>iup yy' then m else m') (-1/0,-1/0)$
(t0 cur,t0 cur):(t1 cur,t1 cur):(snd \$ derivRoots 1e-5 cur)
in
(rnd x,rnd y)

\end{code}