/tags/R2009-05-08/main/image/src/__spatial_filtering__.cc
C++ | 578 lines | 445 code | 71 blank | 62 comment | 173 complexity | 3aaf9a2f2894048a5b3e16f9aa77da60 MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause, LGPL-2.1, GPL-3.0, LGPL-3.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>
- /**
- * 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 * xlog2 (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;
- }