RANSAC (empty) → 0.1.0.0
raw patch · 6 files changed
+268/−0 lines, 6 filesdep +basedep +randomdep +vectorsetup-changed
Dependencies added: base, random, vector
Files
- LICENSE +30/−0
- RANSAC.cabal +33/−0
- Setup.hs +2/−0
- src/Numeric/Ransac.hs +62/−0
- tests/LinearFit.hs +97/−0
- tests/Perf.hs +44/−0
+ LICENSE view
@@ -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.
+ RANSAC.cabal view
@@ -0,0 +1,33 @@+name: RANSAC+version: 0.1.0.0+synopsis: The RANSAC algorithm for parameter estimation.+description: The RANdom SAmple Consensus (RANSAC) algorithm for+ estimating the parameters of a mathematical model+ from a data set. See+ <http://en.wikipedia.org/wiki/RANSAC> for more+ information.+ .+ See @tests/LinearFit.hs@ in the package contents for + an example.+license: BSD3+license-file: LICENSE+author: Anthony Cowley+maintainer: acowley@gmail.com+copyright: (c) Anthony Cowley 2012+category: Math,Numerical+build-type: Simple+cabal-version: >=1.10+extra-source-files: tests/Perf.hs, tests/LinearFit.hs++source-repository head+ type: git+ location: git://github.com/acowley/RANSAC.git++library+ exposed-modules: Numeric.Ransac+ build-depends: base >= 4.6 && < 5, + vector >= 0.10, + random >= 1.0+ hs-source-dirs: src+ default-language: Haskell2010+ ghc-options: -Wall
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Numeric/Ransac.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE BangPatterns #-}+-- | The RANdom SAmple Consensus (RANSAC) algorithm for estimating the+-- parameters of a mathematical model from a data set. See+-- <http://en.wikipedia.org/wiki/RANSAC> for more information.+module Numeric.Ransac (ransac) where+import Control.Applicative+import Control.Monad (replicateM)+import Data.Vector.Generic ((!))+import qualified Data.Vector.Generic as V+import System.Random++--randDistinct :: (Eq a, Random a, Monad m) => Int -> m a -> m [a]+randDistinct :: Int -> IO Int -> IO [Int]+randDistinct n gen = go 0 [] []+ where go !i acc _ | i == n = return acc+ go !i acc [] = replicateM (n-i) gen >>= go i acc+ go !i acc (r:rs) = if r `elem` acc + then go i acc rs + else go (i+1) (r:acc) rs+{- SPECIALIZE randDistinct :: Int -> IO Int -> IO [Int] #-}++untilJust :: Monad m => m (Maybe b) -> m b+untilJust x = go + where go = x >>= maybe go return+{-# INLINE untilJust #-}++-- | @ransac iter sampleSize agreePct fit residual goodFit pts@ draws+-- @iter@ samples of size @sampleSize@ from @pts@. The @fit@ function+-- is used to produce a model from each of these samples. The elements+-- of @pts@ whose residuals pass the @goodFit@ predicate with respect+-- to this model are identified as /inliers/, and used to update the+-- model. The model for which the size of the inliers set is at least+-- @agreePct@ percent of the entire data set and whose error over all+-- points is minimal among all sampled models is returned. If no+-- acceptable model is found (i.e. no model whose inliers were at+-- least @agreePct@ percent of the entire data set), 'Nothing' is+-- returned.+ransac :: (V.Vector v a, V.Vector v d, Num d, Ord d) => + Int -> Int -> Float -> + (v a -> Maybe c) -> (c -> a -> d) -> (d -> Bool) -> + v a -> IO (Maybe (c, v a))+ransac maxIter sampleSize agree fit residual goodFit pts = genModel >>= go 0+ where go i r@(model, bestError, inliers)+ | i == maxIter = if ratioInliers (V.length inliers) < agree+ then return Nothing+ else return (Just (model, inliers))+ | otherwise = do r'@(_, err, inliers') <- genModel+ if ratioInliers (V.length inliers') >= agree+ && err < bestError+ then go (i+1) r'+ else go (i+1) r+ sample = V.fromList . map (pts !) <$> + randDistinct sampleSize ((`rem` n) . abs <$> randomIO)+ genModel = do model <- untilJust (fit <$> sample)+ let !errors = V.map (residual model) pts+ !inliers = V.ifilter (const . goodFit . (errors !)) pts+ Just model' = fit inliers+ err = V.sum $ V.map (residual model') pts+ return (model', err, inliers)+ n = V.length pts+ ratioInliers n' = fromIntegral n' / fromIntegral n+{-# INLINE ransac #-}
+ tests/LinearFit.hs view
@@ -0,0 +1,97 @@+-- | Example use of the RANSAC algorithm to fit a line to some+-- points. We start with points generated by a process defined by the+-- equation of a line in 2D. These points are affected by normally+-- distributed noise, and our data set is further corrupted by a "red+-- herring" cluster of points that we would like to ignore. We use+-- RANSAC to cut through the noise and fit a line to the point data+-- set.+-- +-- The important feature of RANSAC as applied here is that it manages+-- to ignore the spurious (red herring) cluster centered at (0,8).+--+-- The Chart package is used to visualize the data and estimated+-- model.+module Main where+import Control.Applicative+import Control.Lens (view)+import Data.Accessor ((^=))+import Data.Colour (opaque)+import Data.Colour.Names+import qualified Data.Foldable as F+import Data.Random.Normal (normalsIO')+import Data.Vector.Storable (Vector)+import qualified Data.Vector.Storable as V+import Graphics.Rendering.Chart hiding (Vector,Point)+import Linear+import Numeric.Ransac++type Point = V2 Float++-- | Fit a 2D line to a collection of 'Point's.+fitLine :: Vector Point -> Maybe (V2 Float)+fitLine pts = (!* b) <$> inv22 a+ where sx = V.sum $ V.map (view _x) pts+ a = V2 (V2 (V.sum (V.map ((^2).view _x) pts)) sx)+ (V2 sx (fromIntegral (V.length pts)))+ b = V2 (V.sum (V.map F.product pts))+ (V.sum (V.map (view _y) pts))++-- | Compute the error of a 'Point' with respect to a hypothesized+-- linear model.+ptError :: V2 Float -> Point -> Float+ptError (V2 m b) (V2 x y) = sq $ y - (m*x+b)+ where sq x = x * x++-- | Produce a plot of all the points we have to work with. A green+-- dashed line indicates the ground truth linear model, the solid+-- purple line shows the RANSAC model, and the points that are inliers+-- for that model are circled in yellow.+main = do noise <- v2Cast . V.fromList . take (n*2) <$> normalsIO' (0,0.3)+ herring <- V.zipWith V2 + <$> (V.fromList . take 200 <$> normalsIO' (0,0.2))+ <*> (V.fromList . take 200 <$> normalsIO' (8,0.6))+ let pts' = V.zipWith (+) noise pts+ let pts'' = pts' V.++ herring+ res <- ransac 100 2 0.5 fitLine ptError (< 2) pts''+ case res of+ Nothing -> putStrLn "No model found"+ Just (model,inliers) -> + do putStrLn $ "Model "++show model++" with "+++ show (V.length inliers)++" inliers"+ let pp = PlotPoints "data" + (filledCircles 2 (opaque blue))+ (map (toTup . dub) (V.toList pts''))+ ppi = PlotPoints "inliers"+ (hollowCircles 3 2 (opaque yellow))+ (map (toTup . dub) (V.toList inliers))+ lp = PlotLines "truth"+ (dashedLine 3 [10,10] (opaque green))+ [[ toTup $ dub (mkPt 0) + , toTup $ dub (mkPt (n-1)) ]]+ []+ lp' = PlotLines "model"+ (solidLine 5 (opaque purple))+ [[ toTup $ dub (mkPt' model 0)+ , toTup $ dub (mkPt' model (n-1)) ]]+ []+ layout = layout1_title ^="2D Linear Fit"+ $ layout1_background ^= solidFillStyle (opaque white)+ $ layout1_plots ^= [ Left (toPlot pp)+ , Left (toPlot ppi)+ , Left (toPlot lp)+ , Left (toPlot lp') ]+ $ setLayout1Foreground (opaque black)+ $ defaultLayout1+ renderableToPDFFile (toRenderable layout) 600 600 "foo.pdf"+ where n = 1000+ pts = V.generate n mkPt+ mkPt :: Int -> V2 Float+ mkPt i = let x = fromIntegral i / 500+ in V2 x (5*x + 2)+ v2Cast :: Vector Float -> Vector Point+ v2Cast = V.unsafeCast+ toTup (V2 x y) = (x,y)+ dub :: V2 Float -> V2 Double+ dub = fmap realToFrac+ mkPt' (V2 m b) i = let x = fromIntegral i / 500+ in V2 x (x * m + b)
+ tests/Perf.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE BangPatterns #-}+module Main (main) where+import Control.Applicative+import Control.Lens (view)+import Criterion.Main+import Data.Accessor ((^=))+import Data.Colour (opaque)+import Data.Colour.Names+import qualified Data.Foldable as F+import Data.Random.Normal (normalsIO')+import Data.Vector.Storable (Vector)+import qualified Data.Vector.Storable as V+import Linear+import Numeric.Ransac++type Point = V2 Float++fitLine :: Vector Point -> Maybe (V2 Float)+fitLine pts = (!* b) <$> inv22 a+ where sx = V.sum $ V.map (view _x) pts+ a = V2 (V2 (V.sum (V.map ((^2).view _x) pts)) sx)+ (V2 sx (fromIntegral $ V.length pts))+ b = V2 (V.sum (V.map F.product pts))+ (V.sum (V.map (view _y) pts))++ptError :: V2 Float -> Point -> Float+ptError (V2 m b) (V2 x y) = sq $ y - (m*x+b)+ where sq x = x * x++main = do noise <- v2Cast . V.fromList . take (n*2) <$> normalsIO' (0,0.3)+ let !pts' = V.zipWith (+) noise pts+ ran = ransac 100 2 0.6 fitLine ptError (< 1) pts'+ putStr $ "Sanity "+ ran >>= putStrLn . show . fmap fst+ defaultMain [ bench "linear fit" $ fmap (quadrance . fst) <$> ran ]+ where n = 10000+ !pts = V.generate n mkPt+ mkPt :: Int -> V2 Float+ mkPt i = let x = fromIntegral i / 500+ in V2 x (3*x + 2)+ v2Cast :: Vector Float -> Vector Point+ v2Cast = V.unsafeCast++