medianfilter.c File Reference

Core MEX routine implementing a fast version of the classical 1D running median filter of size W, where W is odd. More...

#include "mdb.h"

Go to the source code of this file.

Macros
#define	SWAP(a, b)

Functions
double	quickSelect (double *arr, int n)
void	median_filter (double *x, double *m, long n, long W)
	Applies a median filter to an input signal.

Detailed Description

Core MEX routine implementing a fast version of the classical 1D running median filter of size W, where W is odd.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M.A. Little, N.S. Jones, H. Shang, R. Soliday

Definition in file medianfilter.c.

Macro Definition Documentation

SWAP

#define SWAP	(		a,
			b )

Value:

  {                          \
    register double t = (a); \
    (a) = (b);               \
    (b) = t;                 \
  }

Definition at line 78 of file medianfilter.c.

#define SWAP(a, b)           \
  {                          \
    register double t = (a); \
    (a) = (b);               \
    (b) = t;                 \
  }

Function Documentation

median_filter()

void median_filter	(	double *	x,
		double *	m,
		long	n,
		long	W )

Applies a median filter to an input signal.

This function processes the input signal using a sliding window approach to compute the median value for each position. It handles boundary conditions by replicating the edge values to maintain the window size.

Parameters

x	Input signal array.
m	Output signal buffer array where the median filtered signal will be stored.
n	Size of the input signal.
W	Size of the sliding window (must be an odd number, W = 2*W2 + 1).

Definition at line 157 of file medianfilter.c.

                                                         {
  /* x  -- input signal
     m  -- output signal buffer 
     n  -- size of input signal
     W  -- size of slding window (must be odd number) W = 2*W2 + 1) */
  long i, k, idx, W2;
  double *w;
  if (W % 2 == 0)
    W = W + 1;
  W2 = (W - 1) / 2;
 
  w = calloc(sizeof(*w), W);
  for (i = 0; i < n; i++) {
    for (k = 0; k < W; k++) {
      idx = i - W2 + k;
      if (idx < 0) {
        w[k] = x[0];
      } else if (idx >= n) {
        w[k] = x[n - 1];
      } else {
        w[k] = x[idx];
      }
    }
    m[i] = quickSelect(w, W);
  }
  free(w);
}

quickSelect()

double quickSelect	(	double *	arr,
		int	n )

Definition at line 85 of file medianfilter.c.

                                       {
  long low, high;
  long median;
  long middle, ll, hh;
 
  low = 0;
  high = n - 1;
  median = (low + high) / 2;
  while (1) {
    /* One element only */
    if (high <= low)
      return arr[median];
 
    /* Two elements only */
    if (high == low + 1) {
      if (arr[low] > arr[high])
        SWAP(arr[low], arr[high]);
      return arr[median];
    }
 
    /* Find median of low, middle and high items; swap to low position */
    middle = (low + high) / 2;
    if (arr[middle] > arr[high])
      SWAP(arr[middle], arr[high]);
    if (arr[low] > arr[high])
      SWAP(arr[low], arr[high]);
    if (arr[middle] > arr[low])
      SWAP(arr[middle], arr[low]);
 
    /* Swap low item (now in position middle) into position (low+1) */
    SWAP(arr[middle], arr[low + 1]);
 
    /* Work from each end towards middle, swapping items when stuck */
    ll = low + 1;
    hh = high;
    while (1) {
      do {
        ll++;
      } while (arr[low] > arr[ll]);
 
      do {
        hh--;
      } while (arr[hh] > arr[low]);
 
      if (hh < ll)
        break;
 
      SWAP(arr[ll], arr[hh]);
    }
 
    /* Swap middle item (in position low) back into correct position */
    SWAP(arr[low], arr[hh]);
 
    /* Reset active partition */
    if (hh <= median)
      low = ll;
    if (hh >= median)
      high = hh - 1;
  }
}