Mercurial > hg > octave-image

changeset 350:f8c6b6fa1217

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Support a new set of spatial filters
author hauberg
date Sun, 05 Oct 2008 19:18:25 +0000
parents 327358a87fdc
children 9b16a35cf8d1
files INDEX inst/cordflt2.m inst/entropy.m inst/entropyfilt.m inst/ordfilt2.m inst/ordfiltn.m inst/rangefilt.m inst/stdfilt.m src/Makefile src/__cordfltn__.cc src/__spatial_filtering__.cc
diffstat 11 files changed, 901 insertions(+), 302 deletions(-) [+]
line wrap: on
line diff
--- a/INDEX
+++ b/INDEX
@@ -28,6 +28,7 @@
  imhist 
  mean2 
  std2 
+ entropy
  qtdecomp
  qtgetblk
  qtsetblk
@@ -118,3 +119,6 @@
  nlfilter
  im2col
  col2im
+ rangefilt
+ stdfilt
+ entropyfilt
--- a/inst/cordflt2.m
+++ b/inst/cordflt2.m
@@ -20,8 +20,9 @@
 ## @seealso{ordfilt2}
 ## @end deftypefn
 
-function varargout = cordflt2(varargin)
+function retval = cordflt2(A, nth, domain, S)
   warning(["cordflt2: this function is deprecated and will be removed in upcoming "
            "releases. Use 'ordfilt2' instead."]);
-  [varargout{1:nargout}] = __cordfltn__(varargin{:});
+
+  retval = __spatial_filtering__ (A, domain, "ordered", S, nth);
 endfunction
new file mode 100644
--- /dev/null
+++ b/inst/entropy.m
@@ -0,0 +1,67 @@
+## Copyright (C) 2008 Søren Hauberg
+## 
+## 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, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{E} =} entropy (@var{im})
+## @deftypefnx{Function File} {@var{E} =} entropy (@var{im}, @var{nbins})
+## Computes the entropy of an image.
+##
+## The entropy of the elements of the image @var{im} is computed as
+##
+## @example
+## @var{E} = -sum (@var{P} .* log2 (@var{P})
+## @end example
+##
+## where @var{P} is the distribution of the elements of @var{im}. The distribution
+## is approximated using a histogram with @var{nbins} cells. If @var{im} is
+## @code{logical} then two cells are used by default. For other classes 256 cells
+## are used by default.
+##
+## When the entropy is computed, zero-valued cells of the histogram are ignored.
+##
+## @seealso{entropyfilt}
+## @end deftypefn
+
+function retval = entropy (I, nbins = 0)
+  ## Check input
+  if (nargin == 0)
+    error ("entropy: not enough input arguments");
+  endif
+  
+  if (!ismatrix (I))
+    error ("entropy: first input must be a matrix");
+  endif
+  
+  if (!isscalar (nbins))
+    error ("entropy: second input argument must be a scalar");
+  endif
+  
+  ## Get number of histogram bins
+  if (nbins <= 0)
+    if (islogical (I))
+      nbins = 2;
+    else
+      nbins = 256;
+    endif
+  endif
+  
+  ## Compute histogram
+  P = hist (I (:), nbins, true);
+  
+  ## Compute entropy (ignoring zero-entries of the histogram)
+  P += (P == 0);
+  retval = -sum (P .* log2 (P));
+endfunction
new file mode 100644
--- /dev/null
+++ b/inst/entropyfilt.m
@@ -0,0 +1,100 @@
+## Copyright (C) 2008 Søren Hauberg
+## 
+## 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, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{E} =} entropyfilt (@var{im})
+## @deftypefnx{Function File} {@var{E} =} entropyfilt (@var{im}, @var{domain})
+## @deftypefnx{Function File} {@var{E} =} entropyfilt (@var{im}, @var{domain}, @var{padding}, ...)
+## Computes the local entropy in a neighbourhood around each pixel in an image.
+##
+## The entropy of the elements of the neighbourhood is computed as
+##
+## @example
+## @var{E} = -sum (@var{P} .* log2 (@var{P})
+## @end example
+##
+## where @var{P} is the distribution of the elements of @var{im}. The distribution
+## is approximated using a histogram with @var{nbins} cells. If @var{im} is
+## @code{logical} then two cells are used. For other classes 256 cells
+## are used.
+##
+## When the entropy is computed, zero-valued cells of the histogram are ignored.
+##
+## The neighbourhood is defined by the @var{domain} binary mask. Elements of the
+## mask with a non-zero value are considered part of the neighbourhood. By default
+## a 9 by 9 matrix containing only non-zero values is used.
+##
+## At the border of the image, extrapolation is used. By default symmetric
+## extrapolation is used, but any method supported by the @code{padarray} function
+## can be used. Since extrapolation is used, one can expect a lower entropy near
+## the image border.
+##
+## @seealso{entropy, paddarray, stdfilt}
+## @end deftypefn
+
+function retval = entropyfilt (I, domain = true (9), padding = "symmetric", varargin)
+  ## Check input
+  if (nargin == 0)
+    error ("entropyfilt: not enough input arguments");
+  endif
+  
+  if (!ismatrix (I))
+    error ("entropyfilt: first input must be a matrix");
+  endif
+  
+  if (!ismatrix (domain))
+    error ("entropyfilt: second input argument must be a logical matrix");
+  endif
+  domain = (domain > 0);
+  
+  ## Get number of histogram bins
+  if (islogical (I))
+    nbins = 2;
+  else
+    nbins = 256;
+  endif
+  
+  ## Convert to 8 or 16 bit integers if needed
+  switch (class (I))
+    case {"double", "single", "int16", "int32", "int64", "uint16", "uint32", "uint64"}
+      min_val = double (min (I (:)));
+      max_val = double (max (I (:)));
+      if (min_val == max_val)
+        retval = zeros (size (I));
+        return;
+      endif
+      I = (double (I) - min_val)./(max_val - min_val);
+      I = uint8 (255 * I);
+    case {"logical", "int8", "uint8"}
+      ## Do nothing
+    otherwise
+      error ("entropyfilt: cannot handle images of class '%s'", class (I));
+  endswitch
+  size (I)
+  ## Pad image
+  pad = floor (size (domain)/2);
+  I = padarray (I, pad, padding, varargin {:});
+  even = (round (size (domain)/2) == size (domain)/2);
+  idx = cell (1, ndims (I));
+  for k = 1:ndims (I)
+    idx {k} = (even (k)+1):size (I, k);
+  endfor
+  I = I (idx {:});
+  size (I)
+  ## Perform filtering
+  retval = __spatial_filtering__ (I, domain, "entropy", I, nbins);
+
+endfunction
--- a/inst/ordfilt2.m
+++ b/inst/ordfilt2.m
@@ -63,6 +63,6 @@
 endif;
 
 A = impad(A, xpad, ypad, padding);
-retval = __cordfltn__(A, nth, domain, S);
+retval = __spatial_filtering__ (A, domain, "ordered", S, nth);
 
 endfunction
--- a/inst/ordfiltn.m
+++ b/inst/ordfiltn.m
@@ -95,5 +95,5 @@
   A = A(idx{:});
   
   ## Perform the filtering
-  retval = __cordfltn__(A, nth, domain, S);
+  retval = __spatial_filtering__ (A, domain, "ordered", S, nth);
 endfunction
new file mode 100644
--- /dev/null
+++ b/inst/rangefilt.m
@@ -0,0 +1,71 @@
+## Copyright (C) 2008 Søren Hauberg
+## 
+## 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, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{R} =} rangefilt (@var{im})
+## @deftypefnx{Function File} {@var{R} =} rangefilt (@var{im}, @var{domain})
+## @deftypefnx{Function File} {@var{R} =} rangefilt (@var{im}, @var{domain}, @var{padding}, ...)
+## Computes the local intensity range in a neighbourhood around each pixel in
+## an image.
+##
+## The intensity range of the pixels of a neighbourhood is computed as
+##
+## @example
+## @var{R} = max (@var{x}) - min (@var{x})
+## @end example
+##
+## where @var{x} is the value of the pixels in the neighbourhood,
+##
+## The neighbourhood is defined by the @var{domain} binary mask. Elements of the
+## mask with a non-zero value are considered part of the neighbourhood. By default
+## a 3 by 3 matrix containing only non-zero values is used.
+##
+## At the border of the image, extrapolation is used. By default symmetric
+## extrapolation is used, but any method supported by the @code{padarray} function
+## can be used.
+##
+## @seealso{paddarray, entropyfilt, stdfilt}
+## @end deftypefn
+
+function retval = rangefilt (I, domain = true (3), padding = "symmetric", varargin)
+  ## Check input
+  if (nargin == 0)
+    error ("rangefilt: not enough input arguments");
+  endif
+  
+  if (!ismatrix (I))
+    error ("rangefilt: first input must be a matrix");
+  endif
+  
+  if (!ismatrix (domain))
+    error ("rangefilt: second input argument must be a logical matrix");
+  endif
+  domain = (domain > 0);
+  
+  ## Pad image
+  pad = floor (size (domain)/2);
+  I = padarray (I, pad, padding, varargin {:});
+  even = (round (size (domain)/2) == size (domain)/2);
+  idx = cell (1, ndims (I));
+  for k = 1:ndims (I)
+    idx {k} = (even (k)+1):size (I, k);
+  endfor
+  I = I (idx {:});
+
+  ## Perform filtering
+  retval = __spatial_filtering__ (I, domain, "range", I, 0);
+
+endfunction
new file mode 100644
--- /dev/null
+++ b/inst/stdfilt.m
@@ -0,0 +1,74 @@
+## Copyright (C) 2008 Søren Hauberg
+## 
+## 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, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{S} =} stdfilt (@var{im})
+## @deftypefnx{Function File} {@var{S} =} stdfilt (@var{im}, @var{domain})
+## @deftypefnx{Function File} {@var{S} =} stdfilt (@var{im}, @var{domain}, @var{padding}, ...)
+## Computes the local standard deviation in a neighbourhood around each pixel in
+## an image.
+##
+## The standard deviation of the pixels of a neighbourhood is computed as
+##
+## @example
+## @var{S} = sqrt ((sum (@var{x} - @var{mu}).^2)/(@var{N}-1))
+## @end example
+##
+## where @var{mu} is the mean value of the pixels in the neighbourhood,
+## @var{N} is the number of pixels in the neighbourhood. So, an unbiased estimator
+## is used.
+##
+## The neighbourhood is defined by the @var{domain} binary mask. Elements of the
+## mask with a non-zero value are considered part of the neighbourhood. By default
+## a 3 by 3 matrix containing only non-zero values is used.
+##
+## At the border of the image, extrapolation is used. By default symmetric
+## extrapolation is used, but any method supported by the @code{padarray} function
+## can be used. Since extrapolation is used, one can expect a lower deviation near
+## the image border.
+##
+## @seealso{std2, paddarray, entropyfilt}
+## @end deftypefn
+
+function retval = stdfilt (I, domain = true (3), padding = "symmetric", varargin)
+  ## Check input
+  if (nargin == 0)
+    error ("stdfilt: not enough input arguments");
+  endif
+  
+  if (!ismatrix (I))
+    error ("stdfilt: first input must be a matrix");
+  endif
+  
+  if (!ismatrix (domain))
+    error ("stdfilt: second input argument must be a logical matrix");
+  endif
+  domain = (domain > 0);
+  
+  ## Pad image
+  pad = floor (size (domain)/2);
+  I = padarray (I, pad, padding, varargin {:});
+  even = (round (size (domain)/2) == size (domain)/2);
+  idx = cell (1, ndims (I));
+  for k = 1:ndims (I)
+    idx {k} = (even (k)+1):size (I, k);
+  endfor
+  I = I (idx {:});
+
+  ## Perform filtering
+  retval = __spatial_filtering__ (I, domain, "std", I, 0);
+
+endfunction
--- a/src/Makefile
+++ b/src/Makefile
@@ -8,14 +8,10 @@
 	PNG=pngread.oct pngwrite.oct
 endif
 
-ifdef HAVE_MAGICKXX
-	IMAGEMAGICK=__magick_read__.oct
-endif
-
-all: __cordfltn__.oct __bilateral__.oct __custom_gaussian_smoothing__.oct \
+all: __spatial_filtering__.oct __bilateral__.oct __custom_gaussian_smoothing__.oct \
 	bwlabel.oct bwfill.oct rotate_scale.oct hough_line.oct \
 	graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \
-	$(JPEG) $(PNG) $(IMAGEMAGICK)
+	$(JPEG) $(PNG)
 
 jpgread.oct: jpgread.cc
 	$(MKOCTFILE) $< -ljpeg
@@ -29,7 +25,4 @@
 pngwrite.oct: pngwrite.cc
 	$(MKOCTFILE) $< -lpng
 
-__magick_read__.oct: __magick_read__.cc
-	$(MKOCTFILE) $< `Magick++-config --cppflags` `Magick++-config --ldflags`
-
 clean: ; -$(RM) *.o octave-core core *.oct *~
deleted file mode 100644
--- a/src/__cordfltn__.cc
+++ /dev/null
@@ -1,289 +0,0 @@
-// Copyright (C) 2008 Soren Hauberg
-//
-// 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.
-//      -- Soren Hauberg, March 21st, 2008
-
-#include <octave/oct.h>
-
-// Template function for comparison
-// ET is the type of the matrix element
-template <class ET>
-inline bool compare (const ET a, const ET b)
-{
-    if (a > b)
-      return 1;
-    else
-      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 ();
-    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>
-ET selnth (ET *vals, int 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 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)
-{
-  const int ndims = dom.ndims ();
-  const octave_idx_type dom_numel = dom.numel ();
-  const dim_vector dom_size = dom.dims ();
-  const dim_vector A_size = A.dims ();
-
-  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);
-  const octave_idx_type out_numel = out.numel ();
-
-  // Iterate over every element of 'out'.
-  dim_vector idx_dim (ndims);
-  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);
-  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;
-   
-      for (octave_idx_type j = 0; j < dom_numel; j++)
-        {
-          for (int n = 0; n < ndims; n++)
-            A_idx (n) = out_idx (n) + dom_idx (n);
-       
-          if (dom (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);
-
-      // Prepare for next iteration
-      out.increment_index (out_idx, out_size);
-       
-      OCTAVE_QUIT;
-    }
-    
-  retval (0) = octave_value (out);
-    
-  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, , "\
--*- 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\
-@seealso{cordflt2, ordfilt2}\n\
-@end deftypefn\n\
-")
-{
-  octave_value_list retval;
-  if (args.length () != 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 ();
-  if (error_state)
-    {
-      error ("__cordfltn__: invalid input");
-      return retval;
-    }
-   
-  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");
-      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);
-    } 
-  else
-    {
-      error ("__cordfltn__: first input should be a real, complex, or integer array");
-    }
-
-  return retval;
-}
new file mode 100644
--- /dev/null
+++ b/src/__spatial_filtering__.cc
@@ -0,0 +1,578 @@
+// Copyright (C) 2008 Soren Hauberg
+//
+// 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.
+//      -- Soren Hauberg, March 21st, 2008
+
+#include <octave/oct.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 < 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 ();
+  const dim_vector A_size = A.dims ();
+
+  octave_idx_type len = 0;
+  for (octave_idx_type i = 0; i < dom_numel; i++)
+    len += dom (i);
+
+  dim_vector dim_offset (dom_size);
+  for (int i = 0; i < ndims; i++)
+    dim_offset (i) = (dom_size (i)+1) / 2 - 1;
+
+  // Allocate output
+  dim_vector out_size (dom_size);
+  for (int i = 0; i < ndims; i++)
+    out_size (i) = A_size (i) - dom_size (i) + 1;
+
+  MTout out = MTout (out_size);
+  const octave_idx_type out_numel = out.numel ();
+
+  // Iterate over every element of 'out'.
+  dim_vector idx_dim (ndims);
+  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
+      int l = 0;
+      for (int n = 0; n < ndims; n++)
+        dom_idx (n) = 0;
+   
+      for (octave_idx_type j = 0; j < dom_numel; j++)
+        {
+          for (int n = 0; n < ndims; n++)
+            A_idx (n) = out_idx (n) + dom_idx (n);
+       
+          if (dom (dom_idx))
+            values (l++) = A (A_idx) + S (dom_idx);
+       
+          dom.increment_index (dom_idx, dom_size);
+        }
+            
+      // Compute filter result
+      out (out_idx) = filter_function (values, len, arg4);
+
+      // Prepare for next iteration
+      out.increment_index (out_idx, out_size);
+       
+      OCTAVE_QUIT;
+    }
+    
+  retval (0) = octave_value (out);
+    
+  return retval;
+}
+
+/**
+ * The Octave function
+ */
+ 
+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 \"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;
+  const int nargin = args.length ();
+  if (nargin < 4)
+    {
+        print_usage ();
+        return retval;
+    }
+    
+  const boolNDArray dom = args (1).bool_array_value ();
+  if (error_state)
+    {
+      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 ("__spatial_filtering__: input must be of the same dimension");
+      return retval;
+    }
+  
+  
+  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 ("__spatial_filtering__: unknown method '%s'.", method.c_str ());
+    }
+
+  return retval;
+}