content: -- | Translation of cairo-arc.c
  module Graphics.Rasterific.Arc( Direction( .. ), arcInDirection ) where
  
  import Data.Maybe( fromMaybe )
  import qualified Data.Vector.Unboxed as VU
  
  import Graphics.Rasterific.Transformations
  import Graphics.Rasterific.Linear
  import Graphics.Rasterific.Types
  
  
  errorTable :: VU.Vector (Float, Float)
  errorTable = VU.imap calcAngle errors where
    calcAngle i a = (pi / fromIntegral i, a)
    errors = VU.fromListN 10
      [ 0.0185185185185185036127
      , 0.000272567143730179811158
      , 2.38647043651461047433e-05
      , 4.2455377443222443279e-06
      , 1.11281001494389081528e-06
      , 3.72662000942734705475e-07
      , 1.47783685574284411325e-07
      , 6.63240432022601149057e-08
      , 3.2715520137536980553e-08
      , 1.73863223499021216974e-08
      , 9.81410988043554039085e-09 ]
  
  {- Spline deviation from the circle in radius would be given by:
  
      error = sqrt (x**2 + y**2) - 1
  
     A simpler error function to work with is:
  
      e = x**2 + y**2 - 1
  
     From "Good approximation of circles by curvature-continuous Bezier
     curves", Tor Dokken and Morten Daehlen, Computer Aided Geometric
     Design 8 (1990) 22-41, we learn:
  
      abs (max(e)) = 4/27 * sin**6(angle/4) / cos**2(angle/4)
  
     and
      abs (error) =~ 1/2 * e
  
     Of course, this error value applies only for the particular spline
     approximation that is used in _cairo_gstate_arc_segment.  -}
  fixAngleError :: Int -> Float -> Float
  fixAngleError i tolerance
      | errorNormalized > tolerance = fixAngleError (i + 1) tolerance
      | otherwise = angle
    where
      angle = pi / fromIntegral i
      errorNormalized = 2.0/27.0 * (sin (angle / 4) ** 6) / (cos (angle / 4) ** 2)
  
  arcMaxAngleForToleranceNormalized :: Float -> Float
  arcMaxAngleForToleranceNormalized tolerance = fixAngleError (angleIndex + 1) tolerance
    where
      angleIndex = fromMaybe
        (VU.length errorTable) $
         VU.findIndex ((< tolerance) . snd) errorTable
  
  arcSegmentsNeeded :: Float -> Float -> Transformation -> Float
                    -> Int
  arcSegmentsNeeded angle radius trans tolerance = ceiling (angle / maxAngle) where
    -- the error is amplified by at most the length of the
    -- major axis of the circle; see cairo-pen.c for a more detailed analysis
    -- of this.
    majorAxis = matrixTransformedCircleMajorAxis trans radius
    maxAngle = arcMaxAngleForToleranceNormalized (tolerance / majorAxis)
  
  -- determine the length of the major axis of a circle of the given radius
  -- after applying the transformation matrix.
  matrixTransformedCircleMajorAxis :: Transformation -> Float -> Float
  matrixTransformedCircleMajorAxis (Transformation a c _
                                                   b d _) radius =
       radius * sqrt (f + norm v)
    where
      i = a*a + b*b;
      j = c*c + d*d;
  
      f = 0.5 * (i + j)
      v = V2 (0.5 * (i - j)) (a * c + b * d)
      -- we don't need the minor axis length, which is
      -- double min = radius * sqrt (f - sqrt (g*g+h*h));
  
  -- | Direction of the arc
  data Direction = Forward | Backward
  
  clampAngle :: Float -> Float -> Float
  clampAngle angleMin = go where
    go angleMax
      | angleMax - angleMin > 4 * pi = go $ angleMax - 2 * pi
      | otherwise = angleMax
  
  subdivideAngles :: (Monoid m)
                  => Direction -> (Float -> Float -> m) -> Float -> Float -> m
  subdivideAngles dir f aMin = go aMin . clampAngle aMin where
    go angleMin angleMax | deltaAngle > pi = case dir of
        Forward -> go angleMin angleMid <> go angleMid angleMax
        Backward -> go angleMid angleMax <> go angleMin angleMid
      where
        deltaAngle = angleMax - angleMin
        angleMid = angleMin + deltaAngle / 2
    go angleMin angleMax = f angleMin angleMax
  
  
  {-  We want to draw a single spline approximating a circular arc radius
     R from angle A to angle B. Since we want a symmetric spline that
     matches the endpoints of the arc in position and slope, we know
     that the spline control points must be:
  
  	(R * cos(A), R * sin(A))
  	(R * cos(A) - h * sin(A), R * sin(A) + h * cos (A))
  	(R * cos(B) + h * sin(B), R * sin(B) - h * cos (B))
  	(R * cos(B), R * sin(B))
  
     for some value of h.
  
     "Approximation of circular arcs by cubic poynomials", Michael
     Goldapp, Computer Aided Geometric Design 8 (1991) 227-238, provides
     various values of h along with error analysis for each.
  
     From that paper, a very practical value of h is:
  
  	h = 4/3 * tan(angle/4)
  
     This value does not give the spline with minimal error, but it does
     provide a very good approximation, (6th-order convergence), and the
     error expression is quite simple, (see the comment for
     _arc_error_normalized).
  -}
  arcSegment :: Point -> Float -> Float -> Float -> PathCommand
  arcSegment (V2 xc yc) radius angleA angleB = PathCubicBezierCurveTo p1 p2 p3 where
    rSinA = radius * sin angleA
    rCosA = radius * cos angleA
    rSinB = radius * sin angleB
    rCosB = radius * cos angleB
  
    h = 4.0/3.0 * tan ((angleB - angleA) / 4.0)
  
    p1 = V2 (xc + rCosA - h * rSinA) (yc + rSinA + h * rCosA)
    p2 = V2 (xc + rCosB + h * rSinB) (yc + rSinB - h * rCosB)
    p3 = V2 (xc + rCosB) (yc + rSinB)
  
  -- | Translate an arc with a definition similar to the
  -- one given in Cairo to a list of bezier path command
  arcInDirection :: Point -- ^ center
                 -> Float -- ^ Radius
                 -> Direction
                 -> Float -- ^ Tolerance
                 -> Float -- ^ Angle minimum
                 -> Float -- ^ Angle maximum
                 -> [PathCommand]
  arcInDirection p@(V2 px py) radius dir tolerance 
    | isNaN px || isNaN py || isNaN radius = mempty
    | otherwise = subdivideAngles dir go where
    go angleMin angleMax = commands where
      deltaAngle = angleMax - angleMin
      segmentCount = arcSegmentsNeeded deltaAngle radius mempty tolerance
  
      (angle, angleStep) = case dir of
        Forward -> (angleMin, deltaAngle / fromIntegral segmentCount)
        Backward -> (angleMax, - deltaAngle / fromIntegral segmentCount)
  
      commands =
          [arcSegment p radius a (a + angleStep)
              | i <- [0 .. segmentCount - 1]
              , let a = angle + angleStep * fromIntegral i]
  
  