|
5819
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 1999-2005 Andy Adler |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include <config.h> |
|
|
26 #endif |
|
|
27 |
|
|
28 #include "defun-dld.h" |
|
|
29 #include "error.h" |
|
|
30 #include "oct-obj.h" |
|
|
31 #include "utils.h" |
|
|
32 |
|
|
33 #define MAX(a,b) ((a) > (b) ? (a) : (b)) |
|
|
34 |
|
|
35 enum Shape { SHAPE_FULL, SHAPE_SAME, SHAPE_VALID }; |
|
|
36 |
|
|
37 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
38 extern MArray2<double> |
|
|
39 conv2 (MArray<double>&, MArray<double>&, MArray2<double>&, Shape); |
|
|
40 |
|
|
41 extern MArray2<Complex> |
|
|
42 conv2 (MArray<Complex>&, MArray<Complex>&, MArray2<Complex>&, Shape); |
|
|
43 #endif |
|
|
44 |
|
|
45 template <class T> |
|
|
46 MArray2<T> |
|
|
47 conv2 (MArray<T>& R, MArray<T>& C, MArray2<T>& A, Shape ishape) |
|
|
48 { |
|
|
49 octave_idx_type Rn = R.length (); |
|
|
50 octave_idx_type Cm = C.length (); |
|
|
51 octave_idx_type Am = A.rows (); |
|
|
52 octave_idx_type An = A.columns (); |
|
|
53 |
|
|
54 // Calculate the size of the output matrix: |
|
|
55 // in order to stay Matlab compatible, it is based |
|
|
56 // on the third parameter if it's separable, and the |
|
|
57 // first if it's not |
|
|
58 |
|
|
59 octave_idx_type outM = 0; |
|
|
60 octave_idx_type outN = 0; |
|
|
61 octave_idx_type edgM = 0; |
|
|
62 octave_idx_type edgN = 0; |
|
|
63 |
|
|
64 switch (ishape) |
|
|
65 { |
|
|
66 case SHAPE_FULL: |
|
|
67 outM = Am + Cm - 1; |
|
|
68 outN = An + Rn - 1; |
|
|
69 edgM = Cm - 1; |
|
|
70 edgN = Rn - 1; |
|
|
71 break; |
|
|
72 |
|
|
73 case SHAPE_SAME: |
|
|
74 outM = Am; |
|
|
75 outN = An; |
|
|
76 // Follow the Matlab convention (ie + instead of -) |
|
|
77 edgM = (Cm - 1) /2; |
|
|
78 edgN = (Rn - 1) /2; |
|
|
79 break; |
|
|
80 |
|
|
81 case SHAPE_VALID: |
|
|
82 outM = Am - Cm + 1; |
|
|
83 outN = An - Rn + 1; |
|
|
84 if (outM < 0) |
|
|
85 outM = 0; |
|
|
86 if (outN < 0) |
|
|
87 outN = 0; |
|
|
88 edgM = edgN = 0; |
|
|
89 break; |
|
|
90 |
|
|
91 default: |
|
|
92 error ("conv2: invalid value of parameter ishape"); |
|
|
93 } |
|
|
94 |
|
|
95 MArray2<T> O (outM, outN); |
|
|
96 |
|
|
97 // X accumulates the 1-D conv for each row, before calculating |
|
|
98 // the convolution in the other direction |
|
|
99 // There is no efficiency advantage to doing it in either direction |
|
|
100 // first |
|
|
101 |
|
|
102 MArray<T> X (An); |
|
|
103 |
|
|
104 for (octave_idx_type oi = 0; oi < outM; oi++) |
|
|
105 { |
|
|
106 for (octave_idx_type oj = 0; oj < An; oj++) |
|
|
107 { |
|
|
108 T sum = 0; |
|
|
109 |
|
|
110 octave_idx_type ci = Cm - 1 - MAX(0, edgM-oi); |
|
|
111 octave_idx_type ai = MAX(0, oi-edgM); |
|
|
112 const T* Ad = A.data() + ai + Am*oj; |
|
|
113 const T* Cd = C.data() + ci; |
|
|
114 for ( ; ci >= 0 && ai < Am; ci--, Cd--, ai++, Ad++) |
|
|
115 sum += (*Ad) * (*Cd); |
|
|
116 |
|
|
117 X(oj) = sum; |
|
|
118 } |
|
|
119 |
|
|
120 for (octave_idx_type oj = 0; oj < outN; oj++) |
|
|
121 { |
|
|
122 T sum = 0; |
|
|
123 |
|
|
124 octave_idx_type rj = Rn - 1 - MAX(0, edgN-oj); |
|
|
125 octave_idx_type aj = MAX(0, oj-edgN) ; |
|
|
126 const T* Xd = X.data() + aj; |
|
|
127 const T* Rd = R.data() + rj; |
|
|
128 |
|
|
129 for ( ; rj >= 0 && aj < An; rj--, Rd--, aj++, Xd++) |
|
|
130 sum += (*Xd) * (*Rd); |
|
|
131 |
|
|
132 O(oi,oj)= sum; |
|
|
133 } |
|
|
134 } |
|
|
135 |
|
|
136 return O; |
|
|
137 } |
|
|
138 |
|
|
139 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
140 extern MArray2<double> |
|
|
141 conv2 (MArray2<double>&, MArray2<double>&, Shape); |
|
|
142 |
|
|
143 extern MArray2<Complex> |
|
|
144 conv2 (MArray2<Complex>&, MArray2<Complex>&, Shape); |
|
|
145 #endif |
|
|
146 |
|
|
147 template <class T> |
|
|
148 MArray2<T> |
|
|
149 conv2 (MArray2<T>&A, MArray2<T>&B, Shape ishape) |
|
|
150 { |
|
|
151 // Convolution works fastest if we choose the A matrix to be |
|
|
152 // the largest. |
|
|
153 |
|
|
154 // Here we calculate the size of the output matrix, |
|
|
155 // in order to stay Matlab compatible, it is based |
|
|
156 // on the third parameter if it's separable, and the |
|
|
157 // first if it's not |
|
|
158 |
|
|
159 // NOTE in order to be Matlab compatible, we give argueably |
|
|
160 // wrong sizes for 'valid' if the smallest matrix is first |
|
|
161 |
|
|
162 octave_idx_type Am = A.rows (); |
|
|
163 octave_idx_type An = A.columns (); |
|
|
164 octave_idx_type Bm = B.rows (); |
|
|
165 octave_idx_type Bn = B.columns (); |
|
|
166 |
|
|
167 octave_idx_type outM = 0; |
|
|
168 octave_idx_type outN = 0; |
|
|
169 octave_idx_type edgM = 0; |
|
|
170 octave_idx_type edgN = 0; |
|
|
171 |
|
|
172 switch (ishape) |
|
|
173 { |
|
|
174 case SHAPE_FULL: |
|
|
175 outM = Am + Bm - 1; |
|
|
176 outN = An + Bn - 1; |
|
|
177 edgM = Bm - 1; |
|
|
178 edgN = Bn - 1; |
|
|
179 break; |
|
|
180 |
|
|
181 case SHAPE_SAME: |
|
|
182 outM = Am; |
|
|
183 outN = An; |
|
|
184 edgM = (Bm - 1) /2; |
|
|
185 edgN = (Bn - 1) /2; |
|
|
186 break; |
|
|
187 |
|
|
188 case SHAPE_VALID: |
|
|
189 outM = Am - Bm + 1; |
|
|
190 outN = An - Bn + 1; |
|
|
191 if (outM < 0) |
|
|
192 outM = 0; |
|
|
193 if (outN < 0) |
|
|
194 outN = 0; |
|
|
195 edgM = edgN = 0; |
|
|
196 break; |
|
|
197 } |
|
|
198 |
|
|
199 MArray2<T> O (outM, outN); |
|
|
200 |
|
|
201 for (octave_idx_type oi = 0; oi < outM; oi++) |
|
|
202 { |
|
|
203 for (octave_idx_type oj = 0; oj < outN; oj++) |
|
|
204 { |
|
|
205 T sum = 0; |
|
|
206 |
|
|
207 for (octave_idx_type bj = Bn - 1 - MAX (0, edgN-oj), aj= MAX (0, oj-edgN); |
|
|
208 bj >= 0 && aj < An; bj--, aj++) |
|
|
209 { |
|
|
210 octave_idx_type bi = Bm - 1 - MAX (0, edgM-oi); |
|
|
211 octave_idx_type ai = MAX (0, oi-edgM); |
|
|
212 const T* Ad = A.data () + ai + Am*aj; |
|
|
213 const T* Bd = B.data () + bi + Bm*bj; |
|
|
214 |
|
|
215 for ( ; bi >= 0 && ai < Am; bi--, Bd--, ai++, Ad++) |
|
|
216 { |
|
|
217 sum += (*Ad) * (*Bd); |
|
|
218 // Comment: it seems to be 2.5 x faster than this: |
|
|
219 // sum+= A(ai,aj) * B(bi,bj); |
|
|
220 } |
|
|
221 } |
|
|
222 |
|
|
223 O(oi,oj) = sum; |
|
|
224 } |
|
|
225 } |
|
|
226 |
|
|
227 return O; |
|
|
228 } |
|
|
229 |
|
|
230 /* |
|
|
231 %!test |
|
|
232 %! b = [0,1,2,3;1,8,12,12;4,20,24,21;7,22,25,18]; |
|
|
233 %! assert(conv2([0,1;1,2],[1,2,3;4,5,6;7,8,9]),b); |
|
|
234 */ |
|
|
235 |
|
|
236 DEFUN_DLD (conv2, args, , |
|
|
237 "-*- texinfo -*-\n\ |
|
|
238 @deftypefn {Loadable Function} {y =} conv2 (@var{a}, @var{b}, @var{shape})\n\ |
|
|
239 @deftypefnx {Loadable Function} {y =} conv2 (@var{v1}, @var{v2}, @var{M}, @var{shape})\n\ |
|
|
240 \n\ |
|
|
241 Returns 2D convolution of @var{a} and @var{b} where the size\n\ |
|
|
242 of @var{c} is given by\n\ |
|
|
243 \n\ |
|
|
244 @table @asis\n\ |
|
|
245 @item @var{shape}= 'full'\n\ |
|
|
246 returns full 2-D convolution\n\ |
|
|
247 @item @var{shape}= 'same'\n\ |
|
|
248 same size as a. 'central' part of convolution\n\ |
|
|
249 @item @var{shape}= 'valid'\n\ |
|
|
250 only parts which do not include zero-padded edges\n\ |
|
|
251 @end table\n\ |
|
|
252 \n\ |
|
|
253 By default @var{shape} is 'full'. When the third argument is a matrix\n\ |
|
|
254 returns the convolution of the matrix @var{M} by the vector @var{v1}\n\ |
|
|
255 in the column direction and by vector @var{v2} in the row direction\n\ |
|
|
256 @end deftypefn") |
|
|
257 { |
|
|
258 octave_value retval; |
|
|
259 octave_value tmp; |
|
|
260 int nargin = args.length (); |
|
|
261 std::string shape= "full"; //default |
|
|
262 bool separable= false; |
|
|
263 Shape ishape; |
|
|
264 |
|
|
265 if (nargin < 2) |
|
|
266 { |
|
5823
|
267 print_usage (); |
|
5819
|
268 return retval; |
|
|
269 } |
|
|
270 else if (nargin == 3) |
|
|
271 { |
|
5822
|
272 if (args(2).is_string ()) |
|
5819
|
273 shape = args(2).string_value (); |
|
|
274 else |
|
|
275 separable = true; |
|
|
276 } |
|
|
277 else if (nargin >= 4) |
|
|
278 { |
|
|
279 separable = true; |
|
|
280 shape = args(3).string_value (); |
|
|
281 } |
|
|
282 |
|
|
283 if (shape == "full") |
|
|
284 ishape = SHAPE_FULL; |
|
|
285 else if (shape == "same") |
|
|
286 ishape = SHAPE_SAME; |
|
|
287 else if (shape == "valid") |
|
|
288 ishape = SHAPE_VALID; |
|
|
289 else |
|
|
290 { |
|
5822
|
291 error ("conv2: shape type not valid"); |
|
5823
|
292 print_usage (); |
|
5819
|
293 return retval; |
|
|
294 } |
|
|
295 |
|
|
296 if (separable) |
|
|
297 { |
|
|
298 // If user requests separable, check first two params are vectors |
|
|
299 |
|
|
300 if (! (1 == args(0).rows () || 1 == args(0).columns ()) |
|
|
301 || ! (1 == args(1).rows () || 1 == args(1).columns ())) |
|
|
302 { |
|
5823
|
303 print_usage (); |
|
5819
|
304 return retval; |
|
|
305 } |
|
|
306 |
|
|
307 if (args(0).is_complex_type () |
|
|
308 || args(1).is_complex_type () |
|
|
309 || args(2).is_complex_type ()) |
|
|
310 { |
|
|
311 ComplexColumnVector v1 (args(0).complex_vector_value ()); |
|
|
312 ComplexColumnVector v2 (args(1).complex_vector_value ()); |
|
|
313 ComplexMatrix a (args(2).complex_matrix_value ()); |
|
|
314 ComplexMatrix c (conv2 (v1, v2, a, ishape)); |
|
|
315 if (! error_state) |
|
|
316 retval = c; |
|
|
317 } |
|
|
318 else |
|
|
319 { |
|
|
320 ColumnVector v1 (args(0).vector_value ()); |
|
|
321 ColumnVector v2 (args(1).vector_value ()); |
|
|
322 Matrix a (args(2).matrix_value ()); |
|
|
323 Matrix c (conv2 (v1, v2, a, ishape)); |
|
|
324 if (! error_state) |
|
|
325 retval = c; |
|
|
326 } |
|
|
327 } // if (separable) |
|
|
328 else |
|
|
329 { |
|
|
330 if (args(0).is_complex_type () |
|
|
331 || args(1).is_complex_type ()) |
|
|
332 { |
|
|
333 ComplexMatrix a (args(0).complex_matrix_value ()); |
|
|
334 ComplexMatrix b (args(1).complex_matrix_value ()); |
|
|
335 ComplexMatrix c (conv2 (a, b, ishape)); |
|
|
336 if (! error_state) |
|
|
337 retval = c; |
|
|
338 } |
|
|
339 else |
|
|
340 { |
|
|
341 Matrix a (args(0).matrix_value ()); |
|
|
342 Matrix b (args(1).matrix_value ()); |
|
|
343 Matrix c (conv2 (a, b, ishape)); |
|
|
344 if (! error_state) |
|
|
345 retval = c; |
|
|
346 } |
|
|
347 |
|
|
348 } // if (separable) |
|
|
349 |
|
|
350 return retval; |
|
|
351 } |
|
|
352 |
|
|
353 template MArray2<double> |
|
|
354 conv2 (MArray<double>&, MArray<double>&, MArray2<double>&, Shape); |
|
|
355 |
|
|
356 template MArray2<double> |
|
|
357 conv2 (MArray2<double>&, MArray2<double>&, Shape); |
|
|
358 |
|
|
359 template MArray2<Complex> |
|
|
360 conv2 (MArray<Complex>&, MArray<Complex>&, MArray2<Complex>&, Shape); |
|
|
361 |
|
|
362 template MArray2<Complex> |
|
|
363 conv2 (MArray2<Complex>&, MArray2<Complex>&, Shape); |
|
|
364 |
|
|
365 /* |
|
|
366 ;;; Local Variables: *** |
|
|
367 ;;; mode: C++ *** |
|
|
368 ;;; End: *** |
|
|
369 */ |
