hgeometry-0.8.0.0: src/Data/Geometry/Triangle.hs
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE UndecidableInstances #-}
module Data.Geometry.Triangle where
import Control.Lens
import Data.Bifunctor
import Data.Either (partitionEithers)
import Data.Ext
import Data.Geometry.Ball (Disk, disk)
import Data.Geometry.Boundary
import Data.Geometry.HyperPlane
import Data.Geometry.Line
import Data.Geometry.LineSegment
import Data.Geometry.Point
import Data.Geometry.Properties
import Data.Geometry.Transformation
import Data.Geometry.Vector
import qualified Data.Geometry.Vector as V
import qualified Data.List as List
import Data.Maybe (mapMaybe)
import Data.Vinyl
import Data.Vinyl.CoRec
import GHC.TypeLits
--------------------------------------------------------------------------------
data Triangle d p r = Triangle (Point d r :+ p)
(Point d r :+ p)
(Point d r :+ p)
deriving instance (Arity d, Show r, Show p) => Show (Triangle d p r)
deriving instance (Arity d, Read r, Read p) => Read (Triangle d p r)
deriving instance (Arity d, Eq r, Eq p) => Eq (Triangle d p r)
instance Arity d => Functor (Triangle d p) where
fmap f (Triangle p q r) = let f' = first (fmap f) in Triangle (f' p) (f' q) (f' r)
type instance NumType (Triangle d p r) = r
type instance Dimension (Triangle d p r) = d
instance PointFunctor (Triangle d p) where
pmap f (Triangle p q r) = Triangle (p&core %~ f) (q&core %~ f) (r&core %~ f)
instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (Triangle d p r) where
transformBy = transformPointFunctor
-- | convenience function to construct a triangle without associated data.
triangle' :: Point d r -> Point d r -> Point d r -> Triangle d () r
triangle' p q r = Triangle (ext p) (ext q) (ext r)
sideSegments :: Triangle d p r -> [LineSegment d p r]
sideSegments (Triangle p q r) =
[ClosedLineSegment p q, ClosedLineSegment q r, ClosedLineSegment r p]
-- | Compute the area of a triangle
area :: Fractional r => Triangle 2 p r -> r
area t = doubleArea t / 2
-- | 2*the area of a triangle.
doubleArea :: Num r => Triangle 2 p r -> r
doubleArea (Triangle a b c) = abs $ ax*by - ax*cy
+ bx*cy - bx*ay
+ cx*ay - cx*by
-- Based on determinant of a 3x3 matrix (shoelace formula)
where
Point2 ax ay = a^.core
Point2 bx by = b^.core
Point2 cx cy = c^.core
-- | Checks if the triangle is degenerate, i.e. has zero area.
isDegenerateTriangle :: (Num r, Eq r) => Triangle 2 p r -> Bool
isDegenerateTriangle = (== 0) . doubleArea
-- | get the inscribed disk. Returns Nothing if the triangle is degenerate,
-- i.e. if the points are colinear.
inscribedDisk :: (Eq r, Fractional r)
=> Triangle 2 p r -> Maybe (Disk () r)
inscribedDisk (Triangle p q r) = disk (p^.core) (q^.core) (r^.core)
instance Num r => HasSupportingPlane (Triangle 3 p r) where
supportingPlane (Triangle p q r) = from3Points (p^.core) (q^.core) (r^.core)
-- | Given a point q and a triangle, q inside the triangle, get the baricentric
-- cordinates of q
toBarricentric :: Fractional r
=> Point 2 r -> Triangle 2 p r
-> Vector 3 r
toBarricentric (Point2 qx qy) (Triangle a b c) = Vector3 alpha beta gamma
where
Point2 ax ay = a^.core
Point2 bx by = b^.core
Point2 cx cy = c^.core
dett = (by - cy)*(ax - cx) + (cx - bx)*(ay - cy)
alpha = ((by - cy)*(qx - cx) + (cx - bx)*(qy - cy)) / dett
beta = ((cy - ay)*(qx - cx) + (ax - cx)*(qy - cy)) / dett
gamma = 1 - alpha - beta
-- see https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Conversion_between_barycentric_and_Cartesian_coordinates
-- | Given a vector of barricentric coordinates and a triangle, get the
-- corresponding point in the same coordinate sytsem as the vertices of the
-- triangle.
fromBarricentric :: (Arity d, Num r)
=> Vector 3 r -> Triangle d p r
-> Point d r
fromBarricentric (Vector3 a b c) (Triangle p q r) = let f = view (core.vector) in
Point $ a *^ f p ^+^ b *^ f q ^+^ c *^ f r
-- | Tests if a point lies inside a triangle, on its boundary, or outside the triangle
inTriangle :: (Ord r, Fractional r)
=> Point 2 r -> Triangle 2 p r -> PointLocationResult
inTriangle q t
| all (`inRange` (OpenRange 0 1)) [a,b,c] = Inside
| all (`inRange` (ClosedRange 0 1)) [a,b,c] = OnBoundary
| otherwise = Outside
where
Vector3 a b c = toBarricentric q t
-- | Test if a point lies inside or on the boundary of a triangle
onTriangle :: (Ord r, Fractional r)
=> Point 2 r -> Triangle 2 p r -> Bool
q `onTriangle` t = let Vector3 a b c = toBarricentric q t
in all (`inRange` (ClosedRange 0 1)) [a,b,c]
-- myQ :: Point 2 Rational
-- myQ = read "Point2 [(-5985) % 16,(-14625) % 1]"
-- myTri :: Triangle 2 () Rational
-- myTri = read "Triangle (Point2 [(-15) % 1,0 % 1] :+ ()) (Point2 [225 % 2,0 % 1] :+ ()) (Point2 [135 % 1,0 % 1] :+ ())"
type instance IntersectionOf (Line 2 r) (Triangle 2 p r) =
[ NoIntersection, Point 2 r, LineSegment 2 () r ]
instance (Fractional r, Ord r) => (Line 2 r) `IsIntersectableWith` (Triangle 2 p r) where
nonEmptyIntersection = defaultNonEmptyIntersection
l `intersect` (Triangle p q r) =
case first List.nub . partitionEithers . mapMaybe collect $ sides of
([],[]) -> coRec NoIntersection
(_, [s]) -> coRec $ first (const ()) s
([a],_) -> coRec a
([a,b],_) -> coRec $ ClosedLineSegment (ext a) (ext b)
(_,_) -> error "intersecting a line with a triangle. Triangle is degenerate"
where
sides = [ClosedLineSegment p q, ClosedLineSegment q r, ClosedLineSegment r p]
collect :: LineSegment 2 p r -> Maybe (Either (Point 2 r) (LineSegment 2 p r))
collect s = match (s `intersect` l) $
(H $ \NoIntersection -> Nothing)
:& (H $ \(a :: Point 2 r) -> Just $ Left a)
:& (H $ \(e :: LineSegment 2 p r) -> Just $ Right e)
:& RNil
type instance IntersectionOf (Line 3 r) (Triangle 3 p r) =
[ NoIntersection, Point 3 r, LineSegment 3 () r ]
instance (Fractional r, Ord r) => (Line 3 r) `IsIntersectableWith` (Triangle 3 p r) where
nonEmptyIntersection = defaultNonEmptyIntersection
l@(Line a v) `intersect` t@(Triangle (p :+ _) (q :+ _) (r :+ _)) =
match (l `intersect` h) $
(H $ \NoIntersection -> coRec NoIntersection)
:& (H $ \i@(Point3 _ _ _) -> if onTriangle' i then coRec i else coRec NoIntersection)
:& (H $ \_ -> intersect2d)
:& RNil
where
h@(Plane _ n) = supportingPlane t
-- 2d triangle and the line in terms of 2d-coordinates wr.t. of a
-- coordinate system in the supporting plane of t. The origin of this
-- coordinate system corresponds to the second vertex of t (q)
t' = Triangle (ext $ project p) (ext origin) (ext $ project r)
l' = Line (project a) (project' v)
-- test if the point in terms of its 2d coords lies in side the projected triangle
onTriangle' :: Point 3 r -> Bool
onTriangle' i = (project i) `onTriangle` t'
-- FIXME! these vectors may not be unit vectors. How do we deal with
-- that? (and does that really matter here?)
transf :: Transformation 3 r
transf = let u = p .-. q
in rotateTo (Vector3 u (n `cross` u) n) |.| translation ((-1) *^ (toVec q))
-- inverse of the transformation above.
invTrans :: Transformation 3 r
invTrans = inverseOf transf
project :: Point 3 r -> Point 2 r
project = projectPoint . transformBy transf
project' :: Vector 3 r -> Vector 2 r
project' = toVec . project . Point
lift :: Point 2 r -> Point 3 r
lift = Point . transformBy invTrans . flip V.snoc 0 . toVec
-- lift a 2d point back into plane coordinates
intersect2d :: Intersection (Line 3 r) (Triangle 3 p r)
intersect2d = match (l' `intersect` t') $
(H $ \NoIntersection -> coRec NoIntersection)
:& (H $ \i@(Point2 _ _) -> coRec $ lift i)
:& (H $ \(LineSegment s e) -> coRec $ LineSegment (s&unEndPoint.core %~ lift)
(e&unEndPoint.core %~ lift))
:& RNil