{-# 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