\begin{code}
{-# OPTIONS -XRecordWildCards -XNamedFieldPuns #-}
-- | This module contains the function to approximate a list of curves with
-- degree 3 Bezier curves, using a least squares method.

module Graphics.Typography.Geometry.Approximation(approximate) where

import qualified Data.Vector.Unboxed as UV
import Graphics.Typography.Geometry.Bezier
import Graphics.Typography.Geometry
import Algebra.Polynomials.Bernstein

import Algebra.Polynomials.Numerical
-- import Debug.Trace
rnd::Interval->Double
rnd (Interval a b)=(a+b)/2


-- | Approximates a list of 'Curves' with a list of degree 3 Bernstein curves.
approximate::[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
approximate []=[]
approximate l0@(h0:_)= -- traceShow "starting" $
  let approx::Double->Double->[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
      approx _ _ []=[]
      approx x0 y0 (cc@(Circle {..}):s)= -- traceShow "circle" $
        let theta=abs $ t1-t0 in
        if theta <= pi/2 then
          let x0_=cos $ theta/2
              y0_=sin $ theta/2
              x1_=(4-x0_)/3
              y1_=(1-x0_)*(3-x0_)/(3*y0_)
        
              c0=cos $! theta/2+t0
              s0=sin $! theta/2+t0
        
        
              px0=c0*x0_ - s0*y0_
              py0=s0*x0_ + c0*y0_
        
              px1=c0*x1_ - s0*y1_
              py1=s0*x1_ + c0*y1_
        
              px2=c0*x1_ + s0*y1_
              py2=s0*x1_ - c0*y1_
        
              -- px3=c0*x0_ + s0*y0_
              -- py3=s0*x0_ - c0*y0_
        
              x1=cx0+(a*px0+b*py0)
              y1=cy0+(c*px0+d*py0)
              
              (Matrix2 a b c d)=matrix
              x=Bernsteinp 4 $ UV.fromList
                [ x0, -- cx0+(a*px3+b*py3),
                  cx0+(a*px2+b*py2),
                  cx0+(a*px1+b*py1),
                  x1]
              
              y=Bernsteinp 4 $ UV.fromList
                [ y0, -- cy0+(c*px3+d*py3),
                  cy0+(c*px2+d*py2),
                  cy0+(c*px1+d*py1),
                  y1 ]
          in
           (x,y):(approx x1 y1 s)
     
        else
          let t1'=(t1+t0)/2 in
          approx x0 y0 $ (cc { t1=t1' }):(cc { t0=t1' }):s
{-
      approx x0 y0 (h@(Bezier{}):s)=
        -- incorrect !
        (restriction (cx h) (t0 h) (t1 h),
         restriction (cy h) (t0 h) (t1 h)):
        (approx (UV.last $ coefs $ cx h)
         (UV.last $ coefs $ cy h) s)
-}
      -- Ce qui suit est une methode de moindres carres
      approx x0 y0 (off_:s)= -- traceShow ("offset") $
        -- On commence par chercher les points ou la derivee de la norme
        -- de la tangente est maximale. C'est la qu'on va couper s'il y
        -- a un probleme.
        let bx=restriction (cx off_) (t0 off_) (t1 off_)
            by=restriction (cy off_) (t0 off_) (t1 off_)
            off=off_ { cx=bx,cy=by,t0=0,t1=1 }
            ibx=elevate (bounds by-bounds bx) $ intervalize bx
            iby=elevate (bounds bx-bounds by) $ intervalize by
            points=
              let np=10 in
              map (\x->(x/np,x/np)) [0..np]
            -- Ensuite, moindres carres standard, comme dans Hoschek 88.
      
            vx0=ibx?1-ibx?0
            vy0=iby?1-iby?0
            vx1=ibx?(bounds ibx-2)-ibx?(bounds ibx-1)
            vy1=iby?(bounds iby-2)-iby?(bounds iby-1)
    
            (wx0,wy0)=evalCurve off 0
            (wx1,wy1)=evalCurve off 1

            wx=Bernsteinp 4 $ UV.fromList [wx0,wx0,wx1,wx1] :: Bernsteinp Int Interval
            wy=Bernsteinp 4 $ UV.fromList [wy0,wy0,wy1,wy1] :: Bernsteinp Int Interval

            bern1=Bernsteinp 4 $ UV.fromList [0,1,0,0] :: Bernsteinp Int Interval
            bern2=Bernsteinp 4 $ UV.fromList [0,0,1,0] :: Bernsteinp Int Interval

            sumAll a b c d x1 y1 ((h1,h2):ss)=
        
              let h=Interval h1 h2
                  (xi,yi)=evalCurve off h
            
                  b1=eval bern1 h
                  b2=eval bern2 h
            
                  a'=a + (vx0*vx0+vy0*vy0)*b1*b1
                  b'=b + (vx0*vx1 + vy0*vy1)*b1*b2
                  c'=c + (vx0*vx1 + vy0*vy1)*b1*b2
                  d'=d + (vx1*vx1+vy1*vy1)*b2*b2
            
                  dx=xi-(eval wx h)
                  dy=yi-(eval wy h)
            
                  x1'=x1 + (vx0*dx + vy0*dy)*b1
                  y1'=y1 + (vx1*dx + vy1*dy)*b2
              in
               sumAll a' b' c' d' x1' y1' ss
         
            sumAll a b c d x1 y1 []=(a,b,c,d,x1,y1)
                                    
            (ra,rb,rc,rd,rx1,ry1)=sumAll 0 0 0 0 0 0 points
            
            (Matrix2 a_ b_ c_ d_)=inverse $ Matrix2 ra rb rc rd
            lambda1=a_*rx1+b_*ry1
            lambda2=c_*rx1+d_*ry1
      
            -- On a la courbe optimale. Il faut chercher ou on va couper, maintenant
            xap=Bernsteinp 4 $ UV.fromList [wx0,
                                            wx0+lambda1*vx0,
                                            wx1+lambda2*vx1,
                                            wx1]
            yap=Bernsteinp 4 $ UV.fromList [wy0,
                                            wy0+lambda1*vy0,
                                            wy1+lambda2*vy1,
                                            wy1]
            (err,arg)=foldl (\m (h1,h2)->
                              let (xi,yi)=evalCurve off (Interval h1 h2)
                                  xj=eval xap (Interval h1 h2)
                                  yj=eval yap (Interval h1 h2)
                              in
                               max m (iup $ abs $ (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj), (h1+h2)/2))
                      (0,0) points
        in
         if err<=0.01 then
           (desintervalize xap,desintervalize yap):(approx (rnd wx1) (rnd wy1) s)
         else
           approx x0 y0 $
            (off { cx=restriction (cx off) 0 arg,
                   cy=restriction (cy off) 0 arg }):
            (off { cx=restriction (cx off) arg 1,
                   cy=restriction (cy off) arg 1 }):s
           
      (x0h,y0h)=evalCurve h0 $ Interval (t0 h0) (t0 h0)
      
  in
   approx (rnd x0h) (rnd y0h) l0
           
desintervalize::(Bernsteinp a Interval)->(Bernsteinp a Double)
desintervalize b=b { coefs=UV.map rnd $ coefs b}
  
\end{code}