1 // Copyright (C) 2008 Søren Hauberg <soren@hauberg.org>
2 // Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
3 //
4 // This program is free software; you can redistribute it and/or modify it under
5 // the terms of the GNU General Public License as published by the Free Software
6 // Foundation; either version 3 of the License, or (at your option) any later
7 // version.
8 //
9 // This program is distributed in the hope that it will be useful, but WITHOUT
10 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 // details.
13 //
14 // You should have received a copy of the GNU General Public License along with
15 // this program; if not, see <http://www.gnu.org/licenses/>.
16
17 // The 'compare' and 'selnth' functions are copied from the 'cordflt2' function
18 // developed by Teemu Ikonen. This work was released under GPLv2 or later.
19 // -- Søren Hauberg, March 21st, 2008
20
21 #include <octave/oct.h>
22
23 #define WANTS_OCTAVE_IMAGE_VALUE 1
24 #include "octave-wrappers.h"
25
26 /**
27 * Filter functions for ordered filtering.
28 */
29
30 template <class ET>
compare(const ET a,const ET b)31 inline bool compare (const ET a, const ET b)
32 {
33 if (a > b)
34 return 1;
35 else
36 return 0;
37 }
38
compare(const Complex a,const Complex b)39 template <> inline bool compare<Complex> (const Complex a, const Complex b)
40 {
41 const double anorm2 = a.real () * a.real () + a.imag () * a.imag ();
42 const double bnorm2 = b.real () * b.real () + b.imag () * b.imag ();
43
44 if (anorm2 > bnorm2)
45 return 1;
46 else
47 return 0;
48 }
49
50 // select nth largest member from the array values
51 // Partitioning algorithm, see Numerical recipes chap. 8.5
52 template <class ET, class MT, class ET_OUT>
selnth(MT & vals,octave_idx_type len,int nth)53 ET_OUT selnth (MT &vals, octave_idx_type len, int nth)
54 {
55 ET hinge;
56 int l, r, mid, i, j;
57
58 l = 0;
59 r = len - 1;
60 for (;;)
61 {
62 // if partition size is 1 or two, then sort and return
63 if (r <= l+1)
64 {
65 if (r == l+1 && compare<ET> (vals (l), vals (r)))
66 std::swap (vals (l), vals (r));
67
68 return vals (nth);
69 }
70 else
71 {
72 mid = (l+r) >> 1;
73 std::swap (vals (mid), vals (l+1));
74
75 // choose median of l, mid, r to be the hinge element
76 // and set up sentinels in the borders (order l, l+1 and r)
77 if (compare<ET> (vals (l), vals (r)))
78 std::swap (vals (l), vals (r));
79
80 if (compare<ET> (vals (l+1), vals (r)))
81 std::swap (vals (l+1), vals (r));
82
83 if (compare<ET> (vals (l), vals (l+1)))
84 std::swap (vals (l), vals (l+1));
85
86 i = l + 1;
87 j = r;
88 hinge = vals (l+1);
89 for (;;)
90 {
91 do i++; while (compare<ET> (hinge, vals (i)));
92 do j--; while (compare<ET> (vals (j), hinge));
93 if (i > j)
94 break;
95 std::swap (vals (i), vals (j));
96 }
97 vals (l+1) = vals (j);
98 vals (j) = hinge;
99 if (j >= nth)
100 r = j - 1;
101 if (j <= nth)
102 l = i;
103 }
104 }
105 }
106
107 template <class ET, class MT, class ET_OUT>
min_filt(MT & vals,octave_idx_type len,int not_used)108 ET_OUT min_filt (MT &vals, octave_idx_type len, int not_used)
109 {
110 ET_OUT min_val = vals (0);
111 for (octave_idx_type i = 1; i < len; i++)
112 min_val = compare (min_val, vals (i)) ? vals (i) : min_val;
113
114 return min_val;
115 }
116
117 template <class ET, class MT, class ET_OUT>
max_filt(MT & vals,octave_idx_type len,int not_used)118 ET_OUT max_filt (MT &vals, octave_idx_type len, int not_used)
119 {
120 ET_OUT max_val = vals (0);
121 for (octave_idx_type i = 1; i < len; i++)
122 max_val = compare (max_val, vals (i)) ? max_val : vals (i);
123
124 return max_val;
125 }
126
127 /**
128 * Filter functions for standard deviation filters
129 */
130
131 template <class ET> inline
square(const ET a)132 ET square (const ET a)
133 {
134 return a * a;
135 }
136
137 template <class ET, class MT, class ET_OUT>
std_filt(MT & vals,octave_idx_type len,int norm)138 ET_OUT std_filt (MT &vals, octave_idx_type len, int norm)
139 {
140 // Compute mean
141 ET_OUT mean = 0;
142 for (octave_idx_type i = 0; i < len; i++)
143 mean += (ET_OUT)vals (i);
144 mean /= (ET_OUT)len;
145
146 // Compute sum of square differences from the mean
147 ET_OUT var = 0;
148 for (octave_idx_type i = 0; i < len; i++)
149 var += square ((ET_OUT)vals (i) - mean);
150
151 // Normalise to produce variance
152 var /= (ET_OUT)norm;
153
154 // Compute std. deviation
155 return sqrt (var);
156 }
157
158 /**
159 * Functions for the entropy filter.
160 */
161
162 /* We only need the explicit typed versions */
163 template <class ET>
get_entropy_info(ET & add,int & nbins)164 void get_entropy_info (ET &add, int &nbins)
165 {
166 }
167
168 #define ENTROPY_INFO(TYPE, ADD, NBINS) \
169 template <> \
170 void get_entropy_info<TYPE> (TYPE &add, int &nbins) \
171 { \
172 add = ADD; \
173 if (nbins <= 0) \
174 nbins = NBINS; \
175 }
176
177 ENTROPY_INFO (bool, 0, 2)
178 ENTROPY_INFO (octave_int8, 128, 256)
179 //ENTROPY_INFO (octave_int16, 32768, 65536)
180 ENTROPY_INFO (octave_uint8, 0, 256)
181 //ENTROPY_INFO (octave_uint16, 0, 65536)
182
183 #undef ENTROPY_INFO
184
185 template <class ET, class MT, class ET_OUT>
entropy_filt(MT & vals,octave_idx_type len,int nbins)186 ET_OUT entropy_filt (MT &vals, octave_idx_type len, int nbins)
187 {
188 ET add;
189 get_entropy_info<ET> (add, nbins);
190
191 // Compute histogram from values
192 ColumnVector hist (nbins, 0);
193 for (octave_idx_type i = 0; i < len; i++)
194 hist(vals(i) + add)++;
195 for (octave_idx_type i = 0; i < nbins; i++)
196 hist(i) /= (double)len;
197
198 // Compute entropy
199 double entropy = 0;
200 for (octave_idx_type i = 0; i < nbins; i++)
201 {
202 const double p = hist (i);
203 if (p > 0)
204 entropy -= p * std::log2 (p);
205 }
206
207 return entropy;
208 }
209
210 /**
211 * The function for the range filter
212 */
213 template <class ET, class MT, class ET_OUT>
range_filt(MT & vals,octave_idx_type len,int not_used)214 ET_OUT range_filt (MT &vals, octave_idx_type len, int not_used)
215 {
216 const ET_OUT min_val = min_filt<ET, MT, ET_OUT> (vals, len, not_used);
217 const ET_OUT max_val = max_filt<ET, MT, ET_OUT> (vals, len, not_used);
218
219 return max_val - min_val;
220 }
221
222
223 // The general function for doing the filtering
224 //
225 // We differentiate between MT and MTout for cases such as std and
226 // entropy, where the output will have a different class than the input.
227
228 template <class MT, class ET, class MTout, class ETout>
do_filtering(const MT & in,const boolNDArray & se,ETout (* filter_function)(MT &,octave_idx_type,int),const MT & S,int arg4)229 octave_value do_filtering (const MT &in,
230 const boolNDArray &se,
231 ETout (*filter_function) (MT&, octave_idx_type, int),
232 const MT &S,
233 int arg4)
234 {
235 typedef typename MT::element_type Pin;
236
237 const octave_idx_type ndims = in.ndims ();
238 const octave_idx_type se_nnz = se.nnz ();
239 const dim_vector se_size = se.dims ().redim (ndims);
240 const dim_vector in_size = in.dims ();
241
242 // Create output matrix
243 dim_vector out_size (in_size);
244 for (octave_idx_type i = 0; i < ndims; i++)
245 {
246 out_size (i) = in_size (i) - se_size (i) + 1;
247 }
248 MTout out (out_size);
249
250 // We will loop for all elements of the output matrix. On each iteration
251 // we collect the values from the input matrix as marked by the
252 // structuring element (SE), and pass them to the filtering function.
253 // The value returned is then assigned assigned to the output matrix.
254 // Using dim_vector's and increment_index() allows us to support matrices
255 // with any number of dimensions.
256
257 // Create a 2D array with the subscript indices for each of the
258 // true elements on the SE. Each column has the subscripts for
259 // each true elements, and the rows are the dimensions.
260 // We also don't need the heights in a matrix. We can extract the
261 // ones that matter for us and place them in a vector to access
262 // them easily during the loop.
263 Array<octave_idx_type> se_sub (dim_vector (ndims, 1), 0);
264 Array<octave_idx_type> nnz_sub (dim_vector (ndims, se_nnz));
265 Array<Pin> heights (dim_vector (se_nnz, 1));
266 for (octave_idx_type i = 0, found = 0; found < se_nnz; i++)
267 {
268 if (se(se_sub))
269 {
270 // insert the coordinate vectors on the next column
271 nnz_sub.insert (se_sub, 0, found);
272 // store new height value
273 heights(found) = S(se_sub);
274 found++;
275 }
276 boolNDArray::increment_index (se_sub, se_size);
277 }
278
279 // Create array with subscript indexes for the elements being
280 // evaluated at any given time. We will be using linear indexes
281 // later but need the subscripts to add them.
282 Array<octave_idx_type> in_sub (dim_vector (ndims, 1));
283 Array<octave_idx_type> out_sub (dim_vector (ndims, 1), 0);
284
285 // For each iteration of the output matrix, will store the values from
286 // the input matrix that will be passed to the filtering function.
287 MT values (dim_vector (1, se_nnz));
288
289 // Get pointers to the fortran vector for good performance.
290 ETout* out_fvec = out.fortran_vec ();
291 Pin* values_fvec = values.fortran_vec ();
292 Pin* heights_fvec = heights.fortran_vec ();
293 octave_idx_type* in_sub_fvec = in_sub.fortran_vec ();
294 octave_idx_type* out_sub_fvec = out_sub.fortran_vec ();
295 octave_idx_type* nnz_sub_fvec = nnz_sub.fortran_vec ();
296
297 // We iterate for all elements of the output matrix
298 const octave_idx_type out_numel = out.numel ();
299 for (octave_idx_type out_ind = 0; out_ind < out_numel; out_ind++)
300 {
301 // On each iteration we get the subscript indexes for the output
302 // matrix (obtained with increment_index), and add to it the
303 // subscript indexes of each of the nnz elements in the SE. These
304 // are the subscript indexes for the elements in input matrix that
305 // need to be evaluated for that element in the output matrix
306 octave_idx_type nnz_sub_ind = 0;
307 for (octave_idx_type se_ind = 0; se_ind < se_nnz; se_ind++)
308 {
309 nnz_sub_ind = se_ind * ndims; // move to the next column
310 for (octave_idx_type n = 0; n < ndims; n++, nnz_sub_ind++)
311 {
312 // get subcript indexes for the input matrix
313 in_sub_fvec[n] = out_sub_fvec[n] + nnz_sub_fvec[nnz_sub_ind];
314 }
315 values_fvec[se_ind] = in(in_sub) + heights_fvec[se_ind];
316 }
317
318 // Compute filter result
319 out_fvec[out_ind] = filter_function (values, se_nnz, arg4);
320
321 // Prepare for next iteration
322 boolNDArray::increment_index (out_sub, out_size);
323 OCTAVE_QUIT;
324 }
325
326 return octave_value (out);
327 }
328
329 DEFUN_DLD(__spatial_filtering__, args, , "\
330 -*- texinfo -*-\n\
331 @deftypefn {Loadable Function} __spatial_filtering__(@var{A}, @var{domain},\
332 @var{method}, @var{S}, @var{arg})\n\
333 Implementation of two-dimensional spatial filtering. In general this function\n\
334 should NOT be used -- user interfaces are available in other functions.\n\
335 The function computes local characteristics of the image @var{A} in the domain\n\
336 @var{domain}. The following values of @var{method} are supported.\n\
337 \n\
338 @table @asis\n\
339 @item @qcode{\"ordered\"}\n\
340 Perform ordered filtering. The output in a pixel is the @math{n}th value of a\n\
341 sorted list containing the elements of the neighbourhood. The value of @math{n}\n\
342 is given in the @var{arg} argument. The corresponding user interface is available\n\
343 in @code{ordfilt2} and @code{ordfiltn}.\n\
344 \n\
345 @item @qcode{\"std\"}\n\
346 Compute the local standard deviation. The corresponding user interface is available\n\
347 in @code{stdfilt}.\n\
348 \n\
349 @item @qcode{\"entropy\"}\n\
350 Compute the local entropy. The corresponding user interface is available\n\
351 in @code{entropyfilt}.\n\
352 \n\
353 @item @qcode{\"range\"}\n\
354 Compute the local range of the data. The corresponding user interface is\n\
355 available in @code{rangefilt}.\n\
356 \n\
357 @end table\n\
358 @seealso{ordfilt2}\n\
359 @end deftypefn\n\
360 ")
361 {
362 octave_value_list retval;
363 const octave_idx_type nargin = args.length ();
364 if (nargin < 4)
365 print_usage ();
366
367 const octave_image::value A (args(0));
368 const octave_image::value S (args(3));
369
370 const boolNDArray domain = args(1).bool_array_value ();
371
372 octave_idx_type len = 0;
373 for (octave_idx_type i = 0; i < domain.numel (); i++)
374 len += domain (i);
375
376 const int ndims = domain.ndims ();
377 if (A.ndims () != ndims)
378 error ("__spatial_filtering__: A and DOMAIN must have same dimensions");
379
380 if (S.ndims () != ndims)
381 error ("__spatial_filtering__: DOMAIN and S must have same size");
382 for (octave_idx_type i = 0; i < ndims; i++)
383 if (domain.size (i) != S.dims ()(i))
384 error ("__spatial_filtering__: DOMAIN and S must have same size");
385
386 int arg4 = (nargin == 4) ? 0 : args(4).int_value ();
387
388 const std::string method = args(2).string_value ();
389
390 #define GENERAL_ACTION(MT, FUN, ET, MT_OUT, ET_OUT, FILTER_FUN) \
391 retval = do_filtering<MT, ET, MT_OUT, ET_OUT> (A.FUN (), domain, \
392 FILTER_FUN<ET, MT, ET_OUT>, \
393 S.FUN (), arg4) \
394
395 if (method == "ordered")
396 {
397 // Handle input
398 arg4 -= 1; // convert arg to zero-based index
399 if (arg4 > len - 1)
400 {
401 warning ("__spatial_filtering__: nth should be less than number of non-zero "
402 "values in domain setting nth to largest possible value");
403 arg4 = len - 1;
404 }
405 if (arg4 < 0)
406 {
407 warning ("__spatial_filtering__: nth should be non-negative, setting to 1");
408 arg4 = 0;
409 }
410
411 #define ACTION(MT, FUN, ET) \
412 GENERAL_ACTION(MT, FUN, ET, MT, ET, selnth)
413
414 if (A.is_int8_type ())
415 ACTION (int8NDArray, int8_array_value, octave_int8);
416 else if (A.is_int16_type ())
417 ACTION (int16NDArray, int16_array_value, octave_int16);
418 else if (A.is_int32_type ())
419 ACTION (int32NDArray, int32_array_value, octave_int32);
420 else if (A.is_int64_type ())
421 ACTION (int64NDArray, int64_array_value, octave_int64);
422 else if (A.is_uint8_type ())
423 ACTION (uint8NDArray, uint8_array_value, octave_uint8);
424 else if (A.is_uint16_type ())
425 ACTION (uint16NDArray, uint16_array_value, octave_uint16);
426 else if (A.is_uint32_type ())
427 ACTION (uint32NDArray, uint32_array_value, octave_uint32);
428 else if (A.is_uint64_type ())
429 ACTION (uint64NDArray, uint64_array_value, octave_uint64);
430 else if (A.islogical ())
431 ACTION (boolNDArray, bool_array_value, bool);
432 else if (A.isreal ())
433 {
434 if (A.is_single_type ())
435 ACTION (FloatNDArray, float_array_value, float);
436 else
437 ACTION (NDArray, array_value, double);
438 }
439 else if (A.iscomplex ())
440 {
441 if (A.is_single_type ())
442 ACTION (FloatComplexNDArray, float_complex_array_value, FloatComplex);
443 else
444 ACTION (ComplexNDArray, complex_array_value, Complex);
445 }
446 else
447 error ("__spatial_filtering__: A should be real, complex, or integer");
448
449 #undef ACTION
450 }
451 else if (method == "range")
452 {
453
454 #define ACTION(MT, FUN, ET) \
455 GENERAL_ACTION(MT, FUN, ET, MT, ET, range_filt)
456
457 if (A.is_int8_type ())
458 ACTION (int8NDArray, int8_array_value, octave_int8);
459 else if (A.is_int16_type ())
460 ACTION (int16NDArray, int16_array_value, octave_int16);
461 else if (A.is_int32_type ())
462 ACTION (int32NDArray, int32_array_value, octave_int32);
463 else if (A.is_int64_type ())
464 ACTION (int64NDArray, int64_array_value, octave_int64);
465 else if (A.is_uint8_type ())
466 ACTION (uint8NDArray, uint8_array_value, octave_uint8);
467 else if (A.is_uint16_type ())
468 ACTION (uint16NDArray, uint16_array_value, octave_uint16);
469 else if (A.is_uint32_type ())
470 ACTION (uint32NDArray, uint32_array_value, octave_uint32);
471 else if (A.is_uint64_type ())
472 ACTION (uint64NDArray, uint64_array_value, octave_uint64);
473 else if (A.islogical ())
474 ACTION (boolNDArray, bool_array_value, bool);
475 else if (A.isreal ())
476 {
477 if (A.is_single_type ())
478 ACTION (FloatNDArray, float_array_value, float);
479 else
480 ACTION (NDArray, array_value, double);
481 }
482 else if (A.iscomplex ())
483 {
484 if (A.is_single_type ())
485 ACTION (FloatComplexNDArray, float_complex_array_value, FloatComplex);
486 else
487 ACTION (ComplexNDArray, complex_array_value, Complex);
488 }
489 else
490 error ("__spatial_filtering__: A should be real, complex, or integer");
491
492 #undef ACTION
493 }
494 else if (method == "std")
495 {
496 // Compute normalisation factor
497 if (arg4 == 0)
498 arg4 = len - 1; // unbiased
499 else
500 arg4 = len; // max. likelihood
501
502 #define ACTION(MT, FUN, ET) \
503 GENERAL_ACTION(MT, FUN, ET, NDArray, double, std_filt)
504
505 if (A.is_int8_type ())
506 ACTION (int8NDArray, int8_array_value, octave_int8);
507 else if (A.is_int16_type ())
508 ACTION (int16NDArray, int16_array_value, octave_int16);
509 else if (A.is_int32_type ())
510 ACTION (int32NDArray, int32_array_value, octave_int32);
511 else if (A.is_int64_type ())
512 ACTION (int64NDArray, int64_array_value, octave_int64);
513 else if (A.is_uint8_type ())
514 ACTION (uint8NDArray, uint8_array_value, octave_uint8);
515 else if (A.is_uint16_type ())
516 ACTION (uint16NDArray, uint16_array_value, octave_uint16);
517 else if (A.is_uint32_type ())
518 ACTION (uint32NDArray, uint32_array_value, octave_uint32);
519 else if (A.is_uint64_type ())
520 ACTION (uint64NDArray, uint64_array_value, octave_uint64);
521 else if (A.islogical ())
522 ACTION (boolNDArray, bool_array_value, bool);
523 else if (A.is_real_matrix ())
524 {
525 if (A.is_single_type ())
526 ACTION (FloatNDArray, float_array_value, float);
527 else
528 ACTION (NDArray, array_value, double);
529 }
530 else
531 error ("__spatial_filtering__: A should be real or logical");
532
533 #undef ACTION
534 }
535 else if (method == "entropy")
536 {
537
538 #define ACTION(MT, FUN, ET) \
539 GENERAL_ACTION(MT, FUN, ET, NDArray, double, entropy_filt)
540
541 if (A.islogical ())
542 ACTION (boolNDArray, bool_array_value, bool);
543 else if (A.is_uint8_type ())
544 ACTION (uint8NDArray, uint8_array_value, octave_uint8);
545 else
546 error ("__spatial_filtering__: A should be logical or uint8");
547
548 #undef ACTION
549 }
550 else
551 error ("__spatial_filtering__: unknown method '%s'.", method.c_str ());
552
553 return retval;
554 }
555
556 /*
557 %!error <DOMAIN and S must have same size>
558 %! __spatial_filtering__ (ones (10), ones (3), "std", ones (10), 0)
559 %!error <DOMAIN and S must have same size>
560 %! __spatial_filtering__ (ones (10), ones (3), "std", ones (3, 3, 3), 0)
561 %!error <DOMAIN and S must have same size>
562 %! __spatial_filtering__ (ones (10), ones (3), "std", ones (1, 9), 0)
563
564 %!shared a, domain, s, out
565 %! a = [ 82 2 97 43 79 43 41 65 51 11
566 %! 60 65 21 56 94 77 36 38 75 39
567 %! 32 68 78 1 16 75 76 90 81 56
568 %! 43 90 82 41 36 1 87 19 18 63
569 %! 63 64 2 48 18 43 38 25 22 99
570 %! 12 46 90 79 3 92 39 79 10 22
571 %! 38 98 11 10 40 90 88 38 4 76
572 %! 54 37 9 4 33 98 36 47 53 57
573 %! 38 76 82 50 14 74 64 99 7 33
574 %! 88 96 41 62 84 89 97 23 41 3];
575 %!
576 %! domain = ones (3);
577 %! s = zeros (3);
578 %!
579 %! out = [ 2 1 1 1 16 36 36 11
580 %! 21 1 1 1 1 1 18 18
581 %! 2 1 1 1 1 1 18 18
582 %! 2 2 2 1 1 1 10 10
583 %! 2 2 2 3 3 25 4 4
584 %! 9 4 3 3 3 36 4 4
585 %! 9 4 4 4 14 36 4 4
586 %! 9 4 4 4 14 23 7 3];
587 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
588 %!
589 %! out = [ 97 97 97 94 94 90 90 90
590 %! 90 90 94 94 94 90 90 90
591 %! 90 90 82 75 87 90 90 99
592 %! 90 90 90 92 92 92 87 99
593 %! 98 98 90 92 92 92 88 99
594 %! 98 98 90 98 98 98 88 79
595 %! 98 98 82 98 98 99 99 99
596 %! 96 96 84 98 98 99 99 99];
597 %!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
598 %!
599 %! out = [ 60 43 43 43 43 43 51 51
600 %! 60 56 36 36 36 38 38 39
601 %! 63 48 18 18 36 38 25 25
602 %! 46 48 36 36 36 38 22 22
603 %! 38 46 11 40 39 39 25 22
604 %! 37 11 10 33 39 47 38 38
605 %! 38 11 11 33 40 64 38 38
606 %! 41 41 33 50 64 64 41 33];
607 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
608 %!
609 %! out = [ 31.223 33.788 35.561 31.011 26.096 20.630 20.403 24.712
610 %! 23.428 29.613 32.376 34.002 33.593 32.470 29.605 26.333
611 %! 27.834 32.890 29.903 24.207 30.083 32.497 31.898 32.600
612 %! 32.027 28.995 33.530 31.002 32.241 32.004 27.501 32.070
613 %! 34.682 36.030 33.046 33.745 32.509 27.352 28.607 34.180
614 %! 32.709 37.690 32.992 40.036 34.456 26.656 27.685 26.863
615 %! 30.971 36.227 25.775 34.873 29.917 25.269 32.292 30.410
616 %! 29.135 31.626 30.056 33.594 30.814 28.853 30.917 29.120];
617 %!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
618 %!
619 %! out = [ 95 96 96 93 78 54 54 79
620 %! 69 89 93 93 93 89 72 72
621 %! 88 89 81 74 86 89 72 81
622 %! 88 88 88 91 91 91 77 89
623 %! 96 96 88 89 89 67 84 95
624 %! 89 94 87 95 95 62 84 75
625 %! 89 94 78 94 84 63 95 95
626 %! 87 92 80 94 84 76 92 96];
627 %!assert (__spatial_filtering__ (a, domain, "range", s), out);
628 %!
629 %! domain = [ 1 1 0
630 %! 0 1 1
631 %! 0 1 0];
632 %!
633 %! out = [ 2 2 1 16 36 36 38 39
634 %! 60 1 1 16 1 36 19 18
635 %! 32 2 1 1 1 19 18 18
636 %! 2 2 18 3 1 1 19 10
637 %! 46 2 2 3 18 38 10 4
638 %! 11 9 4 3 3 36 4 4
639 %! 9 4 4 10 36 36 38 4
640 %! 37 9 4 4 33 36 7 7];
641 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
642 %!
643 %! out = [ 82 97 97 94 79 76 90 81
644 %! 90 82 56 94 94 90 90 81
645 %! 90 82 78 36 87 87 90 90
646 %! 90 90 82 43 92 87 87 99
647 %! 98 90 79 92 92 88 79 25
648 %! 98 90 90 90 98 92 79 79
649 %! 98 98 50 98 98 90 99 57
650 %! 96 82 62 84 98 99 99 53];
651 %!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
652 %!
653 %! out = [ 68 78 94 79 77 43 75 75
654 %! 78 78 41 75 77 87 81 75
655 %! 82 78 48 18 75 76 76 81
656 %! 64 90 79 41 43 39 79 22
657 %! 90 79 48 48 90 79 38 22
658 %! 46 46 79 79 92 88 47 76
659 %! 76 82 33 40 90 88 88 53
660 %! 82 50 50 74 89 98 47 47];
661 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
662 %!
663 %! out = [ 34.2389 39.2772 39.6699 31.6812 20.7364 16.5439 22.2419 17.2395
664 %! 11.9248 36.3084 21.6217 30.8350 36.4047 21.6726 30.9144 26.1017
665 %! 22.2980 33.2746 27.5808 14.5017 36.8890 29.0259 34.6020 33.2521
666 %! 32.2490 37.9579 26.9685 17.1959 32.5346 31.3847 33.5976 36.8280
667 %! 21.3354 40.1833 34.0044 33.9882 32.9894 24.1102 25.6613 9.0995
668 %! 35.4641 35.3794 39.0871 35.4753 39.9775 28.7193 26.7451 35.6553
669 %! 35.2179 45.3398 19.3210 35.2987 28.4042 24.0832 26.8421 25.0539
670 %! 23.4307 26.2812 26.3287 35.6959 25.2646 28.1016 34.9829 17.9221];
671 %!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
672 %!
673 %! out = [ 80 95 96 78 43 40 52 42
674 %! 30 81 55 78 93 54 71 63
675 %! 58 80 77 35 86 68 72 72
676 %! 88 88 64 40 91 86 68 89
677 %! 52 88 77 89 74 50 69 21
678 %! 87 81 86 87 95 56 75 75
679 %! 89 94 46 88 62 54 61 53
680 %! 59 73 58 80 65 63 92 46];
681 %!assert (__spatial_filtering__ (a, domain, "range", s), out);
682 %!
683 %! s = [ 1 -3 4
684 %! 6 -7 2
685 %! -1 3 -5];
686 %!
687 %! out = [ -1 3 4 19 38 29 31 41
688 %! 61 3 -6 9 4 33 22 21
689 %! 33 5 -2 2 -6 21 12 11
690 %! 4 -5 20 6 -2 2 16 13
691 %! 39 -1 3 -4 19 32 12 3
692 %! 13 4 3 0 4 36 6 -3
693 %! 11 2 -3 11 38 29 35 1
694 %! 34 6 1 5 34 33 9 0];
695 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
696 %!
697 %! out = [ 83 94 98 87 80 79 93 84
698 %! 93 85 53 91 95 92 83 74
699 %! 84 75 79 29 89 80 87 91
700 %! 87 93 83 45 95 84 88 101
701 %! 101 83 72 94 93 91 72 26
702 %! 91 87 91 92 101 93 76 80
703 %! 95 99 53 100 91 91 102 59
704 %! 99 75 65 87 95 101 92 50];
705 %!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
706 %!
707 %! out = [ 71 81 96 79 78 44 77 68
708 %! 80 71 44 77 78 90 83 72
709 %! 83 75 51 21 72 76 77 78
710 %! 57 91 82 42 40 42 82 20
711 %! 92 81 45 49 85 81 41 24
712 %! 43 47 76 80 90 81 50 78
713 %! 79 85 35 37 87 85 89 46
714 %! 84 52 43 76 92 100 44 48];
715 %!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
716 %!
717 %! out = [ 34.903 40.206 39.885 28.627 20.620 19.248 25.209 17.111
718 %! 14.536 35.865 23.221 32.230 34.903 23.923 28.879 22.621
719 %! 20.635 30.113 29.351 11.610 38.863 25.936 34.608 34.482
720 %! 29.811 40.998 28.279 17.897 34.666 29.978 36.150 38.213
721 %! 25.066 39.240 30.013 37.300 31.856 27.428 22.884 10.281
722 %! 31.890 34.761 39.645 37.526 39.336 27.031 25.648 39.285
723 %! 35.017 47.776 22.764 35.912 25.460 25.636 29.861 24.566
724 %! 25.213 25.000 26.391 38.451 24.631 31.305 31.118 20.611];
725 %!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
726 %!
727 %! out = [ 84 91 94 68 42 50 62 43
728 %! 32 82 59 82 91 59 61 53
729 %! 51 70 81 27 95 59 75 80
730 %! 83 98 63 39 97 82 72 88
731 %! 62 84 69 98 74 59 60 23
732 %! 78 83 88 92 97 57 70 83
733 %! 84 97 56 89 53 62 67 58
734 %! 65 69 64 82 61 68 83 50];
735 %!assert (__spatial_filtering__ (a, domain, "range", s), out);
736 */
737