polygon triangulation (triangle test TODO), fixes #21

ofmooseandmen · ofmooseandmen · Sep 14, 2020 · Sep 8, 2020 · Sep 9, 2020 · Sep 11, 2020
commit c03fa260143342f6df665cf0715e90954f696a12
diff --git a/src/Data/Geo/Jord/Polygon.hs b/src/Data/Geo/Jord/Polygon.hs
@@ -22,6 +22,7 @@ module Data.Geo.Jord.Polygon
     , edges
     , concave
     -- * smart constructors
+    , Error(..)
     , simple
     , circle
     , arc
@@ -30,6 +31,7 @@ module Data.Geo.Jord.Polygon
     , triangulate
     ) where
 
+import Data.List (find)
 import Data.Maybe (isJust, mapMaybe)
 
 import Data.Geo.Jord.Angle (Angle)
@@ -44,6 +46,7 @@ import qualified Data.Geo.Jord.Length as Length
 import qualified Data.Geo.Jord.Math3d as Math3d
 import Data.Geo.Jord.Model (Spherical, surface)
 import Data.Geo.Jord.Triangle (Triangle)
+import qualified Data.Geo.Jord.Triangle as Triangle
 
 -- | A polygon whose vertices are horizontal geodetic positions.
 data Polygon a =
@@ -54,28 +57,35 @@ data Polygon a =
         }
     deriving (Eq, Show)
 
+data Error
+    = NotEnoughVertices -- ^ less than 3 vertices were supplied
+    | InvalidEdge -- ^ 2 consecutives vertices are antipodal or equal
+    | InvalidRadius -- ^ radius of circle or arc is <= 0
+    | EmptyArcRange -- ^ arc start angle == end angle
+    | SeflIntersectingEdge -- ^ 2 edges of the polygon intersect
+    deriving (Eq, Show)
+
 -- | Simple polygon (outer ring only and not self-intersecting) from given vertices. Returns an error ('Left') if:
 --
 --     * less than 3 vertices are given.
 --     * the given vertices defines self-intersecting edges.
 --     * the given vertices contains duplicated positions or antipodal positions.
-simple :: (Spherical a) => [HorizontalPosition a] -> Either String (Polygon a)
+simple :: (Spherical a) => [HorizontalPosition a] -> Either Error (Polygon a)
 simple vs
-    | null vs = Left "no vertex"
+    | null vs = Left NotEnoughVertices
     | head vs == last vs = simple (init vs)
-    | length vs < 3 = Left "not enough vertices"
-    -- FIXME: no equal/antipodal positions
+    | length vs < 3 = Left NotEnoughVertices
     | otherwise = simple' vs
 
 -- | Circle from given centre and radius. The resulting polygon contains @nb@ vertices equally distant from one
 -- another. Returns an error ('Left') if:
 --
 --     * given radius is 0
 --     * given number of positions is less than 3
-circle :: (Spherical a) => HorizontalPosition a -> Length -> Int -> Either String (Polygon a)
+circle :: (Spherical a) => HorizontalPosition a -> Length -> Int -> Either Error (Polygon a)
 circle c r nb
-    | r <= Length.zero = Left "invalid radius"
-    | nb < 3 = Left "invalid number of positions"
+    | r <= Length.zero = Left InvalidRadius
+    | nb < 3 = Left NotEnoughVertices
     | otherwise = Right (discretiseArc c r as)
   where
     n = fromIntegral nb :: Double
@@ -93,11 +103,11 @@ arc :: (Spherical a)
     -> Angle
     -> Angle
     -> Int
-    -> Either String (Polygon a)
+    -> Either Error (Polygon a)
 arc c r sa ea nb
-    | r <= Length.zero = Left "invalid radius"
-    | nb < 3 = Left "invalid number of positions"
-    | range == Angle.zero = Left "invalid range"
+    | r <= Length.zero = Left InvalidRadius
+    | nb < 3 = Left NotEnoughVertices
+    | range == Angle.zero = Left EmptyArcRange
     | otherwise = Right (discretiseArc c r as)
   where
     range = Angle.clockwiseDifference sa ea
@@ -112,9 +122,54 @@ contains :: (Spherical a) => Polygon a -> HorizontalPosition a -> Bool
 contains poly p = GreatCircle.enclosedBy p (vertices poly)
 
 triangulate :: (Spherical a) => Polygon a -> [Triangle a]
-triangulate _ = []
+triangulate p
+    | length vs == 3 = [triangle vs]
+    | otherwise = earClipping vs []
+  where
+    vs = vertices p
 
 -- private
+triangle :: (Spherical a) => [HorizontalPosition a] -> Triangle a
+triangle vs = Triangle.unsafeMake (head vs) (vs !! 1) (vs !! 2)
+
+earClipping :: (Spherical a) => [HorizontalPosition a] -> [Triangle a] -> [Triangle a]
+earClipping vs ts
+    | length vs == 3 = reverse (triangle vs : ts)
+    | otherwise =
+        case (findEar vs) of
+            Nothing -> []
+            (Just (p, e, n)) -> earClipping vs' ts'
+                where ts' = Triangle.unsafeMake p e n : ts
+                      vs' = filter (\v -> v /= e) vs
+
+findEar ::
+       (Spherical a)
+    => [HorizontalPosition a]
+    -> Maybe (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
+findEar ps = find (\c -> isEar c rs) convex
+  where
+    rs = reflices ps
+    t3 = tuple3 ps
+    convex = filter (\(_, v, _) -> (v `notElem` rs)) t3
+
+-- | a convex vertex @c@ is an ear if triangle (prev, c, next) contains no reflex.
+isEar ::
+       (Spherical a)
+    => (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
+    -> [HorizontalPosition a]
+    -> Bool
+isEar (p, c, n) rs = all (\r -> not (GreatCircle.enclosedBy r vs)) rs
+  where
+    vs = [p, c, n]
+
+-- | A reflex is a vertex where the polygon is concave.
+-- a vertex is a reflex if previous vertex is left (assuming a clockwise polygon), otherwise it is a convex vertex
+reflices :: (Spherical a) => [HorizontalPosition a] -> [HorizontalPosition a]
+reflices ps = fmap (\(_, c, _) -> c) rs
+  where
+    t3 = tuple3 ps
+    rs = filter (\(p, c, n) -> GreatCircle.side p c n == GreatCircle.LeftOf) t3
+
 -- | [mAB, mBC, mCD, mDE, mEA]
 -- no intersections:
 -- mAB vs [mCD, mDE]
@@ -146,11 +201,13 @@ makePairs' (xs, ps)
             else tail . tail $ xs
     np = (head xs, versus)
 
-simple' :: (Spherical a) => [HorizontalPosition a] -> Either String (Polygon a)
+simple' :: (Spherical a) => [HorizontalPosition a] -> Either Error (Polygon a)
 simple' vs =
-    if si
-        then Left "self-intersecting edges"
-        else Right (Polygon os es isConcave)
+    if length es /= length vs
+        then Left InvalidEdge
+        else if si
+                 then Left SeflIntersectingEdge
+                 else Right (Polygon os es isConcave)
   where
     zs = tuple3 vs
     clockwise = sum (fmap (\(a, b, c) -> Angle.toRadians (GreatCircle.turn a b c)) zs) < 0.0
@@ -172,7 +229,7 @@ tuple3 ps = zip3 l1 l2 l3
   where
     l1 = last ps : init ps
     l2 = ps
-    l3 = tail ps ++ [head ps] ++ [head ps]
+    l3 = tail ps ++ [head ps]
 
 mkEdges :: (Spherical a) => [HorizontalPosition a] -> [MinorArc a]
 mkEdges ps = mapMaybe (uncurry GreatCircle.minorArc) (zip ps (tail ps ++ [head ps]))

diff --git a/src/Data/Geo/Jord/Triangle.hs b/src/Data/Geo/Jord/Triangle.hs
@@ -6,49 +6,102 @@
 -- Stability:   experimental
 -- Portability: portable
 --
--- TODO
+-- Types and functions for working with triangles at the surface of a __spherical__ celestial body.
+--
+-- In order to use this module you should start with the following imports:
+--
+-- @
+-- import qualified Data.Geo.Jord.Geodetic as Geodetic
+-- import qualified Data.Geo.Jord.Triangle as Triangle
+-- @
 module Data.Geo.Jord.Triangle
     ( Triangle
+    , vertex0
     , vertex1
     , vertex2
-    , vertex3
-    , edge1
-    , edge2
-    , edge3
+    , make
+    , unsafeMake
     , centroid
     , circumcentre
-    , make
     , contains
-    , circumcircleContains
     ) where
 
+import Control.Monad (join)
+import Data.Maybe (fromJust)
+
 import Data.Geo.Jord.Geodetic (HorizontalPosition)
-import Data.Geo.Jord.GreatCircle (MinorArc)
+import qualified Data.Geo.Jord.Geodetic as Geodetic
+import qualified Data.Geo.Jord.GreatCircle as GreatCircle
+import qualified Data.Geo.Jord.Math3d as Math3d
 import Data.Geo.Jord.Model (Spherical)
 
+-- | A triangle whose vertices are horizontal geodetic positions.
 data Triangle a =
-    Triangle
-        { vertex1 :: HorizontalPosition a
-        , vertex2 :: HorizontalPosition a
-        , vertex3 :: HorizontalPosition a
-        , edge1 :: MinorArc a
-        , edge2 :: MinorArc a
-        , edge3 :: MinorArc a
-        , centroid :: HorizontalPosition a
-        , circumcentre :: HorizontalPosition a
-        }
+    Triangle (HorizontalPosition a) (HorizontalPosition a) (HorizontalPosition a)
     deriving (Eq, Show)
 
+-- | First vertex of given triangle.
+vertex0 :: (Spherical a) => Triangle a -> HorizontalPosition a
+vertex0 (Triangle v _ _) = v
+
+-- | Second vertex of given triangle.
+vertex1 :: (Spherical a) => Triangle a -> HorizontalPosition a
+vertex1 (Triangle _ v _) = v
+
+-- | Third vertex of given triangle.
+vertex2 :: (Spherical a) => Triangle a -> HorizontalPosition a
+vertex2 (Triangle _ _ v) = v
+
+-- | Triangle from given vertices. Returns 'Nothing' if some vertices are equal or some are antipodes of others.
 make ::
        (Spherical a)
     => HorizontalPosition a
     -> HorizontalPosition a
     -> HorizontalPosition a
-    -> Either String (Triangle a)
-make _ _ _ = Left "TODO"
+    -> Maybe (Triangle a)
+make v0 v1 v2
+    | v0 == v1 || v1 == v2 || v2 == v0 = Nothing
+    | a0 == v1 || a0 == v2 || a1 == v2 = Nothing
+    | otherwise = Just (Triangle v0 v1 v2)
+  where
+    a0 = Geodetic.antipode v0
+    a1 = Geodetic.antipode v1
 
+-- | Triangle from given vertices. This is unsafe, if any vertices are equal or some are antipodes of others,
+-- the resulting triangle is actually undefined.
+unsafeMake ::
+       (Spherical a)
+    => HorizontalPosition a
+    -> HorizontalPosition a
+    -> HorizontalPosition a
+    -> Triangle a
+unsafeMake = Triangle
+
+-- | @contains t p@ returns 'True' if position @p@ is enclosed by the vertices of triangle
+-- @t@ - see 'GreatCircle.enclosedBy'.
 contains :: (Spherical a) => Triangle a -> HorizontalPosition a -> Bool
-contains _ _ = False
+contains (Triangle v0 v1 v2) p = GreatCircle.enclosedBy p [v0, v1, v2]
 
-circumcircleContains :: (Spherical a) => Triangle a -> HorizontalPosition a -> Bool
-circumcircleContains _ _ = False
+-- | Computes the centroid of the given triangle: the position which is the intersection of the three medians of
+-- the triangle (each median connecting a vertex with the midpoint of the opposite side).
+--
+-- The centroid is always within the triangle.
+centroid :: (Spherical a) => Triangle a -> HorizontalPosition a
+centroid (Triangle v0 v1 v2) = fromJust c -- this is safe unless triangle was created using unsafeMake.
+  where
+    m1 = GreatCircle.mean [v1, v2] >>= GreatCircle.minorArc v0
+    m2 = GreatCircle.mean [v2, v0] >>= GreatCircle.minorArc v1
+    c = join (GreatCircle.intersection <$> m1 <*> m2)
+
+-- | The circumcentre of the triangle: the position which is equidistant from all three vertices.
+--
+-- The circumscribed circle or circumcircle of a triangle is a circle which passes through all
+-- the vertices of the triangle; The circumcentre is not necessarily inside the triangle.
+--
+-- Thanks to STRIPACK: http://orion.math.iastate.edu/burkardt/f_src/stripack/stripack.f90
+circumcentre :: (Spherical a) => Triangle a -> HorizontalPosition a
+circumcentre (Triangle v0 v1 v2) = Geodetic.nvectorPos' cu (Geodetic.model v0)
+  where
+    e0 = Math3d.subtract (Geodetic.nvector v1) (Geodetic.nvector v0)
+    e1 = Math3d.subtract (Geodetic.nvector v2) (Geodetic.nvector v0)
+    cu = Math3d.cross e0 e1