-- | 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 :: Vector (Float, Float)
errorTable = (Int -> Float -> (Float, Float))
-> Vector Float -> Vector (Float, Float)
forall a b.
(Unbox a, Unbox b) =>
(Int -> a -> b) -> Vector a -> Vector b
VU.imap Int -> Float -> (Float, Float)
forall a a b. (Floating a, Integral a) => a -> b -> (a, b)
calcAngle Vector Float
errors where
  calcAngle :: a -> b -> (a, b)
calcAngle a
i b
a = (a
forall a. Floating a => a
pi a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i, b
a)
  errors :: Vector Float
errors = Int -> [Float] -> Vector Float
forall a. Unbox a => Int -> [a] -> Vector a
VU.fromListN Int
10
    [ Float
0.0185185185185185036127
    , Float
0.000272567143730179811158
    , Float
2.38647043651461047433e-05
    , Float
4.2455377443222443279e-06
    , Float
1.11281001494389081528e-06
    , Float
3.72662000942734705475e-07
    , Float
1.47783685574284411325e-07
    , Float
6.63240432022601149057e-08
    , Float
3.2715520137536980553e-08
    , Float
1.73863223499021216974e-08
    , Float
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 :: Int -> Float -> Float
fixAngleError Int
i Float
tolerance
    | Float
errorNormalized Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
> Float
tolerance = Int -> Float -> Float
fixAngleError (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Float
tolerance
    | Bool
otherwise = Float
angle
  where
    angle :: Float
angle = Float
forall a. Floating a => a
pi Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i
    errorNormalized :: Float
errorNormalized = Float
2.0Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/Float
27.0 Float -> Float -> Float
forall a. Num a => a -> a -> a
* (Float -> Float
forall a. Floating a => a -> a
sin (Float
angle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
4) Float -> Float -> Float
forall a. Floating a => a -> a -> a
** Float
6) Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ (Float -> Float
forall a. Floating a => a -> a
cos (Float
angle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
4) Float -> Float -> Float
forall a. Floating a => a -> a -> a
** Float
2)

arcMaxAngleForToleranceNormalized :: Float -> Float
arcMaxAngleForToleranceNormalized :: Float -> Float
arcMaxAngleForToleranceNormalized Float
tolerance = Int -> Float -> Float
fixAngleError (Int
angleIndex Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Float
tolerance
  where
    angleIndex :: Int
angleIndex = Int -> Maybe Int -> Int
forall a. a -> Maybe a -> a
fromMaybe
      (Vector (Float, Float) -> Int
forall a. Unbox a => Vector a -> Int
VU.length Vector (Float, Float)
errorTable) (Maybe Int -> Int) -> Maybe Int -> Int
forall a b. (a -> b) -> a -> b
$
       ((Float, Float) -> Bool) -> Vector (Float, Float) -> Maybe Int
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe Int
VU.findIndex ((Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
< Float
tolerance) (Float -> Bool)
-> ((Float, Float) -> Float) -> (Float, Float) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Float, Float) -> Float
forall a b. (a, b) -> b
snd) Vector (Float, Float)
errorTable

arcSegmentsNeeded :: Float -> Float -> Transformation -> Float
                  -> Int
arcSegmentsNeeded :: Float -> Float -> Transformation -> Float -> Int
arcSegmentsNeeded Float
angle Float
radius Transformation
trans Float
tolerance = Float -> Int
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Float
angle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
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 :: Float
majorAxis = Transformation -> Float -> Float
matrixTransformedCircleMajorAxis Transformation
trans Float
radius
  maxAngle :: Float
maxAngle = Float -> Float
arcMaxAngleForToleranceNormalized (Float
tolerance Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
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 -> Float -> Float
matrixTransformedCircleMajorAxis (Transformation Float
a Float
c Float
_
                                                 Float
b Float
d Float
_) Float
radius =
     Float
radius Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
sqrt (Float
f Float -> Float -> Float
forall a. Num a => a -> a -> a
+ V2 Float -> Float
forall (f :: * -> *) a. (Metric f, Floating a) => f a -> a
norm V2 Float
v)
  where
    i :: Float
i = Float
aFloat -> Float -> Float
forall a. Num a => a -> a -> a
*Float
a Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
bFloat -> Float -> Float
forall a. Num a => a -> a -> a
*Float
b;
    j :: Float
j = Float
cFloat -> Float -> Float
forall a. Num a => a -> a -> a
*Float
c Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
dFloat -> Float -> Float
forall a. Num a => a -> a -> a
*Float
d;

    f :: Float
f = Float
0.5 Float -> Float -> Float
forall a. Num a => a -> a -> a
* (Float
i Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
j)
    v :: V2 Float
v = Float -> Float -> V2 Float
forall a. a -> a -> V2 a
V2 (Float
0.5 Float -> Float -> Float
forall a. Num a => a -> a -> a
* (Float
i Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
j)) (Float
a Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
c Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
b Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
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 :: Float -> Float -> Float
clampAngle Float
angleMin = Float -> Float
go where
  go :: Float -> Float
go Float
angleMax
    | Float
angleMax Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
angleMin Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
> Float
4 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
forall a. Floating a => a
pi = Float -> Float
go (Float -> Float) -> Float -> Float
forall a b. (a -> b) -> a -> b
$ Float
angleMax Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
2 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
forall a. Floating a => a
pi
    | Bool
otherwise = Float
angleMax

subdivideAngles :: (Monoid m)
                => Direction -> (Float -> Float -> m) -> Float -> Float -> m
subdivideAngles :: Direction -> (Float -> Float -> m) -> Float -> Float -> m
subdivideAngles Direction
dir Float -> Float -> m
f Float
aMin = Float -> Float -> m
go Float
aMin (Float -> m) -> (Float -> Float) -> Float -> m
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Float -> Float -> Float
clampAngle Float
aMin where
  go :: Float -> Float -> m
go Float
angleMin Float
angleMax | Float
deltaAngle Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
> Float
forall a. Floating a => a
pi = case Direction
dir of
      Direction
Forward -> Float -> Float -> m
go Float
angleMin Float
angleMid m -> m -> m
forall a. Semigroup a => a -> a -> a
<> Float -> Float -> m
go Float
angleMid Float
angleMax
      Direction
Backward -> Float -> Float -> m
go Float
angleMid Float
angleMax m -> m -> m
forall a. Semigroup a => a -> a -> a
<> Float -> Float -> m
go Float
angleMin Float
angleMid
    where
      deltaAngle :: Float
deltaAngle = Float
angleMax Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
angleMin
      angleMid :: Float
angleMid = Float
angleMin Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
deltaAngle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
2
  go Float
angleMin Float
angleMax = Float -> Float -> m
f Float
angleMin Float
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 Float -> Float -> Float -> Float -> PathCommand
arcSegment (V2 Float
xc Float
yc) Float
radius Float
angleA Float
angleB = V2 Float -> V2 Float -> V2 Float -> PathCommand
PathCubicBezierCurveTo V2 Float
p1 V2 Float
p2 V2 Float
p3 where
  rSinA :: Float
rSinA = Float
radius Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
sin Float
angleA
  rCosA :: Float
rCosA = Float
radius Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
cos Float
angleA
  rSinB :: Float
rSinB = Float
radius Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
sin Float
angleB
  rCosB :: Float
rCosB = Float
radius Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
cos Float
angleB

  h :: Float
h = Float
4.0Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/Float
3.0 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float -> Float
forall a. Floating a => a -> a
tan ((Float
angleB Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
angleA) Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float
4.0)

  p1 :: V2 Float
p1 = Float -> Float -> V2 Float
forall a. a -> a -> V2 a
V2 (Float
xc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
rCosA Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
h Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
rSinA) (Float
yc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
rSinA Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
h Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
rCosA)
  p2 :: V2 Float
p2 = Float -> Float -> V2 Float
forall a. a -> a -> V2 a
V2 (Float
xc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
rCosB Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
h Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
rSinB) (Float
yc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
rSinB Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
h Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
rCosB)
  p3 :: V2 Float
p3 = Float -> Float -> V2 Float
forall a. a -> a -> V2 a
V2 (Float
xc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
rCosB) (Float
yc Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
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 :: V2 Float
-> Float -> Direction -> Float -> Float -> Float -> [PathCommand]
arcInDirection p :: V2 Float
p@(V2 Float
px Float
py) Float
radius Direction
dir Float
tolerance 
  | Float -> Bool
forall a. RealFloat a => a -> Bool
isNaN Float
px Bool -> Bool -> Bool
|| Float -> Bool
forall a. RealFloat a => a -> Bool
isNaN Float
py Bool -> Bool -> Bool
|| Float -> Bool
forall a. RealFloat a => a -> Bool
isNaN Float
radius = Float -> Float -> [PathCommand]
forall a. Monoid a => a
mempty
  | Bool
otherwise = Direction
-> (Float -> Float -> [PathCommand])
-> Float
-> Float
-> [PathCommand]
forall m.
Monoid m =>
Direction -> (Float -> Float -> m) -> Float -> Float -> m
subdivideAngles Direction
dir Float -> Float -> [PathCommand]
go where
  go :: Float -> Float -> [PathCommand]
go Float
angleMin Float
angleMax = [PathCommand]
commands where
    deltaAngle :: Float
deltaAngle = Float
angleMax Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
angleMin
    segmentCount :: Int
segmentCount = Float -> Float -> Transformation -> Float -> Int
arcSegmentsNeeded Float
deltaAngle Float
radius Transformation
forall a. Monoid a => a
mempty Float
tolerance

    (Float
angle, Float
angleStep) = case Direction
dir of
      Direction
Forward -> (Float
angleMin, Float
deltaAngle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
segmentCount)
      Direction
Backward -> (Float
angleMax, - Float
deltaAngle Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
segmentCount)

    commands :: [PathCommand]
commands =
        [V2 Float -> Float -> Float -> Float -> PathCommand
arcSegment V2 Float
p Float
radius Float
a (Float
a Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
angleStep)
            | Int
i <- [Int
0 .. Int
segmentCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            , let a :: Float
a = Float
angle Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
angleStep Float -> Float -> Float
forall a. Num a => a -> a -> a
* Int -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i]