[Haskell] 3d or Nd geometry library?
A.Simon at kent.ac.uk
Mon Feb 16 09:18:58 EST 2004
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.
-------------- next part --------------
module Graham where
--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
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
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)
mp = getHullVertex pts
-- Remove points in the interior.
graham :: [Vertex] -> [Vertex]
graham = gh . remDoubles . preprocess
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
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."
More information about the Haskell