[Haskell] 3d or Nd geometry library?

[Axel Simon]

On Fri, Feb 13, 2004 at 11:31:39AM -0800, Abraham Egnor wrote:
> I was somewhat surprised to see that there's only one geometry library on
> the haskell libraries page, and further dismayed to find that it for the
> most part only does 2d.  It seems like haskell should be a natural fit for
> higher-order geometric libraries - has anyone heard of such?

What kind of geometry algorithms are you searching for? Perhaps the answer 
to you question is that computational geometry is hard to get right and 
thus deserves to be written in C or C++ for better accessiblity. In case 
you are interested in convex polyhedra then I can provide you with a 
Haskell binding to the Parma Polyhedra Library (search for Parma PPL). The 
binding is a bit still buggy, probably because I didn't manage to marshal 
Integers propperly. When I find the flaw, the binding will come with the  
library (like the O'Caml interface).

But you're right, especially due to the built-in arbitrary precision 
integers and rationals, Haskell is ideal for computational geometry. I 
attach the Graham convex hull algorithm, just for the sake of it.

Axel.

-------------- next part --------------
module Graham where

import List
import Ratio
--import Debug.Trace	(trace)

type Vertex = (Rational, Rational)

-- Removes double elements in a sorted list.
remDoubles :: Eq a => [a] -> [a]
remDoubles [] = []
remDoubles (x:xs) = x:rD x xs
  where
    rD _ [] = []
    rD last (x:xs) | x==last = rD last xs
                   | otherwise = x:rD x xs
 
-- Extract the point with the smallest y and largest x coordinate.
getHullVertex :: [Vertex] -> Vertex
getHullVertex [] = error "getHullVertex: empty set of points"
getHullVertex ps = foldl1 maxVertex ps
  where
    maxVertex :: Vertex -> Vertex -> Vertex
    maxVertex p1@(x1, y1) p2@(x2, y2)
      | y1<y2 = p1
      | y2<y1 = p2
      | x1>x2 = p1
      | otherwise = p2


-- Test three points for being ordered clockwise. If the argument simple is
-- false pair this test lexicographically with the distance of the second
-- and the thirds point.
isCW :: Vertex -> Vertex -> Vertex -> Ordering
isCW (a, b) (c, d) (e, f) = let
  dx21 = c-a
  dy21 = d-b
  dx31 = e-a
  dy31 = f-b
  in case compare (dy31*dx21) (dy21*dx31) of
    LT -> LT
    GT -> GT
    EQ -> compare (dx21*dx21+dy21*dy21) (dx31*dx31+dy31*dy31)


  
-- Sort the list of points according to one point on the hull.
preprocess :: [Vertex] -> [Vertex]
preprocess [] = []
preprocess pts = mp:sortBy (isCW mp) (filter ((/=) mp) pts)
 where
   mp = getHullVertex pts

-- Remove points in the interior.
graham :: [Vertex] -> [Vertex]
graham = gh . remDoubles . preprocess
  where
    gh points@[] = points
    gh points@[_] = points
    gh points@[_,_] = points
    gh (p:o:ints) = init $ walk [o,p] (ints++[p])

    -- Unfortunately lists all grow to the left in Haskell. Imagine the
    -- first argument to be written in reverse order. We are traversing
    -- the point set clockwise, so that the result does not have to be
    -- reversed.
    walk :: [Vertex] -> [Vertex] -> [Vertex]
    walk ps [] = ps
    walk (pM:pL:out) (pF:inp) | isCW pL pM pF/=GT 
				          = walk (pF:pM:pL:out) inp
			      | otherwise = walk (pL:out) (pF:inp)
    walk point@[_] (pF:inp) = walk (pF:point) inp
    walk ps ps' = error $ "weird point sets: "++show (length ps)++" and "++
			  show (length ps')++" long."