Commit [f8c6b6] default Maximize Restore History

Support a new set of spatial filters

hauberg hauberg 2008-10-05

added inst/entropy.m
added inst/stdfilt.m
added inst/entropyfilt.m
added inst/rangefilt.m
changed src/Makefile
changed inst/cordflt2.m
changed inst/ordfilt2.m
changed inst/ordfiltn.m
changed INDEX
copied src/__cordfltn__.cc -> src/__spatial_filtering__.cc
inst/entropy.m Diff Switch to side-by-side view
Loading...
inst/stdfilt.m Diff Switch to side-by-side view
Loading...
inst/entropyfilt.m Diff Switch to side-by-side view
Loading...
inst/rangefilt.m Diff Switch to side-by-side view
Loading...
src/Makefile Diff Switch to side-by-side view
Loading...
inst/cordflt2.m Diff Switch to side-by-side view
Loading...
inst/ordfilt2.m Diff Switch to side-by-side view
Loading...
inst/ordfiltn.m Diff Switch to side-by-side view
Loading...
INDEX Diff Switch to side-by-side view
Loading...
src/__cordfltn__.cc to src/__spatial_filtering__.cc
--- a/src/__cordfltn__.cc
+++ b/src/__spatial_filtering__.cc
@@ -19,8 +19,10 @@
 
 #include <octave/oct.h>
 
-// Template function for comparison
-// ET is the type of the matrix element
+/**
+ * Filter functions for ordered filtering.
+ */
+
 template <class ET>
 inline bool compare (const ET a, const ET b)
 {
@@ -30,7 +32,6 @@
       return 0;
 }
 
-// Explicit template function for complex compare
 template <> inline bool compare<Complex> (const Complex a, const Complex b)
 {
     const double anorm2 = a.real () * a.real () + a.imag () * a.imag ();
@@ -44,8 +45,8 @@
 
 // select nth largest member from the array values
 // Partitioning algorithm, see Numerical recipes chap. 8.5
-template <class ET>
-ET selnth (ET *vals, int len, int nth)
+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;
@@ -57,40 +58,40 @@
       // 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];
+          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]);
+          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), 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+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]);
+          if (compare<ET> (vals (l), vals (l+1)))
+            std::swap (vals (l), vals (l+1));
             
           i = l + 1;
           j = r;
-          hinge = vals [l+1];
+          hinge = vals (l+1);
           for (;;)
             {
-              do i++; while (compare<ET> (hinge, vals [i]));
-              do j--; while (compare<ET> (vals [j], hinge));
+              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]);
+              std::swap (vals (i), vals (j));
             }
-          vals [l+1] = vals [j];
-          vals [j] = hinge;
+          vals (l+1) = vals (j);
+          vals (j) = hinge;
           if (j >= nth)
             r = j - 1;
           if (j <= nth)
@@ -99,12 +100,131 @@
     }
 }
 
-// Template function for doing the actual filtering
-// MT is the type of the matrix to be filtered (Matrix or ComplexMatrix)
-// ET is the type of the element of the matrix (double or Complex)
-template <class MT, class ET> 
-octave_value_list do_filtering (MT A, int nth, const boolNDArray dom, MT S)
-{
+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 < len; i++)
+    hist (vals (i) + add) /= (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 * 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.
+ */
+ 
+template <class MT, class ET, class MTout, class ETout> 
+octave_value_list do_filtering (const MT &A, const boolNDArray &dom,
+   ETout (*filter_function) (MT&, octave_idx_type, int), const MT &S, int arg4)
+{
+  octave_value_list retval;
+
   const int ndims = dom.ndims ();
   const octave_idx_type dom_numel = dom.numel ();
   const dim_vector dom_size = dom.dims ();
@@ -113,29 +233,17 @@
   octave_idx_type len = 0;
   for (octave_idx_type i = 0; i < dom_numel; i++)
     len += dom (i);
-  if (nth > len - 1)
-    {
-      warning("__cordfltn__: nth should be less than number of non-zero values "
-              "in domain setting nth to largest possible value");
-      nth = len - 1;
-    }
-  if (nth < 0)
-    {
-      warning("__cordfltn__: nth should be non-negative, setting to 1");
-      nth = 0; // nth is a c-index
-    }
 
   dim_vector dim_offset (dom_size);
   for (int i = 0; i < ndims; i++)
     dim_offset (i) = (dom_size (i)+1) / 2 - 1;
 
   // Allocate output
-  octave_value_list retval;
   dim_vector out_size (dom_size);
   for (int i = 0; i < ndims; i++)
     out_size (i) = A_size (i) - dom_size (i) + 1;
 
-  MT out = MT (out_size);
+  MTout out = MTout (out_size);
   const octave_idx_type out_numel = out.numel ();
 
   // Iterate over every element of 'out'.
@@ -143,10 +251,15 @@
   Array<octave_idx_type> dom_idx (idx_dim);
   Array<octave_idx_type> A_idx (idx_dim);
   Array<octave_idx_type> out_idx (idx_dim, 0);
+  
+  dim_vector values_size (2);
+  values_size (0) = 1;
+  values_size (1) = len;
+  MT values (values_size);
+  
   for (octave_idx_type i = 0; i < out_numel; i++)
     {
       // For each neighbour
-      OCTAVE_LOCAL_BUFFER (ET, values, len);
       int l = 0;
       for (int n = 0; n < ndims; n++)
         dom_idx (n) = 0;
@@ -157,13 +270,13 @@
             A_idx (n) = out_idx (n) + dom_idx (n);
        
           if (dom (dom_idx))
-            values [l++] = A (A_idx) + S (dom_idx);
+            values (l++) = A (A_idx) + S (dom_idx);
        
           dom.increment_index (dom_idx, dom_size);
         }
             
       // Compute filter result
-      out (out_idx) = selnth (values, len, nth);
+      out (out_idx) = filter_function (values, len, arg4);
 
       // Prepare for next iteration
       out.increment_index (out_idx, out_size);
@@ -176,113 +289,289 @@
   return retval;
 }
 
-// instantiate template functions
-//SH template bool compare<double>(const double, const double);
-//SH template double selnth(double *, int, int);
-//SH template Complex selnth(Complex *, int, int);
-//SH template octave_value_list do_filtering<NDArray, double>(NDArray, int, const boolNDArray, NDArray);
-// g++ is broken, explicit instantiation of specialized template function
-// confuses the compiler.
-//template int compare<Complex>(const Complex, const Complex);
-//SH template octave_value_list do_filtering<ComplexNDArray, Complex>(ComplexNDArray, int, const boolNDArray, ComplexNDArray);
-
-DEFUN_DLD(__cordfltn__, args, , "\
+/**
+ * The Octave function
+ */
+ 
+DEFUN_DLD(__spatial_filtering__, args, , "\
 -*- texinfo -*-\n\
-@deftypefn {Loadable Function} __cordfltn__(@var{A}, @var{nth}, @var{domain}, @var{S})\n\
-Implementation of two-dimensional ordered filtering. In general this function\n\
-should NOT be used. Instead use @code{ordfilt2}.\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 \"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\
+@item \"std\"\n\
+Compute the local standard deviation. The correponding user interface is available\n\
+in @code{stdfilt}.\n\
+@item \"entropy\"\n\
+Compute the local entropy. The correponding user interface is available\n\
+in @code{entropyfilt}.\n\
+@item \"range\"\n\
+Compute the local range of the data. The corresponding user interface is\n\
+available in @code{rangefilt}.\n\
+@item \"min\"\n\
+Computes the smallest value in a local neighbourheed.\n\
+@item \"max\"\n\
+Computes the largest value in a local neighbourheed.\n\
+@item \"encoded sign of difference\"\n\
+NOT IMPLEMENTED (local binary patterns style)\n\
+@end table\n\
 @seealso{cordflt2, ordfilt2}\n\
 @end deftypefn\n\
 ")
 {
   octave_value_list retval;
-  if (args.length () != 4)
+  const int nargin = args.length ();
+  if (nargin < 4)
     {
         print_usage ();
         return retval;
     }
     
-  // nth is an index to an array, thus - 1
-  const int nth = (int)args (1).scalar_value () - 1;
-  const boolNDArray dom = args (2).bool_array_value ();
+  const boolNDArray dom = args (1).bool_array_value ();
   if (error_state)
     {
-      error ("__cordfltn__: invalid input");
+      error ("__spatial_filtering__: invalid input");
       return retval;
     }
-   
+    
+  octave_idx_type len = 0;
+  for (octave_idx_type i = 0; i < dom.numel (); i++)
+    len += dom (i);
+
   const int ndims = dom.ndims ();
   const int args0_ndims = args (0).ndims ();
   if (args0_ndims != ndims || args (3).ndims () != ndims)
     {
-      error ("__cordfltn__: input must be of the same dimension");
+      error ("__spatial_filtering__: input must be of the same dimension");
       return retval;
     }
-    
-  // Take action depending on input type
-  if (args (0).is_real_matrix ())
-    {
-      const NDArray A = args (0).array_value ();
-      const NDArray S = args (3).array_value ();
-      retval = do_filtering<NDArray, double> (A, nth, dom, S);
-    } 
-  else if (args (0).is_complex_matrix ())
-    {
-      const ComplexNDArray A = args (0).complex_matrix_value ();
-      const ComplexNDArray S = args (3).complex_matrix_value ();
-      retval = do_filtering<ComplexNDArray, Complex> (A, nth, dom, S);
-    } 
-  else if (args (0).is_int8_type ())
-    {
-      const int8NDArray A = args (0).int8_array_value ();
-      const int8NDArray S = args (3).int8_array_value ();
-      retval = do_filtering<int8NDArray, octave_int8> (A, nth, dom, S);
-    } 
-  else if (args (0).is_int16_type ())
-    {
-      const int16NDArray A = args (0).int16_array_value ();
-      const int16NDArray S = args (3).int16_array_value ();
-      retval = do_filtering<int16NDArray, octave_int16> (A, nth, dom, S);
-    } 
-  else if (args (0).is_int32_type ())
-    {
-      const int32NDArray A = args (0).int32_array_value ();
-      const int32NDArray S = args (3).int32_array_value ();
-      retval = do_filtering<int32NDArray, octave_int32> (A, nth, dom, S);
-    } 
-  else if (args (0).is_int64_type ())
-    {
-      const int64NDArray A = args (0).int64_array_value ();
-      const int64NDArray S = args (3).int64_array_value ();
-      retval = do_filtering<int64NDArray, octave_int64> (A, nth, dom, S);
-    } 
-  else if (args (0).is_uint8_type ())
-    {
-      const uint8NDArray A = args (0).uint8_array_value ();
-      const uint8NDArray S = args (3).uint8_array_value ();
-      retval = do_filtering<uint8NDArray, octave_uint8> (A, nth, dom, S);
-    } 
-  else if (args (0).is_uint16_type ())
-    {
-      const uint16NDArray A = args (0).uint16_array_value ();
-      const uint16NDArray S = args (3).uint16_array_value ();
-      retval = do_filtering<uint16NDArray, octave_uint16> (A, nth, dom, S);
-    } 
-  else if (args (0).is_uint32_type ())
-    {
-      const uint32NDArray A = args (0).uint32_array_value ();
-      const uint32NDArray S = args (3).uint32_array_value ();
-      retval = do_filtering<uint32NDArray, octave_uint32> (A, nth, dom, S);
-    } 
-  else if (args (0).is_uint64_type ())
-    {
-      const uint64NDArray A = args (0).uint64_array_value ();
-      const uint64NDArray S = args (3).uint64_array_value ();
-      retval = do_filtering<uint64NDArray, octave_uint64> (A, nth, dom, S);
-    } 
+  
+  
+  int arg4 = (nargin == 4) ? 0 : args (4).int_value ();
+
+  const std::string method = args (2).string_value ();
+  if (error_state)
+    {
+      error ("__spatial_filtering__: method must be a string");
+      return retval;
+    }
+  
+  #define GENERAL_ACTION(MT, FUN, ET, MT_OUT, ET_OUT, FILTER_FUN) \
+    { \
+      const MT A = args (0).FUN (); \
+      const MT S = args (3).FUN (); \
+      if (error_state) \
+        error ("__spatial_filtering__: invalid input"); \
+      else \
+        retval = do_filtering<MT, ET, MT_OUT, ET_OUT> (A, dom, FILTER_FUN<ET, MT, ET_OUT>, S, 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;
+        }
+
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, MT, ET, selnth)
+      if (args (0).is_real_matrix ())
+        ACTION (NDArray, array_value, double)
+      else if (args (0).is_complex_matrix ())
+        ACTION (ComplexNDArray, complex_array_value, Complex)
+      else if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_int16_type ())
+        ACTION (int16NDArray, int16_array_value, octave_int16)
+      else if (args (0).is_int32_type ())
+        ACTION (int32NDArray, int32_array_value, octave_int32)
+      else if (args (0).is_int64_type ())
+        ACTION (int64NDArray, int64_array_value, octave_int64)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else if (args (0).is_uint16_type ())
+        ACTION (uint16NDArray, uint16_array_value, octave_uint16)
+      else if (args (0).is_uint32_type ())
+        ACTION (uint32NDArray, uint32_array_value, octave_uint32)
+      else if (args (0).is_uint64_type ())
+        ACTION (uint64NDArray, uint64_array_value, octave_uint64)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
+  else if (method == "min")
+    {
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, MT, ET, min_filt)
+      if (args (0).is_real_matrix ())
+        ACTION (NDArray, array_value, double)
+      else if (args (0).is_complex_matrix ())
+        ACTION (ComplexNDArray, complex_array_value, Complex)
+      else if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_int16_type ())
+        ACTION (int16NDArray, int16_array_value, octave_int16)
+      else if (args (0).is_int32_type ())
+        ACTION (int32NDArray, int32_array_value, octave_int32)
+      else if (args (0).is_int64_type ())
+        ACTION (int64NDArray, int64_array_value, octave_int64)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else if (args (0).is_uint16_type ())
+        ACTION (uint16NDArray, uint16_array_value, octave_uint16)
+      else if (args (0).is_uint32_type ())
+        ACTION (uint32NDArray, uint32_array_value, octave_uint32)
+      else if (args (0).is_uint64_type ())
+        ACTION (uint64NDArray, uint64_array_value, octave_uint64)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
+  else if (method == "max")
+    {
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, MT, ET, max_filt)
+      if (args (0).is_real_matrix ())
+        ACTION (NDArray, array_value, double)
+      else if (args (0).is_complex_matrix ())
+        ACTION (ComplexNDArray, complex_array_value, Complex)
+      else if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_int16_type ())
+        ACTION (int16NDArray, int16_array_value, octave_int16)
+      else if (args (0).is_int32_type ())
+        ACTION (int32NDArray, int32_array_value, octave_int32)
+      else if (args (0).is_int64_type ())
+        ACTION (int64NDArray, int64_array_value, octave_int64)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else if (args (0).is_uint16_type ())
+        ACTION (uint16NDArray, uint16_array_value, octave_uint16)
+      else if (args (0).is_uint32_type ())
+        ACTION (uint32NDArray, uint32_array_value, octave_uint32)
+      else if (args (0).is_uint64_type ())
+        ACTION (uint64NDArray, uint64_array_value, octave_uint64)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
+  else if (method == "range")
+    {
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, MT, ET, range_filt)
+      if (args (0).is_real_matrix ())
+        ACTION (NDArray, array_value, double)
+      else if (args (0).is_complex_matrix ())
+        ACTION (ComplexNDArray, complex_array_value, Complex)
+      else if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_int16_type ())
+        ACTION (int16NDArray, int16_array_value, octave_int16)
+      else if (args (0).is_int32_type ())
+        ACTION (int32NDArray, int32_array_value, octave_int32)
+      else if (args (0).is_int64_type ())
+        ACTION (int64NDArray, int64_array_value, octave_int64)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else if (args (0).is_uint16_type ())
+        ACTION (uint16NDArray, uint16_array_value, octave_uint16)
+      else if (args (0).is_uint32_type ())
+        ACTION (uint32NDArray, uint32_array_value, octave_uint32)
+      else if (args (0).is_uint64_type ())
+        ACTION (uint64NDArray, uint64_array_value, octave_uint64)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
+  else if (method == "std")
+    {
+      // Compute normalisation factor
+      if (arg4 == 0)
+        arg4 = len - 1; // unbiased
+      else
+        arg4 = len; // max. likelihood
+      
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, NDArray, double, std_filt)
+      if (args (0).is_real_matrix ())
+        ACTION (NDArray, array_value, double)
+      else if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_int16_type ())
+        ACTION (int16NDArray, int16_array_value, octave_int16)
+      else if (args (0).is_int32_type ())
+        ACTION (int32NDArray, int32_array_value, octave_int32)
+      else if (args (0).is_int64_type ())
+        ACTION (int64NDArray, int64_array_value, octave_int64)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else if (args (0).is_uint16_type ())
+        ACTION (uint16NDArray, uint16_array_value, octave_uint16)
+      else if (args (0).is_uint32_type ())
+        ACTION (uint32NDArray, uint32_array_value, octave_uint32)
+      else if (args (0).is_uint64_type ())
+        ACTION (uint64NDArray, uint64_array_value, octave_uint64)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
+  else if (method == "entropy")
+    {
+      // Do the real work
+      #define ACTION(MT, FUN, ET) \
+              GENERAL_ACTION(MT, FUN, ET, NDArray, double, entropy_filt)
+      if (args (0).is_bool_matrix ())
+        ACTION (boolNDArray, bool_array_value, bool)
+      else if (args (0).is_int8_type ())
+        ACTION (int8NDArray, int8_array_value, octave_int8)
+      else if (args (0).is_uint8_type ())
+        ACTION (uint8NDArray, uint8_array_value, octave_uint8)
+      else
+        error ("__spatial_filtering__: first input should be a real, complex, or integer array");
+        
+      #undef ACTION
+    }
   else
     {
-      error ("__cordfltn__: first input should be a real, complex, or integer array");
+      error ("__spatial_filtering__: unknown method '%s'.", method.c_str ());
     }
 
   return retval;