{-# LANGUAGE MultiParamTypeClasses, TypeFamilies, FlexibleContexts, FlexibleInstances, DeriveGeneric #-}
module Graphics.Rendering.Plot.Light.Internal.Geometry
(
Point(..), mkPoint, setPointX, setPointY,
LabeledPoint(..), mkLabeledPoint, labelPoint, mapLabel,
Frame(..), mkFrame, unitFrame, frameFromPoints, mkFrameOrigin, height, width, xmin, xmax, ymin, ymax, isPointInFrame,
Axis(..), otherAxis,
V2(..), pointFromV2,
Mat2(..), DiagMat2(..), diagMat2,
origin, oneOne, e1, e2,
norm2, normalize2,
v2fromEndpoints, v2fromPoint,
movePoint, moveLabeledPoint, moveLabeledPointV2, moveLabeledPointBwFrames, (-.), pointRange,
frameToFrame, frameToFrameValue,
AdditiveGroup(..), VectorSpace(..), Hermitian(..), LinearMap(..), MultiplicativeSemigroup(..), MatrixGroup(..), Eps(..),
meshGrid, subdivSegment, interpolateBilinear
)
where
import Control.Exception
import Control.Monad.Catch (MonadThrow(..), throwM)
import GHC.Generics
import Data.Semigroup (Semigroup(..))
data Point a = Point { _px :: a,
_py :: a } deriving (Eq, Generic)
instance Ord a => Ord (Point a) where
(Point x1 y1) <= (Point x2 y2) = x1 <= x2 && y1 <= y2
instance Show a => Show (Point a) where
show (Point x y) = show x ++ "," ++ show y
mkPoint :: a -> a -> Point a
mkPoint = Point
lift2Point :: (a -> b -> c) -> Point a -> Point b -> Point c
lift2Point f (Point a b) (Point c d) = Point (f a c) (f b d)
lift1Point :: (a -> a -> b) -> Point a -> b
lift1Point f (Point x y) = f x y
pointInf, pointSup :: (Ord a) => Point a -> Point a -> Point a
pointInf = lift2Point min
pointSup = lift2Point max
origin :: Num a => Point a
origin = Point 0 0
oneOne :: Num a => Point a
oneOne = Point 1 1
norm2fromOrigin :: Floating a => Point a -> a
norm2fromOrigin p = norm2 $ p -. origin
setPointCoord :: Axis -> a -> Point a -> Point a
setPointCoord axis c (Point x y)
| axis == X = Point c y
| otherwise = Point x c
setPointX, setPointY :: a -> Point a -> Point a
setPointX = setPointCoord X
setPointY = setPointCoord Y
data LabeledPoint l a =
LabeledPoint {
_lp :: Point a,
_lplabel :: l
} deriving (Eq, Show)
mkLabeledPoint :: Point a -> l -> LabeledPoint l a
mkLabeledPoint = LabeledPoint
labelPoint :: (Point a -> l) -> Point a -> LabeledPoint l a
labelPoint lf p = LabeledPoint p (lf p)
moveLabeledPoint :: (Point a -> Point b) -> LabeledPoint l a -> LabeledPoint l b
moveLabeledPoint f (LabeledPoint p l) = LabeledPoint (f p) l
mapLabel :: (l1 -> l2) -> LabeledPoint l1 a -> LabeledPoint l2 a
mapLabel f (LabeledPoint p l) = LabeledPoint p (f l)
data Frame a = Frame {
_fpmin :: Point a,
_fpmax :: Point a
} deriving (Eq, Show, Generic)
instance (Ord a) => Semigroup (Frame a) where
(Frame p1min p1max) <> (Frame p2min p2max) = Frame (pointInf p1min p2min) (pointSup p1max p2max)
instance (Ord a, Num a) => Monoid (Frame a) where
mempty = Frame (Point 0 0) (Point 0 0)
mappend = (<>)
isPointInFrame :: Ord a => Frame a -> Point a -> Bool
isPointInFrame (Frame p1 p2) p = p >= p1 && p <= p2
mkFrame :: Point a -> Point a -> Frame a
mkFrame = Frame
mkFrameOrigin :: Num a => a -> a -> Frame a
mkFrameOrigin w h = Frame origin (Point w h)
unitFrame :: Num a => Frame a
unitFrame = mkFrame origin oneOne
frameFromPoints :: (Ord a, Foldable t, Functor t) =>
t (Point a) -> Frame a
frameFromPoints ds = mkFrame (Point mx my) (Point mmx mmy)
where
xcoord = _px <$> ds
ycoord = _py <$> ds
mmx = maximum xcoord
mmy = maximum ycoord
mx = minimum xcoord
my = minimum ycoord
xmin, xmax, ymin, ymax :: Frame a -> a
xmin = _px . _fpmin
xmax = _px . _fpmax
ymin = _py . _fpmin
ymax = _py . _fpmax
width, height :: Num a => Frame a -> a
width f = abs $ xmax f - xmin f
height f = abs $ ymax f - ymin f
interpolateBilinear :: (Ord p, Fractional p, Show p) =>
Frame p -> (Point p -> p) -> Point p -> p
interpolateBilinear fr@(Frame p1 p2) f p
| isPointInFrame fr p = interpolateBilinear' p1 p2 f p
| otherwise = error $ unwords ["Point", show p, "is outside frame", show fr]
interpolateBilinear' :: Fractional a => Point a -> Point a -> (Point a -> a) -> Point a -> a
interpolateBilinear' q11@(Point x1 y1) q22@(Point x2 y2) f (Point x y) =
let
q12 = Point x1 y2
q21 = Point x2 y1
fq11 = f q11
fq22 = f q22
fq12 = f q12
fq21 = f q21
den1 = (x1 - x2) * (y1 - y2)
den2 = (x1 - x2) * (y2 - y1)
c111 = fq11/den1
c112 = fq11/den2
c121 = fq12/den1
c122 = fq12/den2
c211 = fq21/den1
c212 = fq21/den2
c221 = fq22/den1
c222 = fq22/den2
a0 = c111 * x2 * y2 + c122 * x2 * y1 + c212 * x1 * y2 + c221 * x1 * y1
a1 = c112 * y2 + c121 * y1 + c211 * y2 + c222 * y1
a2 = c112 * x2 + c121 * x2 + c211 * x1 + c222 * x1
a3 = c111 + c122 + c212 + c221
in a0 + a1 * x + a2 * y + a3 * x * y
data Axis = X | Y deriving (Eq, Show)
otherAxis :: Axis -> Axis
otherAxis X = Y
otherAxis _ = X
data V2 a = V2 a a deriving (Eq, Show)
instance Num a => Semigroup (V2 a) where
(V2 a b) <> (V2 c d) = V2 (a + c) (b + d)
instance Num a => Monoid (V2 a) where
mempty = V2 0 0
mappend = (<>)
class AdditiveGroup v where
zero :: v
(^+^) :: v -> v -> v
(^-^) :: v -> v -> v
instance Num a => AdditiveGroup (V2 a) where
zero = mempty
(^+^) = mappend
(V2 a b) ^-^ (V2 c d) = V2 (a - c) (b - d)
class AdditiveGroup v => VectorSpace v where
type Scalar v :: *
(.*) :: Scalar v -> v -> v
instance Num a => VectorSpace (V2 a) where
type Scalar (V2 a) = a
n .* (V2 vx vy) = V2 (n*vx) (n*vy)
class VectorSpace v => Hermitian v where
type InnerProduct v :: *
(<.>) :: v -> v -> InnerProduct v
instance Num a => Hermitian (V2 a) where
type InnerProduct (V2 a) = a
(V2 a b) <.> (V2 c d) = (a*c) + (b*d)
norm2 ::
(Hermitian v, Floating n, n ~ (InnerProduct v)) => v -> n
norm2 v = sqrt $ v <.> v
normalize2 :: (InnerProduct v ~ Scalar v, Floating (Scalar v), Hermitian v) =>
v -> v
normalize2 v = (1/norm2 v) .* v
v2fromEndpoints, (-.) :: Num a => Point a -> Point a -> V2 a
v2fromEndpoints (Point px py) (Point qx qy) = V2 (qx-px) (qy-py)
(-.) = v2fromEndpoints
data Mat2 a = Mat2 a a a a deriving (Eq, Show)
class Hermitian v => LinearMap m v where
(#>) :: m -> v -> v
class MultiplicativeSemigroup m where
(##) :: m -> m -> m
instance Num a => MultiplicativeSemigroup (Mat2 a) where
Mat2 a00 a01 a10 a11 ## Mat2 b00 b01 b10 b11 = Mat2 (a00*b00+a01*b10) (a00*b01+a01*b11) (a10*b00+a11*b10) (a10*b01+a11*b11)
instance Num a => LinearMap (Mat2 a) (V2 a) where
(Mat2 a00 a01 a10 a11) #> (V2 vx vy) = V2 (a00 * vx + a01 * vy) (a10 * vx + a11 * vy)
data DiagMat2 a = DMat2 a a deriving (Eq, Show)
instance Num a => Semigroup (DiagMat2 a) where
(<>) = (##)
instance Num a => Monoid (DiagMat2 a) where
mempty = DMat2 1 1
mappend = (<>)
instance Num a => Semigroup (Mat2 a) where
(<>) = (##)
instance Num a => Monoid (Mat2 a) where
mempty = Mat2 1 0 0 1
mappend = (<>)
diagMat2 :: Num a => a -> a -> DiagMat2 a
diagMat2 = DMat2
rotMtx :: Floating a => a -> Mat2 a
rotMtx r = Mat2 (cos r) (- (sin r)) (sin r) (cos r)
class LinearMap m v => MatrixGroup m v where
(<\>) :: m -> v -> v
instance Num a => MultiplicativeSemigroup (DiagMat2 a) where
DMat2 a b ## DMat2 c d = DMat2 (a*c) (b*d)
instance Num a => LinearMap (DiagMat2 a) (V2 a) where
DMat2 d1 d2 #> V2 vx vy = V2 (d1 * vx) (d2 * vy)
instance Fractional a => MatrixGroup (DiagMat2 a) (V2 a) where
DMat2 d1 d2 <\> V2 vx vy = V2 (vx / d1) (vy / d2)
v2fromPoint :: Num a => Point a -> V2 a
v2fromPoint p = origin -. p
pointFromV2 :: V2 a -> Point a
pointFromV2 (V2 x y) = Point x y
movePoint :: Num a => V2 a -> Point a -> Point a
movePoint (V2 vx vy) (Point px py) = Point (px + vx) (py + vy)
moveLabeledPointV2 :: Num a => V2 a -> LabeledPoint l a -> LabeledPoint l a
moveLabeledPointV2 = moveLabeledPoint . movePoint
pointRange :: (Fractional a, Integral n) =>
n -> Point a -> Point a -> [Point a]
pointRange n p q = [ movePoint (fromIntegral x .* vnth) p | x <- [0 .. n]]
where
v = p -. q
vnth = (1/fromIntegral n) .* v
meshGrid
:: (Enum a, RealFrac a) =>
Frame a
-> Int
-> Int
-> [Point a]
meshGrid (Frame (Point xmi ymi) (Point xma yma)) nx ny =
[Point x y |
x <- take nx $ subdivSegment xmi xma nx,
y <- take ny $ subdivSegment ymi yma ny]
data MeshGrid a = MeshGrid (Frame a) Int Int deriving (Eq, Show, Generic)
subdivSegment :: (Real a, Enum b, RealFrac b) => a -> a -> Int -> [b]
subdivSegment x1 x2 n = f <$> [0, 1 ..] where
xmin = min x1 x2
xmax = max x1 x2
l = fromRational $ toRational (xmax - xmin)
f x = x * l/fromIntegral n + fromRational (toRational xmin)
fromFrame :: Fractional a => Frame a -> V2 a -> V2 a
fromFrame from v = mfrom <\> (v ^-^ vfrom) where
vfrom = v2fromPoint (_fpmin from)
mfrom = diagMat2 (width from) (height from)
toFrame :: Num a => Frame a -> V2 a -> V2 a
toFrame to v01 = (mto #> v01) ^+^ vto where
vto = v2fromPoint (_fpmin to)
mto = diagMat2 (width to) (height to)
frameToFrame :: Fractional a =>
Frame a
-> Frame a
-> Bool
-> Bool
-> V2 a
-> V2 a
frameToFrame from to fliplr flipud v = toFrame to v01'
where
v01 = fromFrame from v
v01' | fliplr && flipud = flipLR01 (flipUD01 v01)
| fliplr = flipLR01 v01
| flipud = flipUD01 v01
| otherwise = v01
flipLR01, flipUD01 :: Num a => V2 a -> V2 a
flipLR01 (V2 a b) = V2 (1 - a) b
flipUD01 (V2 a b) = V2 a (1 - b)
frameToFrameValue :: Fractional t =>
Frame t
-> Frame t
-> t
-> t
frameToFrameValue from to x = (x01 * rto) + ymin to where
x01 = (x - ymin from)/rfrom
rfrom = height from
rto = height to
moveLabeledPointBwFrames ::
Fractional a =>
Frame a
-> Frame a
-> Bool
-> Bool
-> LabeledPoint l a
-> LabeledPoint l a
moveLabeledPointBwFrames from to fliplr flipud lp = LabeledPoint p' (_lplabel lp)
where
vlp = v2fromPoint $ _lp lp
vlp' = frameToFrame from to fliplr flipud vlp
p' = pointFromV2 vlp'
e1 :: Num a => V2 a
e1 = V2 1 0
e2 :: Num a => V2 a
e2 = V2 0 1
class Eps a where
(~=) :: a -> a -> Bool
instance Eps Double where
a ~= b = abs (a - b) <= 1e-12
instance Eps Float where
a ~= b = abs (a - b) <= 1e-6
instance Eps (V2 Double) where
v1 ~= v2 = norm2 (v1 ^-^ v2) <= 1e-8
instance Eps (V2 Float) where
v1 ~= v2 = norm2 (v1 ^-^ v2) <= 1e-2