How to easily apply a Gauss filter to a list/array of doubles

1.3k Views Asked by At

My code works already, but IMHO ugly:

public IList<OxyPlot.DataPoint> Points { get; private set; }
public IList<OxyPlot.DataPoint> Points2 { get; private set; }

void ApplyGaussFilter()
{
    var gauss = MathNet.Numerics.Window.Gauss(5, 0.8);
    Points2 = new List<OxyPlot.DataPoint>();
    var inputArray = Points.ToArray();

    for (int i=2; i< inputArray.Length-2; ++i)
    {
        double smoothedVal = (inputArray[i - 2].Y * gauss[0] +
                      inputArray[i - 1].Y * gauss[1] +
                      inputArray[i].Y * gauss[2] +
                      inputArray[i + 1].Y * gauss[3] +
                      inputArray[i + 2].Y * gauss[4]) / gauss.Sum();

        Points2.Add(new DataPoint(inputArray[i].X, smoothedVal));
    }
}

Smoothed measurements

Is there an easier way? Additionally I don't know how to correctly handle the two outer most array values. And for changing the Gauss width from 5 to any other value, the for loop needs to be manually updated as well.

All examples I have found so far are based on images to filter these. E.G. guassian smoothening formula application.
Applying the filter with image libraries with two lines of code looks much nicer compared to my ugly for loop.

GaussianBlur filter = new GaussianBlur(4, 11);
filter.ApplyInPlace(graph);
1

There are 1 best solutions below

1
On

I tried to get MathNet.Numerics.Window.Gauss() to work, but I couldn't. Maybe I am doing something wrong, but I have found other solutions that work for me, for example this one:

private static double[] gaussianKernel1d( int kernelRadius, double sigma )
{
    double[] kernel = new double[kernelRadius + 1 + kernelRadius];
    for( int xx = -kernelRadius; xx <= kernelRadius; xx++ )
        kernel[kernelRadius + xx] = Math.Exp( -(xx * xx) / (2 * sigma * sigma) ) /
                (Math.PI * 2 * sigma * sigma);
    return kernel;
}

Once you have that, one-dimensional gaussian filtering can be achieved as follows:

double[]  array        = (your array)
double    sigma        = (your sigma)
int       kernelRadius = (int)Math.Ceiling( sigma * 2.57 ); // significant radius
double[]  kernel       = gaussianKernel1d( kernelRadius, sigma );
double[]  result       = filter( array, kernel );

static double[] filter( double[] array, double[] kernel )
{
    Assert( kernel.Length % 2 == 1 ); //kernel size must be odd.
    int       kernelRadius = kernel.Length / 2;
    int       width  = array.GetLength( 0 );
    double[] result = new double[width];
    for( int x = 0; x < width; x++ )
    {
        double sumOfValues  = 0;
        double sumOfWeights = 0;
        for( int i = -kernelRadius; i <= kernelRadius; i++ )
        {
            double value = array[clamp( x + i, 0, width - 1 )];
            double weight = kernel[kernelRadius + i];
            sumOfValues  += value * weight;
            sumOfWeights += weight;
        }
        result[x] = sumOfValues / sumOfWeights;
    }
    return result;
}

The clamp() function is this little thing:

private static int clamp( int value, int min, int max )
{
    if( value < min )
        return min;
    if( value > max )
        return max;
    return value;
}

and it ensures that the value at the boundary of your data will be used in place of any values that lie beyond that boundary. This is one way of handling what you are referring to as your "two outer most array values".