module Geom2D.CubicBezier.Intersection
(bezierIntersection, bezierLineIntersections, bezierFindRoot)
where
import Geom2D
import Geom2D.CubicBezier.Basic
import Math.BernsteinPoly
import Data.Maybe
findOuter' upper !dir !p1 l@(p2:rest)
| if upper
then dir `vectorCross` (p2^-^p1) > 0
else dir `vectorCross` (p2^-^p1) < 0 = Left l
| otherwise = case findOuter' upper (p2^-^p1) p2 rest of
Left m -> findOuter' upper dir p1 m
Right m -> Right (p1:m)
findOuter' _ _ p1 p = Right (p1:p)
findOuter upper (p1:p2:rest) =
case findOuter' upper (p2^-^p1) p2 rest of
Right l -> p1:l
Left l -> findOuter upper (p1:l)
findOuter _ l = l
makeHull :: [Double] -> ([Point], [Point])
makeHull ds =
let n = fromIntegral $ length ds 1
points = zipWith Point [i/n | i <- [0..n]] ds
in (findOuter True points,
findOuter False points)
testBelow :: Double -> [Point] -> Maybe Double -> Maybe Double
testBelow dmin [] _ = Nothing
testBelow dmin [_] _ = Nothing
testBelow dmin (p:q:rest) cont
| pointY p >= dmin = cont
| pointY p > pointY q = Nothing
| pointY q < dmin = testBelow dmin (q:rest) cont
| otherwise = Just $ intersectPt dmin p q
testBetween :: Double -> Point -> Maybe Double -> Maybe Double
testBetween dmax (Point x y) cont
| y <= dmax = Just x
| otherwise = cont
testAbove :: Double -> [Point] -> Maybe Double
testAbove dmax [] = Nothing
testAbove dmax [_] = Nothing
testAbove dmax (p:q:rest)
| pointY p < pointY q = Nothing
| pointY q > dmax = testAbove dmax (q:rest)
| otherwise = Just $ intersectPt dmax p q
intersectPt d (Point x1 y1) (Point x2 y2) =
x1 + (d y1) * (x2 x1) / (y2 y1)
chopHull dmin dmax ds = do
let (upper, lower) = makeHull ds
left_t <- testBelow dmin upper $
testBetween dmax (head upper) $
testAbove dmax lower
right_t <- testBelow dmin (reverse upper) $
testBetween dmax (last upper) $
testAbove dmax (reverse lower)
Just (left_t, right_t)
bezierClip p@(CubicBezier !p0 !p1 !p2 !p3) q@(CubicBezier !q0 !q1 !q2 !q3)
tmin tmax umin umax prevClip eps reverse
| isNothing chop_interval = []
| prevClip > 0.8 && newClip > 0.8 =
if new_tmax new_tmin > umax umin
then let
(pl, pr) = splitBezier newP 0.5
half_t = new_tmin + (new_tmax new_tmin) / 2
in bezierClip q pl umin umax new_tmin half_t newClip eps (not reverse) ++
bezierClip q pr umin umax half_t new_tmax newClip eps (not reverse)
else let
(ql, qr) = splitBezier q 0.5
half_t = umin + (umax umin) / 2
in bezierClip ql newP umin half_t new_tmin new_tmax newClip eps (not reverse) ++
bezierClip qr newP half_t umax new_tmin new_tmax newClip eps (not reverse)
| max (umax umin) (new_tmax new_tmin) < eps =
if reverse
then [ (umin + (umaxumin)/2,
new_tmin + (new_tmaxnew_tmin)/2) ]
else [ (new_tmin + (new_tmaxnew_tmin)/2,
umin + (umaxumin)/2) ]
| otherwise =
bezierClip q newP umin umax new_tmin new_tmax newClip eps (not reverse)
where
d = lineDistance (Line q0 q3)
d1 = d q1
d2 = d q2
(dmin, dmax) | d1*d2 > 0 = (3/4 * minimum [0, d1, d2],
3/4 * maximum [0, d1, d2])
| otherwise = (4/9 * minimum [0, d1, d2],
4/9 * maximum [0, d1, d2])
chop_interval = chopHull dmin dmax $
map d [p0, p1, p2, p3]
Just (chop_tmin, chop_tmax) = chop_interval
newP = bezierSubsegment p chop_tmin chop_tmax
newClip = chop_tmax chop_tmin
new_tmin = tmax * chop_tmin + tmin * (1 chop_tmin)
new_tmax = tmax * chop_tmax + tmin * (1 chop_tmax)
bezierIntersection :: CubicBezier -> CubicBezier -> Double -> [(Double, Double)]
bezierIntersection p q eps = bezierClip p q 0 1 0 1 0 eps' False
where
eps' = min (bezierParamTolerance p eps) (bezierParamTolerance q eps)
bezierFindRoot :: BernsteinPoly
-> Double
-> Double
-> Double
-> [Double]
bezierFindRoot p tmin tmax eps
| chop_interval == Nothing = []
| clip > 0.8 =
let (p1, p2) = bernsteinSplit newP 0.5
half_t = new_tmin + (new_tmax new_tmin) / 2
in bezierFindRoot p1 new_tmin half_t eps ++
bezierFindRoot p2 half_t new_tmax eps
| new_tmax new_tmin < eps =
[new_tmin + (new_tmaxnew_tmin)/2]
| otherwise =
bezierFindRoot newP new_tmin new_tmax eps
where
chop_interval = chopHull 0 0 (bernsteinCoeffs p)
Just (chop_tmin, chop_tmax) = chop_interval
newP = bernsteinSubsegment p chop_tmin chop_tmax
clip = chop_tmax chop_tmin
new_tmin = tmax * chop_tmin + tmin * (1 chop_tmin)
new_tmax = tmax * chop_tmax + tmin * (1 chop_tmax)
bezierLineIntersections :: CubicBezier -> Line -> Double -> [Double]
bezierLineIntersections b (Line p q) eps =
bezierFindRoot (listToBernstein $ map pointY [p0, p1, p2, p3]) 0 1 $
bezierParamTolerance b eps
where (CubicBezier p0 p1 p2 p3) =
fromJust (inverse $ translate p $* rotateVec (q ^-^ p)) $* b