diff options
Diffstat (limited to 'packages/base/src/Internal/Convolution.hs')
-rw-r--r-- | packages/base/src/Internal/Convolution.hs | 155 |
1 files changed, 155 insertions, 0 deletions
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 @@ | |||
1 | {-# LANGUAGE FlexibleContexts #-} | ||
2 | ----------------------------------------------------------------------------- | ||
3 | {- | | ||
4 | Module : Internal.Convolution | ||
5 | Copyright : (c) Alberto Ruiz 2012 | ||
6 | License : BSD3 | ||
7 | Maintainer : Alberto Ruiz | ||
8 | Stability : provisional | ||
9 | |||
10 | -} | ||
11 | ----------------------------------------------------------------------------- | ||
12 | {-# OPTIONS_HADDOCK hide #-} | ||
13 | |||
14 | module Internal.Convolution( | ||
15 | corr, conv, corrMin, | ||
16 | corr2, conv2, separable | ||
17 | ) where | ||
18 | |||
19 | import qualified Data.Vector.Storable as SV | ||
20 | import Internal.Vector | ||
21 | import Internal.Matrix hiding (mat) | ||
22 | import Internal.Numeric | ||
23 | import Internal.Element | ||
24 | import Internal.Conversion | ||
25 | import Internal.Container | ||
26 | |||
27 | |||
28 | vectSS :: Element t => Int -> Vector t -> Matrix t | ||
29 | vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ] | ||
30 | |||
31 | |||
32 | corr | ||
33 | :: (Container Vector t, Product t) | ||
34 | => Vector t -- ^ kernel | ||
35 | -> Vector t -- ^ source | ||
36 | -> Vector t | ||
37 | {- ^ correlation | ||
38 | |||
39 | >>> corr (fromList[1,2,3]) (fromList [1..10]) | ||
40 | fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] | ||
41 | |||
42 | -} | ||
43 | corr ker v | ||
44 | | dim ker == 0 = konst 0 (dim v) | ||
45 | | dim ker <= dim v = vectSS (dim ker) v <> ker | ||
46 | | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")" | ||
47 | |||
48 | |||
49 | conv :: (Container Vector t, Product t, Num t) => Vector t -> Vector t -> Vector t | ||
50 | {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) | ||
51 | |||
52 | >>> conv (fromList[1,1]) (fromList [-1,1]) | ||
53 | fromList [-1.0,0.0,1.0] | ||
54 | |||
55 | -} | ||
56 | conv ker v | ||
57 | | dim ker == 0 = konst 0 (dim v) | ||
58 | | otherwise = corr ker' v' | ||
59 | where | ||
60 | ker' = SV.reverse ker | ||
61 | v' = vjoin [z,v,z] | ||
62 | z = konst 0 (dim ker -1) | ||
63 | |||
64 | corrMin :: (Container Vector t, RealElement t, Product t) | ||
65 | => Vector t | ||
66 | -> Vector t | ||
67 | -> Vector t | ||
68 | -- ^ similar to 'corr', using 'min' instead of (*) | ||
69 | corrMin ker v | ||
70 | | dim ker == 0 = error "corrMin: empty kernel" | ||
71 | | otherwise = minEvery ss (asRow ker) <> ones | ||
72 | where | ||
73 | minEvery a b = cond a b a a b | ||
74 | ss = vectSS (dim ker) v | ||
75 | ones = konst 1 (dim ker) | ||
76 | |||
77 | |||
78 | |||
79 | matSS :: Element t => Int -> Matrix t -> [Matrix t] | ||
80 | matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ] | ||
81 | where | ||
82 | v = flatten m | ||
83 | c = cols m | ||
84 | r = rows m | ||
85 | n = dr*c | ||
86 | |||
87 | |||
88 | {- | 2D correlation (without padding) | ||
89 | |||
90 | >>> disp 5 $ corr2 (konst 1 (3,3)) (ident 10 :: Matrix Double) | ||
91 | 8x8 | ||
92 | 3 2 1 0 0 0 0 0 | ||
93 | 2 3 2 1 0 0 0 0 | ||
94 | 1 2 3 2 1 0 0 0 | ||
95 | 0 1 2 3 2 1 0 0 | ||
96 | 0 0 1 2 3 2 1 0 | ||
97 | 0 0 0 1 2 3 2 1 | ||
98 | 0 0 0 0 1 2 3 2 | ||
99 | 0 0 0 0 0 1 2 3 | ||
100 | |||
101 | -} | ||
102 | corr2 :: Product a => Matrix a -> Matrix a -> Matrix a | ||
103 | corr2 ker mat = dims | ||
104 | . concatMap (map (udot ker' . flatten) . matSS c . trans) | ||
105 | . matSS r $ mat | ||
106 | where | ||
107 | r = rows ker | ||
108 | c = cols ker | ||
109 | ker' = flatten (trans ker) | ||
110 | rr = rows mat - r + 1 | ||
111 | rc = cols mat - c + 1 | ||
112 | dims | rr > 0 && rc > 0 = (rr >< rc) | ||
113 | | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")" | ||
114 | sz m = show (rows m)++"x"++show (cols m) | ||
115 | -- TODO check empty kernel | ||
116 | |||
117 | {- | 2D convolution | ||
118 | |||
119 | >>> disp 5 $ conv2 (konst 1 (3,3)) (ident 10 :: Matrix Double) | ||
120 | 12x12 | ||
121 | 1 1 1 0 0 0 0 0 0 0 0 0 | ||
122 | 1 2 2 1 0 0 0 0 0 0 0 0 | ||
123 | 1 2 3 2 1 0 0 0 0 0 0 0 | ||
124 | 0 1 2 3 2 1 0 0 0 0 0 0 | ||
125 | 0 0 1 2 3 2 1 0 0 0 0 0 | ||
126 | 0 0 0 1 2 3 2 1 0 0 0 0 | ||
127 | 0 0 0 0 1 2 3 2 1 0 0 0 | ||
128 | 0 0 0 0 0 1 2 3 2 1 0 0 | ||
129 | 0 0 0 0 0 0 1 2 3 2 1 0 | ||
130 | 0 0 0 0 0 0 0 1 2 3 2 1 | ||
131 | 0 0 0 0 0 0 0 0 1 2 2 1 | ||
132 | 0 0 0 0 0 0 0 0 0 1 1 1 | ||
133 | |||
134 | -} | ||
135 | conv2 | ||
136 | :: (Num (Matrix a), Product a, Container Vector a) | ||
137 | => Matrix a -- ^ kernel | ||
138 | -> Matrix a -> Matrix a | ||
139 | conv2 k m | ||
140 | | empty = konst 0 (rows m + r -1, cols m + c -1) | ||
141 | | otherwise = corr2 (fliprl . flipud $ k) padded | ||
142 | where | ||
143 | padded = fromBlocks [[z,0,0] | ||
144 | ,[0,m,0] | ||
145 | ,[0,0,z]] | ||
146 | r = rows k | ||
147 | c = cols k | ||
148 | z = konst 0 (r-1,c-1) | ||
149 | empty = r == 0 || c == 0 | ||
150 | |||
151 | |||
152 | separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t | ||
153 | -- ^ matrix computation implemented as separated vector operations by rows and columns. | ||
154 | separable f = fromColumns . map f . toColumns . fromRows . map f . toRows | ||
155 | |||