module Wumpus.Basic.Geometry.Base
(
quarter_pi
, half_pi
, two_pi
, Matrix2'2(..)
, DMatrix2'2
, identity2'2
, det2'2
, transpose2'2
, Line(..)
, inclinedLine
, LineEquation(..)
, DLineEquation
, lineEquation
, pointViaX
, pointViaY
, pointLineDistance
, LineSegment(..)
, DLineSegment
, rectangleLineSegments
, polygonLineSegments
, BezierCurve(..)
, DBezierCurve
, bezierLength
, subdivide
, subdividet
, bezierArcPoints
, bezierMinorArc
, affineComb
, midpoint
, lineAngle
)
where
import Wumpus.Basic.Kernel
import Wumpus.Core
import Data.AffineSpace
import Data.VectorSpace
quarter_pi :: Floating u => u
quarter_pi = 0.25 * pi
half_pi :: Floating u => u
half_pi = 0.5 * pi
two_pi :: Floating u => u
two_pi = 2.0 * pi
data Matrix2'2 u = M2'2 !u !u !u !u
deriving (Eq)
type instance DUnit (Matrix2'2 u) = u
type DMatrix2'2 = Matrix2'2 Double
instance Functor Matrix2'2 where
fmap f (M2'2 a b c d) = M2'2 (f a) (f b) (f c) (f d)
instance Show u => Show (Matrix2'2 u) where
show (M2'2 a b c d) = "(M2'2 " ++ body ++ ")"
where
body = show [[a,b],[c,d]]
lift2Matrix2'2 :: (u -> u -> u) -> Matrix2'2 u -> Matrix2'2 u -> Matrix2'2 u
lift2Matrix2'2 op (M2'2 a b c d) (M2'2 m n o p) =
M2'2 (a `op` m) (b `op` n)
(c `op` o) (d `op` p)
instance Num u => Num (Matrix2'2 u) where
(+) = lift2Matrix2'2 (+)
() = lift2Matrix2'2 ()
(*) (M2'2 a b c d) (M2'2 m n o p) =
M2'2 (a*m + b*o) (a*n + b*p)
(c*m + d*o) (c*n + d*p)
abs = fmap abs
negate = fmap negate
signum = fmap signum
fromInteger a = M2'2 a' a' a' a' where a' = fromInteger a
identity2'2 :: Num u => Matrix2'2 u
identity2'2 = M2'2 1 0
0 1
det2'2 :: Num u => Matrix2'2 u -> u
det2'2 (M2'2 a b c d) = a*d b*c
transpose2'2 :: Matrix2'2 u -> Matrix2'2 u
transpose2'2 (M2'2 a b
c d) = M2'2 a c
b d
data Line u = Line (Point2 u) (Point2 u)
deriving (Eq,Show)
type instance DUnit (Line u) = u
inclinedLine :: Floating u => Point2 u -> Radian -> Line u
inclinedLine radial_ogin ang = Line radial_ogin (radial_ogin .+^ avec ang 100)
data LineEquation u = LineEquation
{ line_eqn_A :: !u
, line_eqn_B :: !u
, line_eqn_C :: !u
}
deriving (Eq,Show)
type instance DUnit (LineEquation u) = u
type DLineEquation = LineEquation Double
lineEquation :: Num u => Point2 u -> Point2 u -> LineEquation u
lineEquation (P2 x1 y1) (P2 x2 y2) = LineEquation a b c
where
a = y1 y2
b = x2 x1
c = (x1*y2) (x2*y1)
pointViaX :: Fractional u => u -> LineEquation u -> Point2 u
pointViaX x (LineEquation a b c) = P2 x y
where
y = ((a*x) + c) / (b)
pointViaY :: Fractional u => u -> LineEquation u -> Point2 u
pointViaY y (LineEquation a b c) = P2 x y
where
x = ((b*y) + c) / (a)
pointLineDistance :: Floating u => Point2 u -> LineEquation u -> u
pointLineDistance (P2 u v) (LineEquation a b c) =
((a*u) + (b*v) + c) / base
where
base = sqrt $ (a^two) + (b^two)
two :: Integer
two = 2
data LineSegment u = LineSegment (Point2 u) (Point2 u)
deriving (Eq,Ord,Show)
type instance DUnit (LineSegment u) = u
type DLineSegment = LineSegment Double
rectangleLineSegments :: Num u => u -> u -> Point2 u -> [LineSegment u]
rectangleLineSegments hw hh ctr =
[ LineSegment br tr, LineSegment tr tl, LineSegment tl bl
, LineSegment bl br
]
where
br = ctr .+^ vec hw (hh)
tr = ctr .+^ vec hw hh
tl = ctr .+^ vec (hw) hh
bl = ctr .+^ vec (hw) (hh)
polygonLineSegments :: [Point2 u] -> [LineSegment u]
polygonLineSegments [] = []
polygonLineSegments (x:xs) = step x xs
where
step a [] = [LineSegment a x]
step a (b:bs) = (LineSegment a b) : step b bs
data BezierCurve u = BezierCurve !(Point2 u) !(Point2 u) !(Point2 u) !(Point2 u)
deriving (Eq,Ord,Show)
type instance DUnit (BezierCurve u) = u
type DBezierCurve = BezierCurve Double
bezierLength :: (Floating u, Ord u, Tolerance u)
=> BezierCurve u -> u
bezierLength = gravesenLength length_tolerance
gravesenLength :: (Floating u, Ord u) => u -> BezierCurve u -> u
gravesenLength err_tol crv = step crv
where
step c = let l1 = ctrlPolyLength c
l0 = cordLength c
in if l1l0 > err_tol
then let (a,b) = subdivide c in step a + step b
else 0.5*l0 + 0.5*l1
ctrlPolyLength :: Floating u => BezierCurve u -> u
ctrlPolyLength (BezierCurve p0 p1 p2 p3) = len p0 p1 + len p1 p2 + len p2 p3
where
len pa pb = vlength $ pvec pa pb
cordLength :: Floating u => BezierCurve u -> u
cordLength (BezierCurve p0 _ _ p3) = vlength $ pvec p0 p3
subdivide :: Fractional u
=> BezierCurve u -> (BezierCurve u, BezierCurve u)
subdivide (BezierCurve p0 p1 p2 p3) =
(BezierCurve p0 p01 p012 p0123, BezierCurve p0123 p123 p23 p3)
where
p01 = midpoint p0 p1
p12 = midpoint p1 p2
p23 = midpoint p2 p3
p012 = midpoint p01 p12
p123 = midpoint p12 p23
p0123 = midpoint p012 p123
subdividet :: Real u
=> u -> BezierCurve u -> (BezierCurve u, BezierCurve u)
subdividet t (BezierCurve p0 p1 p2 p3) =
(BezierCurve p0 p01 p012 p0123, BezierCurve p0123 p123 p23 p3)
where
p01 = affineComb t p0 p1
p12 = affineComb t p1 p2
p23 = affineComb t p2 p3
p012 = affineComb t p01 p12
p123 = affineComb t p12 p23
p0123 = affineComb t p012 p123
kappa :: Floating u => u
kappa = 4 * ((sqrt 2 1) / 3)
bezierArcPoints :: Floating u
=> Radian -> u -> Radian -> Point2 u -> [Point2 u]
bezierArcPoints ang radius theta pt = go (circularModulo ang)
where
go a | a <= half_pi = wedge1 a
| a <= pi = wedge2 (a/2)
| a <= 1.5*pi = wedge3 (a/3)
| otherwise = wedge4 (a/4)
wedge1 a =
let (BezierCurve p0 p1 p2 p3) = bezierMinorArc a radius theta pt
in [p0,p1,p2,p3]
wedge2 a =
let (BezierCurve p0 p1 p2 p3) = bezierMinorArc a radius theta pt
(BezierCurve _ p4 p5 p6) = bezierMinorArc a radius (theta+a) pt
in [ p0,p1,p2,p3, p4,p5,p6 ]
wedge3 a =
let (BezierCurve p0 p1 p2 p3) = bezierMinorArc a radius theta pt
(BezierCurve _ p4 p5 p6) = bezierMinorArc a radius (theta+a) pt
(BezierCurve _ p7 p8 p9) = bezierMinorArc a radius (theta+a+a) pt
in [ p0,p1,p2,p3, p4,p5,p6, p7, p8, p9 ]
wedge4 a =
let (BezierCurve p0 p1 p2 p3) = bezierMinorArc a radius theta pt
(BezierCurve _ p4 p5 p6) = bezierMinorArc a radius (theta+a) pt
(BezierCurve _ p7 p8 p9) = bezierMinorArc a radius (theta+a+a) pt
(BezierCurve _ p10 p11 p12) = bezierMinorArc a radius (theta+a+a+a) pt
in [ p0,p1,p2,p3, p4,p5,p6, p7,p8,p9, p10,p11, p12 ]
bezierMinorArc :: Floating u
=> Radian -> u -> Radian -> Point2 u -> BezierCurve u
bezierMinorArc ang radius theta pt = BezierCurve p0 c1 c2 p3
where
kfactor = fromRadian $ ang / (0.5*pi)
rl = kfactor * radius * kappa
totang = circularModulo $ ang + theta
p0 = dispParallel radius theta pt
c1 = dispPerpendicular rl theta p0
c2 = dispPerpendicular (rl) totang p3
p3 = dispParallel radius totang pt
affineComb :: Real u => u -> Point2 u -> Point2 u -> Point2 u
affineComb t p1 p2 = p1 .+^ t *^ (p2 .-. p1)
midpoint :: Fractional u => Point2 u -> Point2 u -> Point2 u
midpoint p0 p1 = p0 .+^ v1 ^/ 2 where v1 = p1 .-. p0
lineAngle :: (Floating u, Real u) => Point2 u -> Point2 u -> Radian
lineAngle (P2 x1 y1) (P2 x2 y2) = step (x2 x1) (y2 y1)
where
step x y | pve x && pve y = toRadian $ atan (y/x)
step x y | pve y = pi (toRadian $ atan (y / abs x))
step x y | pve x = (2*pi) (toRadian $ atan (abs y / x))
step x y = pi + (toRadian $ atan (y/x))
pve a = signum a >= 0