new file mode 100644
--- /dev/null
+++ b/inst/cordflt2.m
@@ -0,0 +1,27 @@
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} cordflt2(@var{A}, @var{nth}, @var{domain}, @var{S})
+## Implementation of two-dimensional ordered filtering. This function has been
+## deprecated and should NOT be used. Instead use @code{ordfilt2}.
+## @seealso{ordfilt2}
+## @end deftypefn
+
+function varargout = cordflt2(varargin)
+  warning(["cordflt2: this function is deprecated and will be removed in upcoming "
+           "releases. Use 'ordfilt2' instead."]);
+  [varargout{1:nargout}] = __cordfltn__(varargin{:});
+endfunction

--- a/src/Makefile
+++ b/src/Makefile
@@ -12,7 +12,7 @@
 	IMAGEMAGICK=__magick_read__.oct
 endif
 
-all: cordflt2.oct __cordfltn__.oct bwlabel.oct bwfill.oct rotate_scale.oct \
+all: __cordfltn__.oct bwlabel.oct bwfill.oct rotate_scale.oct \
 	hough_line.oct graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \
 	$(JPEG) $(PNG) $(IMAGEMAGICK)
 

deleted file mode 100644
--- a/src/cordflt2.cc
+++ /dev/null
@@ -1,211 +0,0 @@
-// Copyright (C) 2000 Teemu Ikonen
-//
-// 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 2
-// 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/>.
-
-#include <octave/oct.h>
-
-#define SWAP(a, b) { SWAP_temp = (a); (a)=(b); (b) = SWAP_temp; }
-
-// 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)
-{
-    double anorm2 = a.real() * a.real() + a.imag() * a.imag();
-    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 SWAP_temp;
-    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])) {
-		SWAP(vals[l], vals[r]);
-	    }
-	    return vals[nth];
-	} else {
-	    mid = (l+r) >> 1;
-	    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])) {
-		SWAP(vals[l], vals[r]);
-	    }
-	    if(compare<ET>(vals[l+1], vals[r])) {
-		SWAP(vals[l+1], vals[r]);
-	    }
-	    if(compare<ET>(vals[l], vals[l+1])) {
-		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;
-		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, boolMatrix dom, MT S)
-{
-    int i, j, c, d;
-    
-    int len = 0;
-    for(j = 0; j < dom.columns(); j++) {
-	for(i = 0; i < dom.rows(); i++) {
-	    if(dom.elem(i,j)) 
-	      len++;
-	}
-    }
-    if(nth > len - 1) {
-	warning("nth should be less than number of non-zero values in domain");
-	warning("setting nth to largest possible value\n");
-	nth = len - 1;
-    }
-    if(nth < 0) {
-	warning("nth should be non-negative, setting to 1\n");
-	nth = 0; // nth is a c-index
-    }
-                
-    int rowoffset = (dom.columns() + 1)/2 - 1;
-    int coloffset = (dom.rows() + 1)/2 - 1;
-
-    //outputs
-    octave_value_list out;
-    const int origx = A.columns() - dom.columns()+1;
-    const int origy = A.rows() - dom.rows()+1;
-    MT retval = MT(origy, origx);
-
-    int *offsets = new int[len];
-    ET *values = new ET[len];
-    ET *adds = new ET[len];
-    
-    c = 0;
-    d = A.rows();
-    for(j = 0; j < dom.columns(); j++) {
-	for(i = 0; i < dom.rows(); i++) {
-	    if(dom.elem(i,j)) {
-		offsets[c] = (i - coloffset) + (j - rowoffset)*d;
-		adds[c] = S.elem(i,j);
-		c++;
-	    }
-	}
-    }
-    
-    ET *data = A.fortran_vec();
-    int base = coloffset + A.rows()*rowoffset;
-    for(j = 0; j < retval.columns(); j++) {
-	for(i = 0; i < retval.rows(); i++) {
-	    for(c = 0; c < len; c++) {
-		values[c] = data[base + offsets[c]] + adds[c];
-	    }
-	    base++;
-	    retval(i, j) = selnth(values, len, nth);
-	}
-	base += dom.rows() - 1;
-    }
-
-    out(0) = octave_value(retval);
-    
-    return out;
-}
-
-// instantiate template functions
-template bool compare<double>(const double, const double);
-template double selnth(double *, int, int);
-template Complex selnth(Complex *, int, int);
-template octave_value_list do_filtering<Matrix, double>(Matrix, int, boolMatrix, Matrix);
-// g++ is broken, explicit instantiation of specialized template function
-// confuses the compiler.
-//template int compare<Complex>(const Complex, const Complex);
-template octave_value_list do_filtering<ComplexMatrix, Complex>(ComplexMatrix, int, boolMatrix, ComplexMatrix);
-
-DEFUN_DLD(cordflt2, args, , "\
--*- texinfo -*-\n\
-@deftypefn {Function File} cordflt2(@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{ordfilt2}\n\
-@end deftypefn\n\
-")
-{
-    if(args.length() != 4) {
-	print_usage ();
-	return octave_value_list();
-    }
-    
-    // nth is an index to an array, thus - 1
-    int nth = (int) (args(1).vector_value())(0) - 1;
-    boolMatrix dom = args(2).bool_matrix_value();
-
-    octave_value_list retval;
-    
-    if(args(0).is_real_matrix()) {
-	Matrix A = args(0).matrix_value();
-	Matrix S = args(3).matrix_value();
-	retval = do_filtering<Matrix, double>(A, nth, dom, S);
-    } 
-    else if(args(0).is_complex_matrix()) {
-	ComplexMatrix A = args(0).complex_matrix_value();
-	ComplexMatrix S = args(3).complex_matrix_value();
-	retval = do_filtering<ComplexMatrix, Complex>(A, nth, dom, S);
-    } 
-    else {
-	error("A should be real or complex matrix\n");
-	return octave_value_list();
-    }
-    
-    return retval;
-     
-}