module Graphics.Rasterific.Arc( Direction( .. ), arcInDirection ) where
import Data.Maybe( fromMaybe )
import Data.Monoid( (<>) )
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 ]
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
majorAxis = matrixTransformedCircleMajorAxis trans radius
maxAngle = arcMaxAngleForToleranceNormalized (tolerance / majorAxis)
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)
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
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)
arcInDirection :: Point
-> Float
-> Direction
-> Float
-> Float
-> Float
-> [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]