{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UndecidableInstances #-}
-- | Triangles in \(d\)-dimensional space.
module Data.Geometry.Triangle where

import           Control.DeepSeq              (NFData)
import           Control.Lens
import           Data.Bifoldable              (Bifoldable (bifoldMap))
import           Data.Bifunctor               (Bifunctor (first))
import           Data.Bitraversable
import           Data.Either                  (partitionEithers)
import           Data.Ext
import           Data.Geometry.Ball           (Disk, disk)
import           Data.Geometry.Boundary       (PointLocationResult (..))
import           Data.Geometry.Box            (IsBoxable (..))
import           Data.Geometry.HyperPlane
import           Data.Geometry.Line           (Line (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.Util                    (Three, pattern Three)
import           Data.Vinyl                   (Rec (RNil, (:&)))
import           Data.Vinyl.CoRec             (Handler (H), match)
import           GHC.Generics                 (Generic)
import           GHC.TypeLits                 (type (+))

--------------------------------------------------------------------------------

-- | A triangle in \(d\)-dimensional space.
data Triangle d p r = Triangle !(Point d r :+ p)
                               !(Point d r :+ p)
                               !(Point d r :+ p)
                      deriving (Generic)

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, NFData r, NFData p) => NFData (Triangle d p r)

instance Arity d => Bifunctor  (Triangle d) where bimap = bimapDefault
instance Arity d => Bifoldable (Triangle d) where bifoldMap = bifoldMapDefault

instance Arity d => Bitraversable (Triangle d) where
  bitraverse f g (Triangle p q r) = let tr = bitraverse (traverse g) f in
    Triangle <$> tr p <*> tr q <*> tr 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)

instance Field1 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where
  _1 = lens (\(Triangle p _ _) -> p) (\(Triangle _ q r) p -> Triangle p q r)
instance Field2 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where
  _2 = lens (\(Triangle _ q _) -> q) (\(Triangle p _ r) q -> Triangle p q r)
instance Field3 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where
  _3 = lens (\(Triangle _ _ r) -> r) (\(Triangle p q _) r -> Triangle p q r)

type instance NumType   (Triangle d p r) = r
type instance Dimension (Triangle d p r) = d

-- | A \(d\)-dimensional triangle is isomorphic to a triple of \(d\)-dimensional points.
_TriangleThreePoints :: Iso' (Triangle d p r) (Three (Point d r :+ p))
_TriangleThreePoints = iso (\(Triangle p q r) -> Three p q r) (\(Three p q r) -> Triangle p q r)

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.
pattern Triangle' :: Point d r -> Point d r -> Point d r -> Triangle d () r
pattern Triangle' p q r <- Triangle (p :+ ()) (q :+ ()) (r :+ ())
  where
    Triangle' p q r = Triangle (ext p) (ext q) (ext r)

-- | Get the three line-segments that make up the sides of a triangle.
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

inTriangleRelaxed     :: (Ord r, Num r)
                 => Point 2 r -> Triangle 2 p r -> PointLocationResult
inTriangleRelaxed q (Triangle a b c)
    | ab == CoLinear && bc == ca = OnBoundary
    | bc == CoLinear && ca == ab = OnBoundary
    | ca == CoLinear && bc == ab = OnBoundary
    | ab == bc && bc == ca       = Inside
    | otherwise                  = Outside
  where
    ab = ccw (a^.core) (b^.core) q
    bc = ccw (b^.core) (c^.core) q
    ca = ccw (c^.core) (a^.core) q

-- | 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]

onTriangleRelaxed :: (Ord r, Num r) => Point 2 r -> Triangle 2 p r -> Bool
q `onTriangleRelaxed` t = inTriangleRelaxed q t /= Outside

-- 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 ]

{- HLINT ignore "Use const" -}
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

instance (Arity d, Ord r) => IsBoxable (Triangle d p r) where
  boundingBox (Triangle a b c) = boundingBox a <> boundingBox b <> boundingBox c