diff --git a/DistanceTransform.cabal b/DistanceTransform.cabal
new file mode 100644
--- /dev/null
+++ b/DistanceTransform.cabal
@@ -0,0 +1,43 @@
+name:                DistanceTransform
+version:             0.1.2
+synopsis:            Distance transform function.
+description:         An n-D distance transform that computes the Euclidean
+                     distance between each element in a discrete field and the nearest cell
+                     containing a zero.
+                     .
+                     The algorithm implemented is based off of
+                     Meijster et al., /"A general algorithm for computing distance/
+                     /transforms in linear time."/ Parallel versions of both the Euclidean
+                     distance transform and squared Euclidean distance transform are also
+                     provided.
+license:             BSD3
+license-file:        LICENSE
+author:              Anthony Cowley
+maintainer:          acowley@gmail.com
+copyright:           (c) Anthony Cowley 2012,2013
+category:            Math
+build-type:          Simple
+cabal-version:       >=1.10
+extra-source-files:  src/tests/Main.hs src/tests/TestPar.hs
+
+source-repository head
+  type: git
+  location: git://github.com/acowley/DistanceTransform.git
+
+library
+  exposed-modules:     DistanceTransform.Euclidean,
+                       DistanceTransform.Internal.Indexer
+  build-depends:       base >=4.5 && < 5, vector >=0.9, primitive
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+  ghc-options:         -O2 -Wall
+
+test-suite tests
+  type: exitcode-stdio-1.0
+  hs-source-dirs: src/tests
+  main-is: Main.hs
+  ghc-options: -Wall -O2 -threaded -rtsopts
+  default-language: Haskell2010
+  build-depends: base >= 4.5 && < 5,
+                 test-framework, test-framework-hunit, HUnit,
+                 vector, DistanceTransform
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2012, Anthony Cowley
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Anthony Cowley nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/src/DistanceTransform/Euclidean.hs b/src/DistanceTransform/Euclidean.hs
new file mode 100644
--- /dev/null
+++ b/src/DistanceTransform/Euclidean.hs
@@ -0,0 +1,203 @@
+{-# LANGUAGE BangPatterns, FlexibleContexts, ScopedTypeVariables, 
+             RankNTypes #-}
+-- |N-dimensional parallel Euclidean distance transform using an
+-- approach derived from: Meijster et al., /"A general algorithm for/
+-- /computing distance transforms in linear time."/
+module DistanceTransform.Euclidean (edt, edtPar, sedt, sedtPar) where
+import Control.Monad (when)
+import Control.Monad.ST (ST)
+import Control.Monad.ST.Unsafe (unsafeIOToST, unsafeSTToIO)
+import Data.Vector.Generic (Vector, (!))
+import qualified Data.Vector.Generic as V
+import qualified Data.Vector.Generic.Mutable as VM
+import qualified Data.Vector.Storable as S
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as UM
+import Data.Word (Word8)
+import DistanceTransform.Internal.Indexer
+
+-- | Higher order function that runs an inner loop across the
+-- indicated dimension.
+type LoopRunner = forall s. Zipper Int -> (Int -> Int -> ST s ()) -> ST s ()
+
+-- This constructs Meijster's G function.
+phase1 :: (Integral a, Vector v a, Vector v Int)
+       => LoopRunner -> Zipper Int -> v a -> v Int
+phase1 runLoop dim p = 
+  V.map (\x -> x*x) $ V.create $
+   do v <- VM.new (product $ fromZipper dim)
+      let pullRight !i = if p ! i == 0
+                         then VM.unsafeWrite v i 0
+                         else VM.unsafeRead v (i-step) >>= 
+                                (VM.unsafeWrite v i $!) . (1+)
+          pushLeft !i = do !prev <- VM.unsafeRead v (i+step)
+                           !curr <- VM.unsafeRead v i
+                           when (prev < curr) 
+                                (VM.unsafeWrite v i $! prev+1)
+          innerLoop !offset _ = 
+            do VM.unsafeWrite v offset $! toInfty offset
+               mapM_ (pullRight . (offset+)) [step,2*step..n' - 1]
+               mapM_ (pushLeft . (offset+)) [n'-2*step,n'-3*step..0]
+      runLoop dim innerLoop
+      return v
+  where toInfty !i = let !dimsum = zipSum dim
+                     in if p ! i == 0 then 0 else dimsum
+        {-# INLINE toInfty #-}
+        step = zipStep dim
+        n = focus dim -- Get the actual dimension size
+        n' = n * step
+{-# SPECIALIZE phase1 :: 
+  LoopRunner -> Zipper Int -> U.Vector Int -> U.Vector Int #-}
+{-# SPECIALIZE phase1 :: 
+  LoopRunner -> Zipper Int -> U.Vector Word8 -> U.Vector Int #-}
+{-# SPECIALIZE phase1 :: 
+  LoopRunner -> Zipper Int -> S.Vector Int -> S.Vector Int #-}
+{-# SPECIALIZE phase1 :: 
+  LoopRunner -> Zipper Int -> S.Vector Word8 -> S.Vector Int #-}
+
+foldMfromStepTo :: (Eq b, Monad m) => 
+                   (a -> b -> m a) -> a -> b -> (b -> b) -> b -> m a
+foldMfromStepTo f z from step to = go from z
+  where to' = step to
+        go !x !acc = if x == to' then return acc else f acc x >>= go (step x)
+{-# INLINE foldMfromStepTo #-}
+
+-- Each phase needs the squared eucilidean distance from the previous
+-- phase.
+phaseN :: Vector v Int => Zipper Int -> v Int -> v Int
+phaseN dim sedt' = 
+  V.create $
+  do v <- VM.new $ V.length sedt'
+     zipFoldMAsYouDo dim (phaseNRow m sedt' v)
+     return v
+  where m = focus dim
+{-# SPECIALIZE phaseN :: Zipper Int -> U.Vector Int -> U.Vector Int #-}
+{-# SPECIALIZE phaseN :: Zipper Int -> S.Vector Int -> S.Vector Int #-}
+
+parPhaseN :: Vector v Int => Zipper Int -> v Int -> v Int
+parPhaseN dim sedt' = 
+  V.create $ 
+  do v <- VM.new $ V.length sedt'
+     unsafeIOToST $ 
+       parZipFoldMAsYouDo dim ((unsafeSTToIO .) . phaseNRow m sedt' v)
+     return v
+  where m = focus dim
+{-# SPECIALIZE parPhaseN :: Zipper Int -> U.Vector Int -> U.Vector Int #-}
+{-# SPECIALIZE parPhaseN :: Zipper Int -> S.Vector Int -> S.Vector Int #-}
+
+phaseNRow :: forall v mv s. (Vector v Int, VM.MVector mv Int)
+          => Int -> v Int -> mv s Int -> Int -> Int -> ST s ()
+phaseNRow m sedt' v offset step = 
+  do s <- UM.new m
+     t <- UM.new m
+     let {-# INLINE fMetric #-}
+         fMetric !x !i = let !d = x - i in d*d + gsq i
+         {-# INLINE sep #-}
+         -- I flipped the order of the arguments from Meijster's paper
+         -- for ease of use in scan3
+         sep !u !i = ((u*u-i*i+gsq u - gsq i) `quot` (2*(u-i))) + 1
+     VM.unsafeWrite s 0 0
+     VM.unsafeWrite t 0 0
+     let {-# INLINE qaux #-}
+         qaux :: Int -> (Int -> ST s Int) -> Int -> ST s Int
+         qaux !u k = goqaux
+           where goqaux !q | q < 0 = k q
+                           | otherwise = do !tq <- VM.unsafeRead t q
+                                            !sq <- VM.unsafeRead s q
+                                            if fMetric tq sq > fMetric tq u
+                                            then let !q' = q-1 in goqaux q'
+                                            else k q
+         scan3 !q0 !u = let {-# INLINE aux #-}
+                            aux !q = 
+                              if q < 0 
+                              then VM.unsafeWrite s 0 u >> return 0
+                              else do !w <- (sep u $!) `fmap` VM.unsafeRead s q
+                                      if w < m
+                                      then let !q' = q+1
+                                           in do VM.unsafeWrite s q' u
+                                                 VM.unsafeWrite t q' w
+                                                 return q'
+                                      else return q
+                        in qaux u aux q0
+         scan4 !q !u = do !sq <- VM.unsafeRead s q
+                          let !i = offset + u * step
+                          VM.unsafeWrite v i $! fMetric u sq
+                          !tq <- VM.unsafeRead t q
+                          if u == tq then let !q' = q-1 in return q' 
+                                     else return q
+     q <- foldMfromStepTo scan3 (0::Int) 1 (+1) (m-1)
+     _ <- foldMfromStepTo scan4 q (m-1) (subtract 1) (0::Int)
+     return ()
+  where gsq !i = sedt' ! (offset+step*i)
+        {-# INLINE gsq #-}
+{-# SPECIALIZE phaseNRow :: Int -> U.Vector Int -> U.MVector s Int
+                         -> Int -> Int -> ST s () #-}
+{-# SPECIALIZE phaseNRow :: Int -> S.Vector Int -> S.MVector s Int
+                         -> Int -> Int -> ST s () #-}
+
+-- |Compute the squared Euclidean distance transform of an
+-- N-dimensional array. Dimensions given as
+-- @[width,height,depth...]@. The left-most dimension is the
+-- inner-most.
+sedt :: (Vector v a, Vector v Int, Integral a) => [Int] -> v a -> v Int
+sedt dims p = go (left dim0) (phase1 zipFoldMAsYouDo dim0 p)
+  where dim0 = rightmost . unsafeToZipper $ reverse dims
+        go Nothing sedt' = sedt'
+        go (Just dim) sedt' = go (left dim) (phaseN dim sedt')
+{-# SPECIALIZE sedtPar :: [Int] -> U.Vector Int -> U.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> U.Vector Word8 -> U.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> S.Vector Int -> S.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> S.Vector Word8 -> S.Vector Int #-}
+
+-- |Compute the Euclidean distance transform of an N-dimensional
+-- array. Dimensions given as @[width,height,depth...]@. The left-most
+-- dimension is the inner-most. For an array representing a 2D
+-- collection in row-major format, we would give @[width,height]@ or
+-- @[columns,rows]@.
+edt :: (Integral a, Floating b, Vector v a, Vector v b, Vector v Int)
+    => [Int] -> v a -> v b
+edt dims v = V.map aux $ sedt dims v
+  where aux = sqrt . fromIntegral
+{-# SPECIALIZE edt :: [Int] -> U.Vector Int -> U.Vector Float #-}
+{-# SPECIALIZE edt :: [Int] -> U.Vector Int -> U.Vector Double #-}
+{-# SPECIALIZE edt :: [Int] -> U.Vector Word8 -> U.Vector Float #-}
+{-# SPECIALIZE edt :: [Int] -> U.Vector Word8 -> U.Vector Double #-}
+{-# SPECIALIZE edt :: [Int] -> S.Vector Int -> S.Vector Float #-}
+{-# SPECIALIZE edt :: [Int] -> S.Vector Int -> S.Vector Double #-}
+{-# SPECIALIZE edt :: [Int] -> S.Vector Word8 -> S.Vector Float #-}
+{-# SPECIALIZE edt :: [Int] -> S.Vector Word8 -> S.Vector Double #-}
+
+-- |Compute the Euclidean distance transform of an N-dimensional array
+-- using multiple processor cores. Dimensions given as
+-- @[width,height,depth...]@. The left-most dimension is the
+-- inner-most. For an array representing a 2D collection in row-major
+-- format, we would give @[width,height]@ or @[columns,rows]@.
+edtPar :: (Integral a, Floating b, Vector v a, Vector v b, Vector v Int)
+       => [Int] -> v a -> v b
+edtPar dims v = V.map aux $ sedtPar dims v
+  where aux = sqrt . fromIntegral
+{-# SPECIALIZE edtPar :: [Int] -> U.Vector Int -> U.Vector Float #-}
+{-# SPECIALIZE edtPar :: [Int] -> U.Vector Int -> U.Vector Double #-}
+{-# SPECIALIZE edtPar :: [Int] -> U.Vector Word8 -> U.Vector Float #-}
+{-# SPECIALIZE edtPar :: [Int] -> U.Vector Word8 -> U.Vector Double #-}
+{-# SPECIALIZE edtPar :: [Int] -> S.Vector Int -> S.Vector Float #-}
+{-# SPECIALIZE edtPar :: [Int] -> S.Vector Int -> S.Vector Double #-}
+{-# SPECIALIZE edtPar :: [Int] -> S.Vector Word8 -> S.Vector Float #-}
+{-# SPECIALIZE edtPar :: [Int] -> S.Vector Word8 -> S.Vector Double #-}
+
+-- |Compute the squared Euclidean distance transform of an
+-- N-dimensional array using multiple processor cores. Dimensions
+-- given as @[width,height,depth...]@. The left-most dimension is the
+-- inner-most.
+sedtPar :: (Vector v a, Vector v Int, Integral a) => [Int] -> v a -> v Int
+sedtPar dims p = go (left dim0) (phase1 parZipFoldMAsYouDo' dim0 p)
+  where dim0 = rightmost . unsafeToZipper $ reverse dims
+        go Nothing sedt' = sedt'
+        go (Just dim) sedt' = go (left dim) (parPhaseN dim sedt')
+        parZipFoldMAsYouDo' :: Zipper Int -> (Int -> Int -> ST s ()) -> ST s ()
+        parZipFoldMAsYouDo' z f = unsafeIOToST $ 
+                                  parZipFoldMAsYouDo z ((unsafeSTToIO .) .f)
+{-# SPECIALIZE sedtPar :: [Int] -> U.Vector Int -> U.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> U.Vector Word8 -> U.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> S.Vector Int -> S.Vector Int #-}
+{-# SPECIALIZE sedtPar :: [Int] -> S.Vector Word8 -> S.Vector Int #-}
diff --git a/src/DistanceTransform/Internal/Indexer.hs b/src/DistanceTransform/Internal/Indexer.hs
new file mode 100644
--- /dev/null
+++ b/src/DistanceTransform/Internal/Indexer.hs
@@ -0,0 +1,133 @@
+-- |Helpers for performing nested loop iteration. Includes variants
+-- for parallel computation.
+module DistanceTransform.Internal.Indexer where
+import Control.Monad (foldM_)
+import Control.Concurrent (forkIO, getNumCapabilities, 
+                           newEmptyMVar, putMVar, takeMVar)
+import Data.Maybe (fromMaybe)
+
+-- | We use a zipper on list to walk over dimensions of an array.
+data Zipper a = Zip [a] a [a]
+
+-- | Create a 'Zipper' from a non-empty list, with the cursor on the
+-- leftmost element.
+toZipper :: a -> [a] -> Zipper a
+toZipper = Zip []
+
+-- | Create a 'Zipper' from a non-empty list, with the cursor on the
+-- leftmost element. An exception is thrown if the given list is
+-- empty.
+unsafeToZipper :: [a] -> Zipper a
+unsafeToZipper [] = error "A comonad can't be empty!"
+unsafeToZipper (x:xs) = Zip [] x xs
+
+-- | Convert a 'Zipper' to a list.
+fromZipper :: Zipper a -> [a]
+fromZipper (Zip l x r) = reverse l ++ x : r
+
+-- | Move a 'Zipper' to the left.
+left :: Zipper a -> Maybe (Zipper a)
+left (Zip [] _ _) = Nothing
+left (Zip (l:ls) x r) = Just $ Zip ls l (x:r)
+
+unsafeLeft :: Zipper a -> Zipper a
+unsafeLeft z = fromMaybe z $ left z
+
+right :: Zipper a -> Maybe (Zipper a)
+right (Zip _ _ []) = Nothing
+right (Zip ls x (r:rs)) = Just $ Zip (x:ls) r rs
+
+-- | Comonadic coreturn: produce the value a 'Zipper' is currently
+-- focused upon.
+focus :: Zipper a -> a
+focus (Zip _ x _) = x
+
+-- | Slide a 'Zipper' over until focused on its rightmost element.
+rightmost :: Zipper a -> Zipper a
+rightmost z@(Zip _ _ []) = z
+rightmost (Zip ls x (r:rs)) = rightmost $ Zip (x:ls) r rs
+
+zipSum, zipStride, zipStep :: Num a => Zipper a -> a
+-- | Since we are using 'Zipper's to track the size of
+-- multidemensional arrays, the sum of all zipper elements gives the
+-- size of the entire array.
+zipSum = sum . fromZipper
+
+-- | Computes the stride between rows at the currently focused
+-- dimension. This involves stepping over the rest of the current row
+-- along all nested dimensions.
+zipStride (Zip _ x rs) = product $ x:rs
+
+-- | Computes the step between consective elements at the currently
+-- focused dimension. This involves stepping over all nested
+-- dimensions.
+zipStep (Zip _ _ rs) = product rs
+
+-- Each inner loop is stateful.
+zipFoldM :: Monad m => Zipper Int -> (a -> Int -> m a) -> a -> [Int] -> m ()
+zipFoldM (Zip ls x rs) f z indices = gol 0 (reverse ls)
+  where innerDimStride = x * product rs
+        gol offset [] = gor offset rs
+        gol offset (d:ds) = mapM_ (\i -> gol (offset + i*stride) ds) [0..d-1]
+          where stride = product ds * innerDimStride
+        gor offset [] = foldM_ (\s i -> f s (offset + i*stride)) z indices
+          where stride = product rs
+        gor offset (d:ds) = mapM_ (\i -> gor (offset + i*stride) ds) [0..d-1]
+          where stride = product ds
+{-# INLINE zipFoldM #-}
+
+parChunkMapM_ :: (a -> IO ()) -> [a] -> IO ()
+parChunkMapM_ f xs0 = do caps <- getNumCapabilities
+                         let sz = length xs0 `quot` caps
+                         let chunk ts [] = sequence_ ts
+                             chunk ts xs = let (c,xs') = splitAt sz xs
+                                           in do m <- newEmptyMVar
+                                                 _ <- forkIO $ mapM_ f c >> 
+                                                               putMVar m ()
+                                                 chunk (takeMVar m:ts) xs'
+                         chunk [] xs0
+
+parZipFoldM :: Zipper Int -> (a -> Int -> IO a) -> a -> [Int] -> IO ()
+parZipFoldM (Zip ls x rs) f z indices = golPar $ reverse ls
+  where innerDimStride = x * product rs
+        golPar [] = case rs of
+                      [] -> gor 0 []
+                      r:rs' -> let stride = product rs'
+                               in parChunkMapM_ (\i -> gor (i*stride) rs')
+                                                [0..r-1]
+        golPar (d:ds) = parChunkMapM_ (\i -> gol (i*stride) ds) [0..d-1]
+          where stride = product ds * innerDimStride
+        gol offset [] = gor offset rs
+        gol offset (d:ds) = mapM_ (\i -> gol (offset + i*stride) ds) [0..d-1]
+          where stride = product ds * innerDimStride
+        gor offset [] = foldM_ (\s i -> f s (offset + i*stride)) z indices
+          where stride = product rs
+        gor offset (d:ds) = mapM_ (\i -> gor (offset + i*stride) ds) [0..d-1]
+          where stride = product ds
+{-# INLINE parZipFoldM #-}
+
+zipMapM_ :: Monad m => Zipper Int -> (Int -> m ()) -> [Int] -> m ()
+zipMapM_ z f = zipFoldM z (const f) ()
+{-# INLINE zipMapM_ #-}
+
+-- Give a function an offset to the start of its indices and the step
+-- between indices. This lets you walk along *any* dimension within a
+-- packed array. The idea is to let the caller do whatever the heck
+-- they want in that traversal.
+zipFoldMAsYouDo :: Monad m => Zipper Int -> (Int -> Int -> m ()) -> m ()
+zipFoldMAsYouDo z f = zipFoldM z auxOffset Nothing [0,1]
+  where auxOffset Nothing offset = return $ Just offset
+        auxOffset (Just offset) step' = f offset (step' - offset) >>
+                                        return Nothing
+{-# INLINE zipFoldMAsYouDo #-}
+
+-- Give a function an offset to the start of its indices and the step
+-- between indices. This lets you walk along *any* dimension within a
+-- packed array. The idea is to let the caller do whatever the heck
+-- they want in that traversal.
+parZipFoldMAsYouDo :: Zipper Int -> (Int -> Int -> IO ()) -> IO ()
+parZipFoldMAsYouDo z f = parZipFoldM z auxOffset Nothing [0,1]
+  where auxOffset Nothing offset = return $ Just offset
+        auxOffset (Just offset) step' = f offset (step' - offset) >>
+                                        return Nothing
+{-# INLINE parZipFoldMAsYouDo #-}
diff --git a/src/tests/Main.hs b/src/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/src/tests/Main.hs
@@ -0,0 +1,50 @@
+-- | Uhit test that seeds a 3D grid with a few points, computes the
+-- Euclidean distance transform of that grid, then checks a few points
+-- to see if the distance transformed grid agrees with an exhaustive
+-- nearest-neighbor search.
+module Main (main) where
+import qualified Data.Vector.Unboxed as V
+import qualified Data.Vector.Unboxed.Mutable as VM
+import Test.Framework (defaultMain)
+import Test.Framework.Providers.HUnit (testCase)
+import Test.HUnit (assert)
+import DistanceTransform.Euclidean
+
+testRes :: Int
+testRes = 64
+
+-- A 3D point.
+data Point = Point !Int !Int !Int deriving Show
+
+pointToI :: Point -> Int
+pointToI (Point x y z) = z * testRes * testRes + y * testRes + x
+
+distance :: Point -> Point -> Float
+distance (Point x1 y1 z1) (Point x2 y2 z2) = sqrt . fromIntegral $ 
+                                             dx*dx + dy*dy + dz*dz
+  where dx = x2 - x1
+        dy = y2 - y1
+        dz = z2 - z1
+
+mkGrid :: [Point] -> V.Vector Int
+mkGrid pts = V.create $ do v <- VM.replicate (testRes^(3::Int)) 1
+                           mapM_ (flip (VM.write v) 0 . pointToI) pts
+                           return v
+
+main :: IO ()
+main = defaultMain $ map (testPoint g1) probes ++ map (testPoint g2) probes
+  where probes = [ Point 48 32 32
+                 , Point 32 54 35
+                 , Point 0 62 54
+                 , Point 35 35 35 ]
+        pts = Point mid mid mid : 
+              [Point x y z | x <- [0,hi], y <- [0,hi], z <- [0,hi]]
+        mid = testRes `quot` 2
+        hi = testRes - 1
+        rawGrid = mkGrid pts
+        g1 = edt (replicate 3 testRes) rawGrid
+        g2 = edtPar (replicate 3 testRes) $ rawGrid
+        testPoint g probe = let x = minimum $ map (distance probe) pts
+                                y = g V.! pointToI probe
+                            in testCase ("Probing "++show probe)
+                                        (assert (abs (x - y) < 0.0001))
diff --git a/src/tests/TestPar.hs b/src/tests/TestPar.hs
new file mode 100644
--- /dev/null
+++ b/src/tests/TestPar.hs
@@ -0,0 +1,25 @@
+module Main (main) where
+import Criterion.Main
+import qualified Data.Vector.Storable as V
+import qualified Data.Vector.Storable.Mutable as VM
+import qualified Data.Vector.Unboxed as U
+import DistanceTransform.Euclidean
+
+testRes :: Int
+testRes = 64
+
+-- A cube of ones with a zero at the center.
+testData :: V.Vector Int
+testData = V.create $ do v <- VM.replicate (testRes^3) 1
+                         VM.write v (4*testRes^2+4*testRes+4) 0
+                         return v
+
+main = do putStr "I am sane, true or false? "
+          print (edt' dims testData == edtPar' dims testData)
+          defaultMain [ bench "serial" $ whnf (edt' dims) testData
+                      , bench "parallel" $ whnf (edtPar' dims) testData ]
+  where dims = replicate 3 testRes
+        edt' :: [Int] -> V.Vector Int -> V.Vector Float
+        edt' = edt
+        edtPar' :: [Int] -> V.Vector Int -> V.Vector Float
+        edtPar' = edtPar