My old 2D Delaunay/Voronoi tessellation tools

Delaunay.hs

{-# LANGUAGE TupleSections #-}

module Delaunay where

import qualified Data.Function         as F
import qualified Data.List             as L
import           Data.Map.Strict       (Map)
import qualified Data.Map.Strict       as Map
import           Data.Set              (Set)
import qualified Data.Set              as Set
import qualified Numeric.LinearAlgebra as LA

type Coord = Double

-- | An angle, in radians
newtype Angle = Angle { unAngle :: Double }
              deriving (Eq, Ord, Show)

-- | A single point in 2D Euclidean space
data Point = Point { pointX :: Coord
                   , pointY :: Coord
                   } deriving (Eq, Ord)

-- | A line-segment consisting of two points
data Edge = Edge Point Point
          deriving (Eq, Ord, Show)

-- | A triangle, defined b three points
data Tri = Tri Point Point Point
         deriving (Eq, Ord, Show)

-- | A polygon, consisting of an arbitrary number of points.  This
-- polygon should form a convex hull or a non-intersecting star.
data Poly = Poly (Set Point)
          deriving (Eq, Ord, Show)

-- | A non-tesselated cloud of points, which may have some properties
-- such as the centroid.
data PointCloud = PointCloud (Set Point)
                deriving (Eq, Ord, Show)

-- | Coerce a Poly to Tri, presuming it has exactly 3 points.
triPoly :: Poly -> Tri
triPoly p =
  case ccwVertexList p of
   [a, b, c] -> Tri a b c
   _ -> error "triangles must be exactly 3 points"

-- | Coerce a Tri to Poly
polyTri :: Tri -> Poly
polyTri (Tri a b c) = Poly $ Set.fromList [a, b, c]

-- | A class describing similarities in computable geometric
-- structures
class Shape s where
  centroid :: s -> Point -- ^ The barycenter of the shape
  vertices :: s -> Set Point -- ^ Set of the vertices
  vertexCount :: s -> Int -- ^ The count of vertices
  edges :: s -> Set Edge -- ^ Set of the edges
  edgeCount :: s -> Int -- ^ The count of edges
  edgesAbout :: s -> Point -> [Edge] -- ^ CCW ordered set of edges seen from point
  cenX :: s -> Coord -- ^ The x coordinate of the centroid
  cenY :: s -> Coord -- ^ The y coordinate of the centroid
  angleTo :: (Shape a) => s -> a -> Angle -- ^ Angle from the centroid of a to b
  ccwVertexList :: s -> [Point] -- ^ List of vertices, sorted counter-clockwise
  ccwVertexListAbout :: s -> Point -> [Point]
                       -- ^ List of vertices, sorted counter-clockwise
  edgeAdjacent :: (Shape a) => s -> a -> Bool
                 -- ^ Determine if shape a shares at least two vertices
  area :: s -> Coord
         -- ^ The area covered by this shape

  -- minimal complete is just vertices
  vertexCount s = Set.size $ vertices s
  edgeCount s = Set.size $ edges s
  centroid s | (length vl) > 0 = Point (meanCoord pointX) (meanCoord pointY)
    where
      vl = Set.toList $ vertices s
      meanCoord accessor = (sum $ map accessor vl) / (fromIntegral $ length vl)
  centroid _ = error "shapes with no vertices have no computable centroid"
  cenX s = pointX $ centroid s
  cenY s = pointY $ centroid s
  angleTo a b = Angle $ atan2 ((cenY b)-(cenY a)) ((cenX b)-(cenX a))
  edges s = Set.fromList $ s `edgesAbout` (centroid s)
  edgesAbout s _ | (Set.size $ vertices s) < 2 = []
  edgesAbout s c = take (vertexCount s) $ zipWith Edge vertexRing $ tail vertexRing
    where
      vertexRing = cycle $ s `ccwVertexListAbout` c
  ccwVertexList s = ccwVertexListAbout s (centroid s)
  ccwVertexListAbout s c =
    map snd
    $ L.sortBy (compare `F.on` fst)
    $ map (\p-> (c `angleTo` p, p))
    $ Set.toList
    $ vertices s
  edgeAdjacent a b =
    (Set.size $ Set.intersection (vertices a) (vertices b)) >= 2
  area _ = 0

instance Show Point where
  showsPrec _ (Point px py) suf = '<':(shows px $ ',': (shows py $ '>':suf))

instance Shape Point where
  centroid = id
  vertices = Set.singleton

instance Num Point where
  a + b = Point ((pointX a) + (pointX b)) ((pointY a) + (pointY b))
  a * b = Point ((pointX a) * (pointX b) + (pointY a) * (pointY b)) 0 -- dot product in x
  negate a = Point (negate $ pointX a) (negate $ pointY a)
  abs a = Point (sqrt((pointX a)^(2::Int) + (pointY a)^(2::Int))) 0
  signum _ = Point 0 0
  fromInteger a = Point (fromInteger a) 0

instance Shape Edge where
  vertices (Edge a b) = Set.fromList [a, b]

edgeLen :: Edge -> Coord
edgeLen (Edge a b) = pointX $ abs $ a-b

instance Shape Tri where
  vertices (Tri a b c) = Set.fromList [a, b, c]
  -- https://en.wikipedia.org/wiki/Heron%27s_formula
  area (Tri a b c) = sqrt $ s * (s-ab) * (s-bc) * (s-ca)
    where
      s = (ab + bc + ca) / 2
      ab = edgeLen $ Edge a b
      bc = edgeLen $ Edge b c
      ca = edgeLen $ Edge c a
  centroid (Tri a b c) =
    let midBC = centroid $ PointCloud $ Set.fromList [b, c]
    in centroid $ PointCloud $ Set.fromList [a, midBC, midBC]

instance Shape Poly where
  vertices (Poly ps) = ps
  area (Poly ps) = sum $ map area $ Set.toList $ starTriangulate (centroid $ PointCloud ps) (Poly ps)
  centroid p@(Poly ps) =
    case Set.toList $ vertices p of
      [] -> Point 0 0
      [a] -> a
      [a, b] -> centroid $ Edge a b
      [a, b, c] -> centroid $ Tri a b c
      _ -> case (sum $ map (\t->case (centroid t, area t) of
                                (Point cx cy, at) -> Point (cx*at) (cy*at)) starTris) of
            Point sumx sumy -> Point (sumx/polyArea) (sumy/polyArea)
    where
      polyArea = area $ Poly ps
      starTris = Set.toList $ starTriangulate (centroid $ PointCloud ps) (Poly ps)

-- https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
-- a http://paulbourke.net/geometry/polygonmesh/ in particular

instance Shape PointCloud where
  vertices (PointCloud ps) = ps
  edgesAbout _ _ = []

-- | See if Point is contained in the circumcircle of Tri. Colinearly
-- degenerate triangles are defined to include all points.
circContains :: Tri -> Point -> Bool
circContains t _ | triColinear t = True
circContains t (Point dx dy) =
  case ccwVertexList t of
   [Point ax ay, Point bx by, Point cx cy] ->
     (LA.det $ LA.fromLists [[ax, ay, ax*ax+ay*ay, 1],
                             [bx, by, bx*bx+by*by, 1],
                             [cx, cy, cx*cx+cy*cy, 1],
                             [dx, dy, dx*dx+dy*dy, 1]]) > 0
   _ -> True --  error "circContains received a degenerate triangle"

-- | See if the points of the triangle are colinear.
triColinear :: Tri -> Bool
triColinear tri =
  (/=3) $ Set.size $ Set.map edgeAngle $ edges tri
  where
    edgeAngle (Edge ea eb) = ea `angleTo` eb

-- | Identify triangles that no longer satisfy the Delaunay constraint
-- because their circumcircles contain the point.
badTriangles :: Point -> (Set Tri) -> (Set Tri)
badTriangles p tris = Set.filter (`circContains` p) tris

-- | Trace the border of a set of connected polygons
traceBorder :: (Shape s) => (Set s) -> Poly
traceBorder tris =
  Poly $ Set.fromList $ concatMap (Set.toList . vertices)
  $ Map.keys
  $ Map.filter (==(1::Int))
  $ Map.fromListWith (+)
  $ map (,1)
  $ concatMap (Set.toList . edges)
  $ Set.toList
  $ tris

-- | Radially triangulate a star-polygon with triangles sharing vertex
-- Point.
starTriangulate :: Point -> Poly -> Set Tri
starTriangulate p poly = Set.fromList $ map makeTri $ poly `edgesAbout` p
  where
    makeTri (Edge a b) = Tri a b p

-- | Calculate a triangle large enough to include all points in the
-- pointcloud.
constructSuperTri :: PointCloud -> Tri
constructSuperTri pc@(PointCloud ps) =
  Tri a b c
  where
    cen = centroid pc
    distCen p = pointX $ abs $ p - cen
    r = 8 + (Set.findMax $ Set.map distCen ps)
    a = cen + (Point (-20*r) (30*r))
    b = cen + (Point (-20*r) (-30*r))
    c = cen + (Point (30*r) 0)

-- | Calculate the center and radius of the circumcircle for a
-- triangle.
circumcircle :: Tri -> (Point, Coord)
circumcircle t | triColinear t =
                   error $ "colinear points from degenerate triangle "
                   ++ (show t) ++ " have no circumcircle"
circumcircle (Tri a b c) =
  case LA.toLists $ matAinv LA.<> matB of
   [[ccx],[ccy]] -> ((Point ccx ccy), pointX $ abs $ (Point ccx ccy)-a)
   x -> error $ "\n\ncircumcircle fault in " ++ (show $ Tri a b c)
       ++ "\n\nproduced " ++ (show x)
       ++ "\n\nfrom mat A " ++ (show matA)
       ++ "\n\ninv " ++ (show matAinv)
       ++ "\n\nmatB " ++ (show matB)
  where
    row :: Edge -> ([Double], [Double])
    row e@(Edge (Point ax ay) (Point bx by)) =
      if ay==by then ([0, 1], [ay])
      else if ax==bx then ([1, 0], [ax])
           else ([m, -1], [m*cx-cy])
      where
        (Point cx cy) = centroid e
        m = -(bx-ax)/(by-ay)
    rows = map row [Edge a b, Edge b c, Edge c a]
    pzrows = take 2 $ L.nub $ (filter (\([ra,rb],_)->ra==0||rb==0) rows)
             ++ (reverse rows)
    matA = LA.fromLists $ map fst pzrows
    matB = LA.fromLists $ map snd pzrows
    matAinv = LA.inv matA

-- | Add a point to an existing triangulation set, shattering and
-- reforming triangles as needed.
addPoint :: Point -> (Set Tri) -> (Set Tri)
addPoint p tris =
  (tris `Set.difference` badTris) `Set.union` newTris
  where
    badTris = badTriangles p tris
    newTris = starTriangulate p $ traceBorder badTris

-- | Calculate a <https://en.wikipedia.org/wiki/Delaunay_triangulation Delaunay triangulation>
-- of the pointcloud using the
-- <https://en.wikipedia.org/wiki/Bowyer%E2%80%93Watson_algorithm Bowyer-Watson algorithm>.
bowyerWatsonTriangulate :: PointCloud -> Set Tri
bowyerWatsonTriangulate pc@(PointCloud ps) =
  Set.filter notSuper $
  Set.foldr addPoint (Set.singleton superTri) ps
  where
    superTri = constructSuperTri pc
    notSuper :: Tri -> Bool
    notSuper t = Set.null $ (vertices t) `Set.intersection` (vertices superTri)

-- | Identify the adjacent triangles in a Delaunay triangulation given
-- a triangle.
delaunayNeighbors :: Set Tri -> Tri -> Set Tri
delaunayNeighbors tris tri = Set.filter (edgeAdjacent tri) tris

-- | Produce a Delaunay connectivity map.
delaunayConnectivity :: Set Tri -> Map Tri (Set Tri)
delaunayConnectivity tris = Map.fromSet (delaunayNeighbors tris) tris

-- | Produce a map of points identifying which triangles contain them.
vertexMembers :: Set Tri -> Map Point (Set Tri)
vertexMembers dTris =
  let vs = Set.unions $ map vertices $ Set.toList dTris
  in Map.fromSet (\v->Set.filter (\t->Set.member v $ vertices t) dTris) vs

-- | Produce the Voronoi diagram of a Delaunay triangulation.
voronoi :: Set Tri -> Set Poly
voronoi dTris =
  Set.fromList $ map (Poly . (Set.map (fst . circumcircle))) $ Map.elems $ vertexMembers dTris

-- | Perform Lloyd-relaxation on a set of polygons using their
-- centroid.
lloydRelax :: Set Poly -> PointCloud
lloydRelax ps = PointCloud $ Set.map centroid ps