TimedomainMedianFilter

TimedomainMedianFilter.hpp

#pragma once

#include "ofMain.h"

template <typename T>
class TimedomainMedianFilter
{
public:
  
  TimedomainMedianFilter() : num_frame_buffer(5), current_frame_index(0) {}
	
	void setup(int w, int h)
	{
		pixels.resize(w * h * num_frame_buffer, 0);
		result.allocate(w, h, OF_IMAGE_GRAYSCALE);
		result.set(0);
	}
	
	const ofPixels_<T>& update(ofPixels_<T> &pix)
	{
		const int W = pix.getWidth();
		const int H = pix.getHeight();
		const int N = W * H;
		
		T *pixels_ptr = &pixels[0];
		
		memcpy(pixels_ptr + current_frame_index * N, pix.getPixels(), N);
		
		{
			T *result_ptr = result.getPixels();
			T median[num_frame_buffer];
			
			for (int n = 0; n < N; n++)
			{
				T *ptr = pixels_ptr + n;
				
				for (int i = 0; i < num_frame_buffer; i++)
				{
					median[i] = *ptr;
					ptr += N;
				}
				
				// result_ptr[n] = opt_med3(median);
				result_ptr[n] = opt_med5(median);
				// result_ptr[n] = torben(median, num_frame_buffer);
			}
		}
		
		current_frame_index++;
		current_frame_index %= num_frame_buffer;
		
		return result;
	}
	
protected:
	
	const int num_frame_buffer;
	int current_frame_index;
	
	vector<T> pixels;
	ofPixels_<T> result;

	// via http://ndevilla.free.fr/median/median.pdf
	
#define PIX_SORT(a,b) { if ((a)>(b)) std::swap((a),(b)); }
	
	inline static T opt_med3(T * p)
	{
		PIX_SORT(p[0],p[1]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[0],p[1]) ;
		return(p[1]) ;
	}
	
	inline static T opt_med5(T * p)
	{
		PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ;
		PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ;
		PIX_SORT(p[1],p[2]) ; return(p[2]) ;
	}
	
#undef PIX_SORT
	
	inline static T torben(T *m, int n)
	{
		int i, less, greater, equal;
		T min, max, guess, maxltguess, mingtguess;
		min = max = m[0] ;
		for (i=1 ; i<n ; i++) {
			if (m[i]<min) min=m[i];
			if (m[i]>max) max=m[i];
		}
		while (1) {
			guess = (min+max)/2;
			less = 0; greater = 0; equal = 0;
			maxltguess = min ;
			mingtguess = max ;
			for (i=0; i<n; i++) {
				if (m[i]<guess) {
					less++;
					if (m[i]>maxltguess) maxltguess = m[i] ;
				} else if (m[i]>guess) {
					greater++;
					if (m[i]<mingtguess) mingtguess = m[i] ;
				} else equal++;
			}
			if (less <= (n+1)/2 && greater <= (n+1)/2) break ;
			else if (less>greater) max = maxltguess ;
			else min = mingtguess;
		}
		if (less >= (n+1)/2) return maxltguess;
		else if (less+equal >= (n+1)/2) return guess;
		else return mingtguess;
	}
};