File: __spatial_filtering__.cc

package info (click to toggle)
octave-image 2.12.0-10
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 2,524 kB
  • sloc: cpp: 4,191; objc: 1,014; makefile: 30; xml: 30
file content (736 lines) | stat: -rw-r--r-- 27,325 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
// Copyright (C) 2008 Søren Hauberg <soren@hauberg.org>
// Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
//
// This program is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free Software
// Foundation; either version 3 of the License, or (at your option) any later
// version.
//
// This program is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License along with
// this program; if not, see <http://www.gnu.org/licenses/>.

// The 'compare' and 'selnth' functions are copied from the 'cordflt2' function
// developed by Teemu Ikonen. This work was released under GPLv2 or later.
//      -- Søren Hauberg, March 21st, 2008

#include <octave/oct.h>

#define WANTS_OCTAVE_IMAGE_VALUE 1
#include "octave-wrappers.h"

/**
 * Filter functions for ordered filtering.
 */

template <class ET>
inline bool compare (const ET a, const ET b)
{
    if (a > b)
      return 1;
    else
      return 0;
}

template <> inline bool compare<Complex> (const Complex a, const Complex b)
{
    const double anorm2 = a.real () * a.real () + a.imag () * a.imag ();
    const double bnorm2 = b.real () * b.real () + b.imag () * b.imag ();

    if (anorm2 > bnorm2)
      return 1;
    else
      return 0;
}

// select nth largest member from the array values
// Partitioning algorithm, see Numerical recipes chap. 8.5
template <class ET, class MT, class ET_OUT>
ET_OUT selnth (MT &vals, octave_idx_type len, int nth)
{
  ET hinge;
  int l, r, mid, i, j;

  l = 0;
  r = len - 1;
  for (;;)
    {
      // if partition size is 1 or two, then sort and return
      if (r <= l+1)
        {
          if (r == l+1 && compare<ET> (vals (l), vals (r)))
            std::swap (vals (l), vals (r));

          return vals (nth);
        }
      else
        {
          mid = (l+r) >> 1;
          std::swap (vals (mid), vals (l+1));

          // choose median of l, mid, r to be the hinge element
          // and set up sentinels in the borders (order l, l+1 and r)
          if (compare<ET> (vals (l), vals (r)))
            std::swap (vals (l), vals (r));

          if (compare<ET> (vals (l+1), vals (r)))
            std::swap (vals (l+1), vals (r));

          if (compare<ET> (vals (l), vals (l+1)))
            std::swap (vals (l), vals (l+1));

          i = l + 1;
          j = r;
          hinge = vals (l+1);
          for (;;)
            {
              do i++; while (compare<ET> (hinge, vals (i)));
              do j--; while (compare<ET> (vals (j), hinge));
              if (i > j)
                break;
              std::swap (vals (i), vals (j));
            }
          vals (l+1) = vals (j);
          vals (j) = hinge;
          if (j >= nth)
            r = j - 1;
          if (j <= nth)
            l = i;
        }
    }
}

template <class ET, class MT, class ET_OUT>
ET_OUT min_filt (MT &vals, octave_idx_type len, int not_used)
{
  ET_OUT min_val = vals (0);
  for (octave_idx_type i = 1; i < len; i++)
    min_val = compare (min_val, vals (i)) ? vals (i) : min_val;

  return min_val;
}

template <class ET, class MT, class ET_OUT>
ET_OUT max_filt (MT &vals, octave_idx_type len, int not_used)
{
  ET_OUT max_val = vals (0);
  for (octave_idx_type i = 1; i < len; i++)
    max_val = compare (max_val, vals (i)) ? max_val : vals (i);

  return max_val;
}

/**
 * Filter functions for standard deviation filters
 */

template <class ET> inline
ET square (const ET a)
{
  return a * a;
}

template <class ET, class MT, class ET_OUT>
ET_OUT std_filt (MT &vals, octave_idx_type len, int norm)
{
  // Compute mean
  ET_OUT mean = 0;
  for (octave_idx_type i = 0; i < len; i++)
    mean += (ET_OUT)vals (i);
  mean /= (ET_OUT)len;

  // Compute sum of square differences from the mean
  ET_OUT var = 0;
  for (octave_idx_type i = 0; i < len; i++)
    var += square ((ET_OUT)vals (i) - mean);

  // Normalise to produce variance
  var /= (ET_OUT)norm;

  // Compute std. deviation
  return sqrt (var);
}

/**
 * Functions for the entropy filter.
 */

/* We only need the explicit typed versions */
template <class ET>
void get_entropy_info (ET &add, int &nbins)
{
}

#define ENTROPY_INFO(TYPE, ADD, NBINS) \
  template <> \
  void get_entropy_info<TYPE> (TYPE &add, int &nbins) \
  { \
    add = ADD; \
    if (nbins <= 0) \
      nbins = NBINS; \
  }

ENTROPY_INFO (bool, 0, 2)
ENTROPY_INFO (octave_int8, 128, 256)
//ENTROPY_INFO (octave_int16, 32768, 65536)
ENTROPY_INFO (octave_uint8, 0, 256)
//ENTROPY_INFO (octave_uint16, 0, 65536)

#undef ENTROPY_INFO

template <class ET, class MT, class ET_OUT>
ET_OUT entropy_filt (MT &vals, octave_idx_type len, int nbins)
{
  ET add;
  get_entropy_info<ET> (add, nbins);

  // Compute histogram from values
  ColumnVector hist (nbins, 0);
  for (octave_idx_type i = 0; i < len; i++)
    hist(vals(i) + add)++;
  for (octave_idx_type i = 0; i < nbins; i++)
    hist(i) /= (double)len;

  // Compute entropy
  double entropy = 0;
  for (octave_idx_type i = 0; i < nbins; i++)
    {
      const double p = hist (i);
      if (p > 0)
        entropy -= p * std::log2 (p);
    }

  return entropy;
}

/**
 * The function for the range filter
 */
template <class ET, class MT, class ET_OUT>
ET_OUT range_filt (MT &vals, octave_idx_type len, int not_used)
{
  const ET_OUT min_val = min_filt<ET, MT, ET_OUT> (vals, len, not_used);
  const ET_OUT max_val = max_filt<ET, MT, ET_OUT> (vals, len, not_used);

  return max_val - min_val;
}


//      The general function for doing the filtering
//
// We differentiate between MT and MTout for cases such as std and
// entropy, where the output will have a different class than the input.

template <class MT, class ET, class MTout, class ETout>
octave_value do_filtering (const MT &in,
                           const boolNDArray &se,
                           ETout (*filter_function) (MT&, octave_idx_type, int),
                           const MT &S,
                           int arg4)
{
  typedef typename MT::element_type Pin;

  const octave_idx_type ndims     = in.ndims ();
  const octave_idx_type se_nnz    = se.nnz ();
  const dim_vector se_size        = se.dims ().redim (ndims);
  const dim_vector in_size        = in.dims ();

  // Create output matrix
  dim_vector out_size (in_size);
  for (octave_idx_type i = 0; i < ndims; i++)
    {
      out_size (i) = in_size (i) - se_size (i) + 1;
    }
  MTout out (out_size);

  // We will loop for all elements of the output matrix. On each iteration
  // we collect the values from the input matrix as marked by the
  // structuring element (SE), and pass them to the filtering function.
  // The value returned is then assigned assigned to the output matrix.
  // Using dim_vector's and increment_index() allows us to support matrices
  // with any number of dimensions.

  // Create a 2D array with the subscript indices for each of the
  // true elements on the SE. Each column has the subscripts for
  // each true elements, and the rows are the dimensions.
  // We also don't need the heights in a matrix. We can extract the
  // ones that matter for us and place them in a vector to access
  // them easily during the loop.
  Array<octave_idx_type> se_sub   (dim_vector (ndims, 1), 0);
  Array<octave_idx_type> nnz_sub  (dim_vector (ndims, se_nnz));
  Array<Pin>             heights  (dim_vector (se_nnz, 1));
  for (octave_idx_type i = 0, found = 0; found < se_nnz; i++)
    {
      if (se(se_sub))
        {
          // insert the coordinate vectors on the next column
          nnz_sub.insert (se_sub, 0, found);
          // store new height value
          heights(found) = S(se_sub);
          found++;
        }
      boolNDArray::increment_index (se_sub, se_size);
    }

  // Create array with subscript indexes for the elements being
  // evaluated at any given time. We will be using linear indexes
  // later but need the subscripts to add them.
  Array<octave_idx_type> in_sub  (dim_vector (ndims, 1));
  Array<octave_idx_type> out_sub (dim_vector (ndims, 1), 0);

  // For each iteration of the output matrix, will store the values from
  // the input matrix that will be passed to the filtering function.
  MT values (dim_vector (1, se_nnz));

  // Get pointers to the fortran vector for good performance.
  ETout* out_fvec               = out.fortran_vec ();
  Pin* values_fvec              = values.fortran_vec ();
  Pin* heights_fvec             = heights.fortran_vec ();
  octave_idx_type* in_sub_fvec  = in_sub.fortran_vec ();
  octave_idx_type* out_sub_fvec = out_sub.fortran_vec ();
  octave_idx_type* nnz_sub_fvec = nnz_sub.fortran_vec ();

  // We iterate for all elements of the output matrix
  const octave_idx_type out_numel = out.numel ();
  for (octave_idx_type out_ind = 0; out_ind < out_numel; out_ind++)
    {
      // On each iteration we get the subscript indexes for the output
      // matrix (obtained with increment_index), and add to it the
      // subscript indexes of each of the nnz elements in the SE. These
      // are the subscript indexes for the elements in input matrix that
      // need to be evaluated for that element in the output matrix
      octave_idx_type nnz_sub_ind = 0;
      for (octave_idx_type se_ind = 0; se_ind < se_nnz; se_ind++)
        {
          nnz_sub_ind = se_ind * ndims; // move to the next column
          for (octave_idx_type n = 0; n < ndims; n++, nnz_sub_ind++)
            {
              // get subcript indexes for the input matrix
              in_sub_fvec[n] = out_sub_fvec[n] + nnz_sub_fvec[nnz_sub_ind];
            }
          values_fvec[se_ind] = in(in_sub) + heights_fvec[se_ind];
        }

      // Compute filter result
      out_fvec[out_ind] = filter_function (values, se_nnz, arg4);

      // Prepare for next iteration
      boolNDArray::increment_index (out_sub, out_size);
      OCTAVE_QUIT;
    }

  return octave_value (out);
}

DEFUN_DLD(__spatial_filtering__, args, , "\
-*- texinfo -*-\n\
@deftypefn {Loadable Function} __spatial_filtering__(@var{A}, @var{domain},\
@var{method}, @var{S}, @var{arg})\n\
Implementation of two-dimensional spatial filtering. In general this function\n\
should NOT be used -- user interfaces are available in other functions.\n\
The function computes local characteristics of the image @var{A} in the domain\n\
@var{domain}. The following values of @var{method} are supported.\n\
\n\
@table @asis\n\
@item @qcode{\"ordered\"}\n\
Perform ordered filtering. The output in a pixel is the @math{n}th value of a\n\
sorted list containing the elements of the neighbourhood. The value of @math{n}\n\
is given in the @var{arg} argument. The corresponding user interface is available\n\
in @code{ordfilt2} and @code{ordfiltn}.\n\
\n\
@item @qcode{\"std\"}\n\
Compute the local standard deviation. The corresponding user interface is available\n\
in @code{stdfilt}.\n\
\n\
@item @qcode{\"entropy\"}\n\
Compute the local entropy. The corresponding user interface is available\n\
in @code{entropyfilt}.\n\
\n\
@item @qcode{\"range\"}\n\
Compute the local range of the data. The corresponding user interface is\n\
available in @code{rangefilt}.\n\
\n\
@end table\n\
@seealso{ordfilt2}\n\
@end deftypefn\n\
")
{
  octave_value_list retval;
  const octave_idx_type nargin = args.length ();
  if (nargin < 4)
    print_usage ();

  const octave_image::value A (args(0));
  const octave_image::value S (args(3));

  const boolNDArray domain = args(1).bool_array_value ();

  octave_idx_type len = 0;
  for (octave_idx_type i = 0; i < domain.numel (); i++)
    len += domain (i);

  const int ndims = domain.ndims ();
  if (A.ndims () != ndims)
    error ("__spatial_filtering__: A and DOMAIN must have same dimensions");

  if (S.ndims () != ndims)
    error ("__spatial_filtering__: DOMAIN and S must have same size");
  for (octave_idx_type i = 0; i < ndims; i++)
    if (domain.size (i) != S.dims ()(i))
      error ("__spatial_filtering__: DOMAIN and S must have same size");

  int arg4 = (nargin == 4) ? 0 : args(4).int_value ();

  const std::string method = args(2).string_value ();

  #define GENERAL_ACTION(MT, FUN, ET, MT_OUT, ET_OUT, FILTER_FUN) \
    retval = do_filtering<MT, ET, MT_OUT, ET_OUT> (A.FUN (), domain, \
                                                   FILTER_FUN<ET, MT, ET_OUT>, \
                                                   S.FUN (), arg4) \

  if (method == "ordered")
    {
      // Handle input
      arg4 -= 1; // convert arg to zero-based index
      if (arg4 > len - 1)
        {
          warning ("__spatial_filtering__: nth should be less than number of non-zero "
                   "values in domain setting nth to largest possible value");
          arg4 = len - 1;
        }
      if (arg4 < 0)
        {
          warning ("__spatial_filtering__: nth should be non-negative, setting to 1");
          arg4 = 0;
        }

      #define ACTION(MT, FUN, ET) \
        GENERAL_ACTION(MT, FUN, ET, MT, ET, selnth)

      if (A.is_int8_type ())
        ACTION (int8NDArray, int8_array_value, octave_int8);
      else if (A.is_int16_type ())
        ACTION (int16NDArray, int16_array_value, octave_int16);
      else if (A.is_int32_type ())
        ACTION (int32NDArray, int32_array_value, octave_int32);
      else if (A.is_int64_type ())
        ACTION (int64NDArray, int64_array_value, octave_int64);
      else if (A.is_uint8_type ())
        ACTION (uint8NDArray, uint8_array_value, octave_uint8);
      else if (A.is_uint16_type ())
        ACTION (uint16NDArray, uint16_array_value, octave_uint16);
      else if (A.is_uint32_type ())
        ACTION (uint32NDArray, uint32_array_value, octave_uint32);
      else if (A.is_uint64_type ())
        ACTION (uint64NDArray, uint64_array_value, octave_uint64);
      else if (A.islogical ())
        ACTION (boolNDArray, bool_array_value, bool);
      else if (A.isreal ())
        {
          if (A.is_single_type ())
            ACTION (FloatNDArray, float_array_value, float);
          else
            ACTION (NDArray, array_value, double);
        }
      else if (A.iscomplex ())
        {
          if (A.is_single_type ())
            ACTION (FloatComplexNDArray, float_complex_array_value, FloatComplex);
          else
            ACTION (ComplexNDArray, complex_array_value, Complex);
        }
      else
        error ("__spatial_filtering__: A should be real, complex, or integer");

      #undef ACTION
    }
  else if (method == "range")
    {

      #define ACTION(MT, FUN, ET) \
        GENERAL_ACTION(MT, FUN, ET, MT, ET, range_filt)

      if (A.is_int8_type ())
        ACTION (int8NDArray, int8_array_value, octave_int8);
      else if (A.is_int16_type ())
        ACTION (int16NDArray, int16_array_value, octave_int16);
      else if (A.is_int32_type ())
        ACTION (int32NDArray, int32_array_value, octave_int32);
      else if (A.is_int64_type ())
        ACTION (int64NDArray, int64_array_value, octave_int64);
      else if (A.is_uint8_type ())
        ACTION (uint8NDArray, uint8_array_value, octave_uint8);
      else if (A.is_uint16_type ())
        ACTION (uint16NDArray, uint16_array_value, octave_uint16);
      else if (A.is_uint32_type ())
        ACTION (uint32NDArray, uint32_array_value, octave_uint32);
      else if (A.is_uint64_type ())
        ACTION (uint64NDArray, uint64_array_value, octave_uint64);
      else if (A.islogical ())
        ACTION (boolNDArray, bool_array_value, bool);
      else if (A.isreal ())
        {
          if (A.is_single_type ())
            ACTION (FloatNDArray, float_array_value, float);
          else
            ACTION (NDArray, array_value, double);
        }
      else if (A.iscomplex ())
        {
          if (A.is_single_type ())
            ACTION (FloatComplexNDArray, float_complex_array_value, FloatComplex);
          else
            ACTION (ComplexNDArray, complex_array_value, Complex);
        }
      else
        error ("__spatial_filtering__: A should be real, complex, or integer");

      #undef ACTION
    }
  else if (method == "std")
    {
      // Compute normalisation factor
      if (arg4 == 0)
        arg4 = len - 1; // unbiased
      else
        arg4 = len; // max. likelihood

      #define ACTION(MT, FUN, ET) \
        GENERAL_ACTION(MT, FUN, ET, NDArray, double, std_filt)

      if (A.is_int8_type ())
        ACTION (int8NDArray, int8_array_value, octave_int8);
      else if (A.is_int16_type ())
        ACTION (int16NDArray, int16_array_value, octave_int16);
      else if (A.is_int32_type ())
        ACTION (int32NDArray, int32_array_value, octave_int32);
      else if (A.is_int64_type ())
        ACTION (int64NDArray, int64_array_value, octave_int64);
      else if (A.is_uint8_type ())
        ACTION (uint8NDArray, uint8_array_value, octave_uint8);
      else if (A.is_uint16_type ())
        ACTION (uint16NDArray, uint16_array_value, octave_uint16);
      else if (A.is_uint32_type ())
        ACTION (uint32NDArray, uint32_array_value, octave_uint32);
      else if (A.is_uint64_type ())
        ACTION (uint64NDArray, uint64_array_value, octave_uint64);
      else if (A.islogical ())
        ACTION (boolNDArray, bool_array_value, bool);
      else if (A.is_real_matrix ())
        {
          if (A.is_single_type ())
            ACTION (FloatNDArray, float_array_value, float);
          else
            ACTION (NDArray, array_value, double);
        }
      else
        error ("__spatial_filtering__: A should be real or logical");

      #undef ACTION
    }
  else if (method == "entropy")
    {

      #define ACTION(MT, FUN, ET) \
        GENERAL_ACTION(MT, FUN, ET, NDArray, double, entropy_filt)

      if (A.islogical ())
        ACTION (boolNDArray, bool_array_value, bool);
      else if (A.is_uint8_type ())
        ACTION (uint8NDArray, uint8_array_value, octave_uint8);
      else
        error ("__spatial_filtering__: A should be logical or uint8");

      #undef ACTION
    }
  else
    error ("__spatial_filtering__: unknown method '%s'.", method.c_str ());

  return retval;
}

/*
%!error <DOMAIN and S must have same size>
%!  __spatial_filtering__ (ones (10), ones (3), "std", ones (10), 0)
%!error <DOMAIN and S must have same size>
%!  __spatial_filtering__ (ones (10), ones (3), "std", ones (3, 3, 3), 0)
%!error <DOMAIN and S must have same size>
%!  __spatial_filtering__ (ones (10), ones (3), "std", ones (1, 9), 0)

%!shared a, domain, s, out
%! a = [ 82    2   97   43   79   43   41   65   51   11
%!       60   65   21   56   94   77   36   38   75   39
%!       32   68   78    1   16   75   76   90   81   56
%!       43   90   82   41   36    1   87   19   18   63
%!       63   64    2   48   18   43   38   25   22   99
%!       12   46   90   79    3   92   39   79   10   22
%!       38   98   11   10   40   90   88   38    4   76
%!       54   37    9    4   33   98   36   47   53   57
%!       38   76   82   50   14   74   64   99    7   33
%!       88   96   41   62   84   89   97   23   41    3];
%!
%! domain = ones  (3);
%! s      = zeros (3);
%!
%! out = [  2    1    1    1   16   36   36   11
%!         21    1    1    1    1    1   18   18
%!          2    1    1    1    1    1   18   18
%!          2    2    2    1    1    1   10   10
%!          2    2    2    3    3   25    4    4
%!          9    4    3    3    3   36    4    4
%!          9    4    4    4   14   36    4    4
%!          9    4    4    4   14   23    7    3];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
%!
%! out = [ 97   97   97   94   94   90   90   90
%!         90   90   94   94   94   90   90   90
%!         90   90   82   75   87   90   90   99
%!         90   90   90   92   92   92   87   99
%!         98   98   90   92   92   92   88   99
%!         98   98   90   98   98   98   88   79
%!         98   98   82   98   98   99   99   99
%!         96   96   84   98   98   99   99   99];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
%!
%! out = [ 60   43   43   43   43   43   51   51
%!         60   56   36   36   36   38   38   39
%!         63   48   18   18   36   38   25   25
%!         46   48   36   36   36   38   22   22
%!         38   46   11   40   39   39   25   22
%!         37   11   10   33   39   47   38   38
%!         38   11   11   33   40   64   38   38
%!         41   41   33   50   64   64   41   33];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
%!
%! out = [ 31.223   33.788   35.561   31.011   26.096   20.630   20.403   24.712
%!         23.428   29.613   32.376   34.002   33.593   32.470   29.605   26.333
%!         27.834   32.890   29.903   24.207   30.083   32.497   31.898   32.600
%!         32.027   28.995   33.530   31.002   32.241   32.004   27.501   32.070
%!         34.682   36.030   33.046   33.745   32.509   27.352   28.607   34.180
%!         32.709   37.690   32.992   40.036   34.456   26.656   27.685   26.863
%!         30.971   36.227   25.775   34.873   29.917   25.269   32.292   30.410
%!         29.135   31.626   30.056   33.594   30.814   28.853   30.917   29.120];
%!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
%!
%! out = [ 95   96   96   93   78   54   54   79
%!         69   89   93   93   93   89   72   72
%!         88   89   81   74   86   89   72   81
%!         88   88   88   91   91   91   77   89
%!         96   96   88   89   89   67   84   95
%!         89   94   87   95   95   62   84   75
%!         89   94   78   94   84   63   95   95
%!         87   92   80   94   84   76   92   96];
%!assert (__spatial_filtering__ (a, domain, "range", s), out);
%!
%! domain = [ 1 1 0
%!            0 1 1
%!            0 1 0];
%!
%! out = [  2    2    1   16   36   36   38   39
%!         60    1    1   16    1   36   19   18
%!         32    2    1    1    1   19   18   18
%!          2    2   18    3    1    1   19   10
%!         46    2    2    3   18   38   10    4
%!         11    9    4    3    3   36    4    4
%!          9    4    4   10   36   36   38    4
%!         37    9    4    4   33   36    7    7];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
%!
%! out = [ 82   97   97   94   79   76   90   81
%!         90   82   56   94   94   90   90   81
%!         90   82   78   36   87   87   90   90
%!         90   90   82   43   92   87   87   99
%!         98   90   79   92   92   88   79   25
%!         98   90   90   90   98   92   79   79
%!         98   98   50   98   98   90   99   57
%!         96   82   62   84   98   99   99   53];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
%!
%! out = [ 68   78   94   79   77   43   75   75
%!         78   78   41   75   77   87   81   75
%!         82   78   48   18   75   76   76   81
%!         64   90   79   41   43   39   79   22
%!         90   79   48   48   90   79   38   22
%!         46   46   79   79   92   88   47   76
%!         76   82   33   40   90   88   88   53
%!         82   50   50   74   89   98   47   47];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
%!
%! out = [ 34.2389   39.2772   39.6699   31.6812   20.7364   16.5439   22.2419   17.2395
%!         11.9248   36.3084   21.6217   30.8350   36.4047   21.6726   30.9144   26.1017
%!         22.2980   33.2746   27.5808   14.5017   36.8890   29.0259   34.6020   33.2521
%!         32.2490   37.9579   26.9685   17.1959   32.5346   31.3847   33.5976   36.8280
%!         21.3354   40.1833   34.0044   33.9882   32.9894   24.1102   25.6613    9.0995
%!         35.4641   35.3794   39.0871   35.4753   39.9775   28.7193   26.7451   35.6553
%!         35.2179   45.3398   19.3210   35.2987   28.4042   24.0832   26.8421   25.0539
%!         23.4307   26.2812   26.3287   35.6959   25.2646   28.1016   34.9829   17.9221];
%!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
%!
%! out = [ 80   95   96   78   43   40   52   42
%!         30   81   55   78   93   54   71   63
%!         58   80   77   35   86   68   72   72
%!         88   88   64   40   91   86   68   89
%!         52   88   77   89   74   50   69   21
%!         87   81   86   87   95   56   75   75
%!         89   94   46   88   62   54   61   53
%!         59   73   58   80   65   63   92   46];
%!assert (__spatial_filtering__ (a, domain, "range", s), out);
%!
%! s = [  1  -3   4
%!        6  -7   2
%!       -1   3  -5];
%!
%! out = [ -1    3    4   19   38   29   31   41
%!         61    3   -6    9    4   33   22   21
%!         33    5   -2    2   -6   21   12   11
%!          4   -5   20    6   -2    2   16   13
%!         39   -1    3   -4   19   32   12    3
%!         13    4    3    0    4   36    6   -3
%!         11    2   -3   11   38   29   35    1
%!         34    6    1    5   34   33    9    0];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 1), out);
%!
%! out = [  83    94    98    87    80    79    93    84
%!          93    85    53    91    95    92    83    74
%!          84    75    79    29    89    80    87    91
%!          87    93    83    45    95    84    88   101
%!         101    83    72    94    93    91    72    26
%!          91    87    91    92   101    93    76    80
%!          95    99    53   100    91    91   102    59
%!          99    75    65    87    95   101    92    50];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, nnz (domain)), out);
%!
%! out = [  71    81    96    79    78    44    77    68
%!          80    71    44    77    78    90    83    72
%!          83    75    51    21    72    76    77    78
%!          57    91    82    42    40    42    82    20
%!          92    81    45    49    85    81    41    24
%!          43    47    76    80    90    81    50    78
%!          79    85    35    37    87    85    89    46
%!          84    52    43    76    92   100    44    48];
%!assert (__spatial_filtering__ (a, domain, "ordered", s, 4), out);
%!
%! out = [ 34.903   40.206   39.885   28.627   20.620   19.248   25.209   17.111
%!         14.536   35.865   23.221   32.230   34.903   23.923   28.879   22.621
%!         20.635   30.113   29.351   11.610   38.863   25.936   34.608   34.482
%!         29.811   40.998   28.279   17.897   34.666   29.978   36.150   38.213
%!         25.066   39.240   30.013   37.300   31.856   27.428   22.884   10.281
%!         31.890   34.761   39.645   37.526   39.336   27.031   25.648   39.285
%!         35.017   47.776   22.764   35.912   25.460   25.636   29.861   24.566
%!         25.213   25.000   26.391   38.451   24.631   31.305   31.118   20.611];
%!assert (__spatial_filtering__ (a, domain, "std", s), out, 0.001);
%!
%! out = [ 84   91   94   68   42   50   62   43
%!         32   82   59   82   91   59   61   53
%!         51   70   81   27   95   59   75   80
%!         83   98   63   39   97   82   72   88
%!         62   84   69   98   74   59   60   23
%!         78   83   88   92   97   57   70   83
%!         84   97   56   89   53   62   67   58
%!         65   69   64   82   61   68   83   50];
%!assert (__spatial_filtering__ (a, domain, "range", s), out);
*/