--- /dev/null
+## 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/>.
+##
+## This function is based on 'ordfilt2' by Teemu Ikonen which is released under
+## GPLv2 or later.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} ordfiltn(@var{A}, @var{nth}, @var{domain}, [@var{S}, @var{padding}])
+## Two dimensional ordered filtering.
+##
+## Ordered filter replaces an element of @var{A} with the @var{nth}
+## element of the sorted set of neighbours defined by the logical
+## (boolean) matrix @var{domain}.
+## Neighbour elements are selected to the sort if the corresponding
+## element in the @var{domain} matrix is true.
+##
+## The optional variable @var{S} is a matrix of size(@var{domain}).
+## Values of @var{S} corresponding to nonzero values of domain are
+## added to values obtained from @var{A} when doing the sorting.
+##
+## Optional variable @var{padding} determines how the matrix @var{A}
+## is padded from the edges. See @code{padarray} for details.
+##
+## @seealso{ordfilt2, padarray}
+## @end deftypefn
+
+
+## Author: Teemu Ikonen <tpikonen@pcu.helsinki.fi>
+## Created: 5.5.2000
+## Keywords: image processing filtering
+
+function retval = ordfiltn(A, nth, domain, varargin)
+ ## Check input
+ if (nargin < 3)
+ error("ordfiltn: not enough input arguments");
+ endif
+ if (!ismatrix(A))
+ error("ordfiltn: first input must be an array");
+ endif
+ if (!isscalar(nth) || nth <= 0 || nth != round(nth))
+ error("ordfiltn: second input argument must be a positive integer");
+ endif
+ if (!ismatrix(domain) && !isscalar(domain))
+ error("ordfiltn: third input argument must be an array or a scalar");
+ endif
+ if (isscalar(domain) && (domain <= 0 || domain != round(domain)))
+ error("ordfiltn: third input argument must be a positive integer, when it is a scalar");
+ endif
+ if (isscalar(domain))
+ domain = ones(repmat(domain, 1, ndims(A)), "logical");
+ endif
+
+ if (ndims(A) != ndims(domain))
+ error("ordfiltn: first and second argument must have same dimensionality");
+ endif
+ if (any(size(A) < size(domain)))
+ error("ordfiltn: domain array cannot be larger than the data array");
+ endif
+
+ ## Parse varargin
+ S = zeros(size(domain));
+ padding = 0;
+ for i=1:length(varargin)
+ a = varargin{:};
+ if (ischar(a) || isscalar(a))
+ padding = a;
+ elseif (ismatrix(a) && size_equal(a, domain))
+ S = a;
+ endif
+ endfor
+
+ ## Make sure 'domain' is logical. The C++ code assumes this.
+ domain = logical(domain);
+
+ ## Pad array
+ pad = floor(size(domain)/2);
+ A = padarray(A, pad, padding);
+ even = ( round(size(domain)/2) == size(domain)/2 );
+ idx = cell(1, ndims(A));
+ for k = 1:ndims(A)
+ idx{k} = (even(k)+1):size(A,k);
+ endfor
+ A = A(idx{:});
+
+ ## Perform the filtering
+ retval = __spatial_filtering__ (A, domain, "ordered", S, nth);
+endfunction