Naive.hs

--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.ConvexHull.Naive
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--------------------------------------------------------------------------------
module Algorithms.Geometry.ConvexHull.Naive( ConvexHull
                                           , lowerHull', lowerHullAll

                                           , isValidTriangle, upperHalfSpaceOf
                                           ) where

import           Control.Lens
import           Data.Ext
import           Data.Foldable (toList)
import           Data.Geometry.HalfSpace
import           Data.Geometry.HyperPlane
import           Data.Geometry.Line
import           Data.Geometry.Point
import           Data.Geometry.Triangle
import           Data.Geometry.Vector
import           Data.Intersection(intersects)
import           Data.List.NonEmpty (NonEmpty(..))
import           Data.List (find)
import           Data.Maybe (isNothing)
import           Data.Util
--------------------------------------------------------------------------------

type ConvexHull d p r = [Triangle 3 p r]

-- | Computes the lower hull without its vertical triangles.
--
-- pre: The points are in general position. In particular, no four
-- points should be coplanar.
--
-- running time: \(O(n^4)\)
lowerHull' :: forall r p. (Ord r, Fractional r, Show r)
           => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r
lowerHull' = filter (not . isVertical) . lowerHullAll
  where
    isVertical (Triangle p q r) =
      ccw' (p&core %~ projectPoint) (q&core %~ projectPoint) (r&core %~ projectPoint) == CoLinear

-- | Generates a set of triangles to be used to construct a complete
-- convex hull. In particular, it may contain vertical triangles.
--
-- pre: The points are in general position. In particular, no four
-- points should be coplanar.
--
-- running time: \(O(n^4)\)
lowerHullAll                 :: forall r p. (Ord r, Fractional r, Show r)
                             => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r
lowerHullAll (toList -> pts) = let mkT (Three p q r) = Triangle p q r in
    [ t | t <- mkT <$> uniqueTriplets pts, isNothing (isValidTriangle t pts) ]



_killOverlapping :: (Ord r, Fractional r) => [Triangle 3 p r] -> [Triangle 3 p r]
_killOverlapping = foldr keepIfNotOverlaps []
  where
    keepIfNotOverlaps t ts | any (t `overlaps`) ts = ts
                           | otherwise             = t:ts

overlaps :: (Fractional r, Ord r) => Triangle 3 p1 r -> Triangle 3 p2 r -> Bool
t1 `overlaps` t2 = upperHalfSpaceOf t1 == upperHalfSpaceOf t2 && False



-- | Tests if this is a valid triangle for the lower envelope. That
-- is, if all point lie above the plane through these points. Returns
-- a Maybe; if the result is a Nothing the triangle is valid, if not
-- it returns a counter example.
--
-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0))
-- >>> isValidTriangle t [ext $ Point3 5 5 0]
-- Nothing
-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0))
-- >>> isValidTriangle t [ext $ Point3 5 5 (-10)]
-- Just (Point3 5 5 (-10) :+ ())
isValidTriangle   :: (Num r, Ord r)
                  => Triangle 3 p r -> [Point 3 r :+ q] -> Maybe (Point 3 r :+ q)
isValidTriangle t = find (\a -> not $ (a^.core) `intersects` h)
  where
    h = upperHalfSpaceOf t


-- | Computes the halfspace above the triangle.
--
-- >>> upperHalfSpaceOf (Triangle (ext $ origin) (ext $ Point3 10 0 0) (ext $ Point3 0 10 0))
-- HalfSpace {_boundingPlane = HyperPlane {_inPlane = Point3 0 0 0, _normalVec = Vector3 0 0 100}}
upperHalfSpaceOf                  :: (Ord r, Num r) => Triangle 3 p r -> HalfSpace 3 r
upperHalfSpaceOf (Triangle p q r) = HalfSpace h
  where
    h' = from3Points (p^.core) (q^.core) (r^.core)
    c  = p&core.zCoord -~ 1
    h  = if (c^.core) `liesBelow` h' then h' else h'&normalVec %~ ((-1) *^)
    a `liesBelow` plane = (a `onSideUpDown` plane) == Below