From octave-sources-request at bevo dot che dot wisc dot edu  Wed Sep  8 07:55:18 1999
Subject: Median filtering
From: Olli Saarela <olli dot saarela at kcl dot fi>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Wed, 08 Sep 1999 15:54:25 +0300

This is a multi-part message in MIME format.

--------------7D99217968F6
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Median filtering is an important operation in signal processing. An *.m
implementation requires each overlapping data window to be sorted
separately, which makes calculation relatively inefficient.

As far as I know, the enclosed medfilt1.cc (one-dimensional median
filtering) implementation is compatible with the medfilt1.m in the
Matlab signal processing toolbox, with the exception of values for
signal edges (values before and after the scope of the input signal).
This implementation assumes that the values of the first and the last
data points in the signal are repeated before and after
(correspondingly) the data sequence available. In the Matlab version,
zero values are used. (Any desired edge handling can be easily
implemented in an *.m wrapper.)

Please consider this implementation for inclusion in future releases of
Octave. This is the first time I have written any C++ code for Octave,
so I appreciate any comments you might have.

Best regards,
	Olli
-- 
Olli Saarela                       KCL Development Oy
Olli dot Saarela at kcl dot fi                tel. +358-9-4371538 (office)
Tekniikantie 2, Espoo-Otaniemi     fax. +358-9-464305
P.O. Box 70, FIN-02151 Espoo, Finland

--------------7D99217968F6
Content-Type: text/plain; charset=us-ascii; name="medfilt1.cc"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline; filename="medfilt1.cc"

/*

Copyright (C) 1999 Olli Saarela

This file is not (yet?) part of Octave.

Octave 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, or (at your option) any
later version.

Octave 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 Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

// medfilt1: one-dimensional median filtering
// (like moving average, but median instead of average)

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "defun-dld.h"
#include "error.h"
#include "oct-obj.h"
#include "help.h"
#include "utils.h"
#include "mappers.h"
#include "gripes.h"

typedef int (*comparefun) (const void *, const void *);

static int
medfilt_compare(const double *x, const double *y) {
  // NaN values are sorted to the end of the array
  // (test for NaN values before the inequalities because
  // some compilers have difficulties with NaN inequalities)
  if (xisnan(*x))
    return xisnan(*y) ? 0 : 1;
  if (xisnan(*y))
    return -1;

  if (*x > *y)
    return 1;
  if (*y > *x)
    return -1;
  return 0;
}

static int
medfilt_compare(const Complex *x, const Complex *y) {
  double absx = std::abs(*x);
  double absy = std::abs(*y);
  return medfilt_compare(&absx, &absy);
}

static inline double*
medfilt_search(double *window, const long wlen, const double out) {
  // select the proper overload
  int (*fun) (const double *, const double *) = medfilt_compare;
  return (double*)bsearch(&out, window, wlen, sizeof(double), (comparefun)fun);
}

static inline Complex*
medfilt_search(Complex *window, const long wlen, const Complex out) {
  // ambiguous sort order for complex values 
  // - based on abs()
  // - no proper handling for "semi-NaN" values, e.g., (3,NaN) vs. (NaN,3)
  // ---> perform a full search
  for (long j=0; j<wlen; j++)
    if (0 == memcmp(&out, window+j, sizeof(Complex)))
      return window+j;
  return NULL;
}

template <class T>
static void
medfilt_update_sorted(T *window, const long wlen, const T out, const T in) {
  T *inp, *outp, *endp;

  outp = medfilt_search(window, wlen, out);
  assert (outp != NULL);
  inp  = outp;

  switch (medfilt_compare(&in, &out)) {
  case 0:
    // replace in place
    break;

  case 1:
    // replace after
    endp = window + wlen;
    while (inp+1 < endp && medfilt_compare(&in, inp+1) == 1) {
      *inp = *(inp+1);
      inp++;
    }
    break;

  case -1:
    // replace before
    while (inp > window && medfilt_compare(&in, inp-1) == -1) {
      *inp = *(inp-1);
      inp--;
    }

  }

  *inp = in;
}

template static void
medfilt_update_sorted(double *window, const long wlen, const double out, const double in);

template static void
medfilt_update_sorted(Complex *window, const long wlen, const Complex out, const Complex in);

template <class T>
static void
medfilt1_filter(const T *inarr, T *outarr, const long arrlen, const long wlen) {
  int     odd = wlen%2;
  long    j, wlen21, inj, outj;
  T       *window = new T [wlen];
  int (*fun) (const T *, const T *) = medfilt_compare; // select the proper overload

  if (odd)
    wlen21 = (wlen-1)/2; // odd
  else
    wlen21 = wlen/2; // even

  // initialise the median window
  inj = -wlen21-1;
  for (j=0; j<wlen; j++) {
    if (inj < 0)
      window[j] = inarr[0];
    else if (inj >= arrlen)
      window[j] = inarr[arrlen-1];
    else
      window[j] = inarr[inj];
    inj++;
  }
  qsort(window, wlen, sizeof(T), (comparefun)fun);
  
  // filtering
  if (odd)
    inj = wlen21; // odd
  else
    inj = wlen21-1; // even
  outj = -wlen21-1;
  for (j=0; j<arrlen; j++) {
    medfilt_update_sorted(window, wlen,
                          (outj<0) ? inarr[0] : inarr[outj],
                          (inj<arrlen) ? inarr[inj] : inarr[arrlen-1]);
    // a NaN value pollutes the output
    if (xisnan(window[wlen-1]))
      outarr[j] = window[wlen-1];
    else if (odd)
      outarr[j] = window[wlen21];
    else
      outarr[j] = (window[wlen21-1] + window[wlen21]) / 2.0;
    //
    inj++;
    outj++;
  }
  delete [] window;
}

template static void
medfilt1_filter(const double *inarr, double *outarr, const long arrlen, const long wlen);

template static void
medfilt1_filter(const Complex *inarr, Complex *outarr, const long arrlen, const long wlen);


static void
medfilt1(const Matrix &x, Matrix &y, const long wlen) {
  long nr = x.rows();
  long nc = x.columns();

  if (nr == 1 || nc == 1) {
    double *py = y.fortran_vec ();
    const double *px = x.data ();

    medfilt1_filter(px, py, nr*nc, wlen);

  } else {
    ColumnVector ycol(nr);
    double *py = y.fortran_vec ();

    for (long j=0; j<nc; j++) {
      const ColumnVector xcol = x.column(j);
      const double *px = xcol.data ();

      medfilt1_filter(px, py, nr, wlen);
      y.insert(ycol, 0, j);
    }
  }
}

static void
medfilt1(const ComplexMatrix &x, ComplexMatrix &y, const long wlen) {
  long nr = x.rows();
  long nc = x.columns();

  if (nr == 1 || nc == 1) {
    Complex *py = y.fortran_vec ();
    const Complex *px = x.data ();

    medfilt1_filter(px, py, nr*nc, wlen);

  } else {
    ComplexColumnVector ycol(nr);
    Complex *py = y.fortran_vec ();

    for (long j=0; j<nc; j++) {
      const ComplexColumnVector xcol = x.column(j);
      const Complex *px = xcol.data ();

      medfilt1_filter(px, py, nr, wlen);
      y.insert(ycol, 0, j);
    }
  }
}

DEFUN_DLD (medfilt1, args, nargout,
  "medfilt1 (X, N): Median filter\n\
\n\
If X is a vector, it is filtered with moving median of length N.\n\
If X is a matrix, each column is filtered separately. The value\n\
returned is the same size as X. The default value for N is 3.\n\
\n\
For the edge points, constant values X(1) and X(length(X)) are\n\
assumed before and after the scope of X, correspondingly.\n\
\n\
An optional third argument is accepted (and ignored) for\n
compatibility with other implementations.\n")
{
  long wlen = 3;
  octave_value retval;

  int nargin  = args.length ();

  // Allow a dummy third argument for compatibility with some
  // m-file implementations (which use that argument for
  // performance tuning between time and space)
  if (nargin < 1 || nargin > 3 || nargout > 1) {
    print_usage ("medfilt1");
    return retval;
  }

  // Filter window length
  if (nargin == 1) {
    wlen = 3;
  } else if (args(1).is_real_scalar ()) {
    if (xisnan(args(1).double_value ()))
      error ("medfilt1: the median window length must be > 0");
    wlen = NINT(args(1).double_value());
  } else {
    print_usage ("medfilt1");
    return retval;
  }
  if (wlen < 1) {
    error ("medfilt1: the median window length must be > 0");
    return retval;
  }

  // Data to be filtered
  if (!(args(0).is_numeric_type() || args(0).is_range())) {
    gripe_wrong_type_arg ("medfilt1", args(0));

  } else if (wlen <= 1 || args(0).is_scalar_type() || args(0).is_empty()) {
    retval = args(0);
    retval.make_unique();

  } else if (args(0).is_real_type()) {
    Matrix x = args(0).matrix_value();
    if (error_state) {
      error ("medfilt1: the first argument must be a vector or a matrix");
    } else {
      Matrix y (x.rows(), x.columns());
      medfilt1(x, y, wlen);
      retval = y;
    }

  } else if (args(0).is_complex_type()) {
    ComplexMatrix x = args(0).complex_matrix_value();
    if (error_state) {
      error ("medfilt1: the first argument must be a vector or a matrix");
    } else {
      ComplexMatrix y (x.rows(), x.columns());
      medfilt1(x, y, wlen);
      retval = y;
    }

  } else {
    gripe_wrong_type_arg ("medfilt1", args(0));
  }

  return retval;
}


/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/  

--------------7D99217968F6--