From 8b093ecca2b4e200ff191b84cb0b56a12312867b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 5 Jun 2015 16:38:09 +0200 Subject: move convolution --- packages/base/src/Internal/Convolution.hs | 155 ++++++++++++++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 packages/base/src/Internal/Convolution.hs (limited to 'packages/base/src/Internal') diff --git a/packages/base/src/Internal/Convolution.hs b/packages/base/src/Internal/Convolution.hs new file mode 100644 index 0000000..1a70011 --- /dev/null +++ b/packages/base/src/Internal/Convolution.hs @@ -0,0 +1,155 @@ +{-# LANGUAGE FlexibleContexts #-} +----------------------------------------------------------------------------- +{- | +Module : Internal.Convolution +Copyright : (c) Alberto Ruiz 2012 +License : BSD3 +Maintainer : Alberto Ruiz +Stability : provisional + +-} +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Internal.Convolution( + corr, conv, corrMin, + corr2, conv2, separable +) where + +import qualified Data.Vector.Storable as SV +import Internal.Vector +import Internal.Matrix hiding (mat) +import Internal.Numeric +import Internal.Element +import Internal.Conversion +import Internal.Container + + +vectSS :: Element t => Int -> Vector t -> Matrix t +vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ] + + +corr + :: (Container Vector t, Product t) + => Vector t -- ^ kernel + -> Vector t -- ^ source + -> Vector t +{- ^ correlation + +>>> corr (fromList[1,2,3]) (fromList [1..10]) +fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] + +-} +corr ker v + | dim ker == 0 = konst 0 (dim v) + | dim ker <= dim v = vectSS (dim ker) v <> ker + | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")" + + +conv :: (Container Vector t, Product t, Num t) => Vector t -> Vector t -> Vector t +{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) + +>>> conv (fromList[1,1]) (fromList [-1,1]) +fromList [-1.0,0.0,1.0] + +-} +conv ker v + | dim ker == 0 = konst 0 (dim v) + | otherwise = corr ker' v' + where + ker' = SV.reverse ker + v' = vjoin [z,v,z] + z = konst 0 (dim ker -1) + +corrMin :: (Container Vector t, RealElement t, Product t) + => Vector t + -> Vector t + -> Vector t +-- ^ similar to 'corr', using 'min' instead of (*) +corrMin ker v + | dim ker == 0 = error "corrMin: empty kernel" + | otherwise = minEvery ss (asRow ker) <> ones + where + minEvery a b = cond a b a a b + ss = vectSS (dim ker) v + ones = konst 1 (dim ker) + + + +matSS :: Element t => Int -> Matrix t -> [Matrix t] +matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ] + where + v = flatten m + c = cols m + r = rows m + n = dr*c + + +{- | 2D correlation (without padding) + +>>> disp 5 $ corr2 (konst 1 (3,3)) (ident 10 :: Matrix Double) +8x8 +3 2 1 0 0 0 0 0 +2 3 2 1 0 0 0 0 +1 2 3 2 1 0 0 0 +0 1 2 3 2 1 0 0 +0 0 1 2 3 2 1 0 +0 0 0 1 2 3 2 1 +0 0 0 0 1 2 3 2 +0 0 0 0 0 1 2 3 + +-} +corr2 :: Product a => Matrix a -> Matrix a -> Matrix a +corr2 ker mat = dims + . concatMap (map (udot ker' . flatten) . matSS c . trans) + . matSS r $ mat + where + r = rows ker + c = cols ker + ker' = flatten (trans ker) + rr = rows mat - r + 1 + rc = cols mat - c + 1 + dims | rr > 0 && rc > 0 = (rr >< rc) + | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")" + sz m = show (rows m)++"x"++show (cols m) +-- TODO check empty kernel + +{- | 2D convolution + +>>> disp 5 $ conv2 (konst 1 (3,3)) (ident 10 :: Matrix Double) +12x12 +1 1 1 0 0 0 0 0 0 0 0 0 +1 2 2 1 0 0 0 0 0 0 0 0 +1 2 3 2 1 0 0 0 0 0 0 0 +0 1 2 3 2 1 0 0 0 0 0 0 +0 0 1 2 3 2 1 0 0 0 0 0 +0 0 0 1 2 3 2 1 0 0 0 0 +0 0 0 0 1 2 3 2 1 0 0 0 +0 0 0 0 0 1 2 3 2 1 0 0 +0 0 0 0 0 0 1 2 3 2 1 0 +0 0 0 0 0 0 0 1 2 3 2 1 +0 0 0 0 0 0 0 0 1 2 2 1 +0 0 0 0 0 0 0 0 0 1 1 1 + +-} +conv2 + :: (Num (Matrix a), Product a, Container Vector a) + => Matrix a -- ^ kernel + -> Matrix a -> Matrix a +conv2 k m + | empty = konst 0 (rows m + r -1, cols m + c -1) + | otherwise = corr2 (fliprl . flipud $ k) padded + where + padded = fromBlocks [[z,0,0] + ,[0,m,0] + ,[0,0,z]] + r = rows k + c = cols k + z = konst 0 (r-1,c-1) + empty = r == 0 || c == 0 + + +separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t +-- ^ matrix computation implemented as separated vector operations by rows and columns. +separable f = fromColumns . map f . toColumns . fromRows . map f . toRows + -- cgit v1.2.3