diff --git a/Examples/MeanShiftFilter.hs b/Examples/MeanShiftFilter.hs
new file mode 100644
--- /dev/null
+++ b/Examples/MeanShiftFilter.hs
@@ -0,0 +1,56 @@
+{-#LANGUAGE BangPatterns#-}
+module Main where
+import Data.List hiding (sum)
+import qualified Data.Vector.Unboxed as V
+import Prelude hiding (sum)
+import Control.Arrow
+import MeanShift
+import CV.Image
+import CV.Pixelwise
+import CV.Transforms
+import System.Environment
+import Data.Maybe
+import Control.Monad
+
+-- | Extract a pixel from an image.
+sample :: (Int,Int) -> Image RGB D32 -> (Int,Int) -> Vector
+sample (w,h) image (!x,!y)
+   | x>0 && y>0 && x < w && y < h = let (r,g,b) = getPixel (x,y) image
+                                    in V.fromList [fi x, fi y, colorScale*realToFrac r, colorScale*realToFrac g, colorScale*realToFrac b]
+   | otherwise = V.fromList [fi x, fi y, 0,0,0]
+
+-- | Extract color and coordinate information from the sample (see above.) 
+coord v = let (!x):(!y):_ = V.toList v in (round x,round y)
+color v = let _:_:r:g:b:_ = V.toList v in (realToFrac (b/colorScale), realToFrac (g/colorScale), realToFrac (r/colorScale))
+
+-- | Weighting term to balance the effect between spatial and color domains.
+colorScale = 50
+
+main :: IO ()
+main = do
+   [fn] <- getArgs
+   image <- readFromFile fn
+   let -- Reading the test-data. You might want to scale it to smaller size to avoid waiting too much. 
+       testData :: Image RGB D32
+       testData =  image 
+       (w,h) = getSize testData
+
+       -- A windowing function that is used to extract image patches. 
+       window :: Window
+       window (x,s) = [sample (w,h) testData (u+j,v+i)
+                      | let size  = round s
+                      , i <- [-size,-size+1..size]
+                      , j <- [-size,-size+1..size]
+                      , let (u,v) = coord x
+                      ]
+
+       -- Application of meanshift to the image data.
+       shift   = meanShiftWindow 5 window 3
+       process = last . take 12 . fixedPointE 0.1 shift
+       
+       -- The resulting image. This will be evaluated in parallel if you compile with `-rtsopt -threaded`-options and
+       -- execute the resulting program with `-RTS -N` flags
+       r :: Image RGB D32
+       r = toImagePar 8 $ MkP (getSize testData) (color . process . sample (w,h) testData)
+   saveImage "filtered.png" r
+   return ()
diff --git a/Examples/meanshift-filter.png b/Examples/meanshift-filter.png
new file mode 100644
Binary files /dev/null and b/Examples/meanshift-filter.png differ
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2012, Ville Tirrone
+
+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 Ville Tirrone 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/MeanShift.cabal b/MeanShift.cabal
new file mode 100644
--- /dev/null
+++ b/MeanShift.cabal
@@ -0,0 +1,27 @@
+Name:                MeanShift
+Version:             0.1
+Synopsis:            Mean shift algorithm
+Description:         Mean shift is a general, non-parametric feature-space analysis tool. It can be used
+                     for clustering, segmentation, filtering, object tracking, and even optimization. This package aims to
+                     provide a basic, easy to use version of the method.
+License:             BSD3
+License-file:        LICENSE
+Author:              Ville Tirronen
+Maintainer:          aleator@gmail.com
+Category:            Math
+Build-type:          Simple
+Cabal-version:       >=1.6
+
+Extra-source-files: Examples/MeanShiftFilter.hs
+                        ,Examples/meanshift-filter.png
+
+source-repository head
+    type: git
+    location: https://github.com/aleator/Meanshift
+
+Library
+    hs-source-dirs: src
+    ghc-options: -O2
+    exposed-modules: Math.Meanshift
+    build-depends:   base > 4 && < 5
+                    ,vector > 0.9 && < 0.10 
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/Math/Meanshift.hs b/src/Math/Meanshift.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Meanshift.hs
@@ -0,0 +1,138 @@
+{-#LANGUAGE BangPatterns, ParallelListComp#-}
+-- | This module presents a basic version of the meanshift algorithm for
+-- feature-space analysis. Mean shifting is an iterative process with
+-- fixed points that correspond to 
+-- modes of kernel density estimate performed
+-- with the same bandwidth (first parameter). This 
+-- can be used to, for example, to partition the data by
+-- determining which fixed point each of the samples belongs to.
+-- 
+-- Usage example:
+--  > fixedPointE 0.001 (meanShift 0.1 points) (V.fromList [1,1,1])
+--
+--  More examples can be found in the Examples directory of this package.
+
+module Math.Meanshift 
+    (
+     -- * Basic Meanshift routines
+      meanShift,meanShiftWindow
+     -- * Auxiliary functions for iterating the meanshift steps.
+     ,fixedPoint, fixedPointE
+     -- * Types
+     ,Window,Support
+     -- * (multidimensional) Kernel Density Estimates
+     ,kde
+    ) where
+import qualified Data.Vector.Unboxed as V
+import Data.List hiding (sum)
+import qualified Data.List as L
+import Prelude hiding (sum)
+
+-- The project cabal file <url:../../MeanShift.cabal>
+
+-- | Euclidian norm
+norm² :: Vector -> Double
+norm² = V.sum . V.map (**2)
+
+
+-- | One dimensional normal kernel and its derivative.
+normalKernel,normalKernel' :: Double -> Double
+normalKernel x = exp(-0.5 * x)
+normalKernel' x = -2 *  exp(-0.5 * x)
+
+-- | Kernel density estimate of given points. Uses a normal kernel.
+kde :: Double -> [Vector] -> (Vector -> Double)
+kde h vs x = (1 / (n*((2*π)**(d/2))*(h**d)))
+             * (L.sum $ map (\xi -> normalKernel (norm² ((x ^- xi) ./ h))) vs)
+   where
+      n = fi . length $ vs
+      d = fi . V.length . head $ vs
+
+
+-- | Calculate the Mean shift for a point in a dataset. This is
+-- efficient only when we cannot make an a priori estimate on which
+-- points contribute to the mean shift at given location.
+--
+meanShift :: Double -> [Vector] -> (Vector -> Vector)
+meanShift h vs x = sumW d vs dists (1/V.sum dists) 
+   where
+    d = V.length (head vs)
+    dists = V.fromList $ map (\xi -> normalKernel' $ distPerH x xi) vs
+    distPerH :: Vector -> Vector  -> Double
+    distPerH !a !b = V.sum (V.zipWith (\u v -> ((u-v) / h)^(2::Int)) a b)
+
+type Window = Support -> [Vector]
+type Support = (Vector,Double)
+
+-- | Mean shift with a windowing function. Performing mean shift is more
+--   efficient if we can index and calculate only those points that are in
+--   the support of our kernel.
+{-#INLINEABLE meanShiftWindow#-}
+meanShiftWindow :: Int -> Window -> Double -> (Vector -> Vector)
+meanShiftWindow d window h x
+    = sumW d w dists (1/V.sum dists) 
+   where
+    dists = V.fromList $ map (\xi -> normalKernel' $ distPerH x xi) w
+    w = window (x,h*2) -- TODO: Think this through
+    distPerH :: Vector -> Vector -> Double
+    distPerH !a !b = V.sum  (V.zipWith (\u v -> ((u-v) / h)^(2::Int)) a b)
+
+-- | Find a path to the fixed point of a function.
+{-#INLINEABLE fixedPoint#-}
+fixedPoint :: Eq a => (a -> a) -> a -> [a]
+fixedPoint f x = x:let x' = f x in if x'/=x then fixedPoint f x' else [x']
+
+fixedPointE :: Double -> (Vector -> Vector) -> Vector -> [Vector]
+fixedPointE e f x = x:let x' = f x
+                    in if V.sum (V.map abs $ x' ^- x) > e then fixedPointE e f x' else [x']
+
+
+-- * Auxiliary functions, and shorthands
+
+type Vector = V.Vector Double
+
+v :: [Double] -> Vector
+v = V.fromList
+
+{-#INLINE (^+)#-}
+{-#INLINE (^-)#-}
+{-#INLINE (^/)#-}
+(^+),(^-),(^/) :: Vector -> Vector -> Vector
+(^-) = V.zipWith (-)
+(^+) = V.zipWith (+)
+(^/) = V.zipWith (/)
+a .+ b = V.map (+b) a
+(.+),(./),(.*) :: Vector -> Double -> Vector
+a ./ b = V.map (/b) a
+a .* b = V.map (*b) a
+infixl 7 ^/
+infixl 6 ^+ , ^-
+box :: x -> [x]
+box x = [x]
+box2 :: x -> x -> [x]
+box2 x y = [x,y]
+
+sv :: Double -> Vector
+sv = V.singleton
+fs :: Vector -> Double
+fs = V.head
+
+{-#INLINEABLE sumD#-}
+sumD :: Int -> [Vector] -> Vector
+sumD d xs = V.generate d (\i -> L.sum (map (`V.unsafeIndex` i) xs) )
+
+
+{-#INLINEABLE sumW#-}
+sumW :: Int -> [Vector] -> Vector -> Double -> Vector
+sumW d es ws n = V.generate d (\i -> go i 0 es 0)
+   where
+      go i j (x:xs) acc = go i (j+1) xs $ acc + n*(x V.! i)*(ws V.! j)
+      go _ _ []     acc = acc
+
+
+π :: Double
+π = pi
+
+fi :: Int -> Double
+fi = fromIntegral
+