MeanShift (empty) → 0.1

raw patch · 6 files changed

+253/−0 lines, 6 filesdep +basedep +vectorsetup-changedbinary-added

Dependencies added: base, vector

Files

• Examples/MeanShiftFilter.hs +56/−0
• Examples/meanshift-filter.png binary
• LICENSE +30/−0
• MeanShift.cabal +27/−0
• Setup.hs +2/−0
• src/Math/Meanshift.hs +138/−0

@@ -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 ()

binary file changed (absent → 286577 bytes)

@@ -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.

@@ -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 

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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+