Filter bank taking too long for kernels to be created

Question

Take a look at this question.

public partial class Form1 : Form
{
    public Form1()
    {
        InitializeComponent();

        Bitmap image = DataConverter2d.ReadGray(StandardImage.LenaGray);
        Array2d<double> dImage = DataConverter2d.ToDouble(image);            

        int newWidth = Tools.ToNextPowerOfTwo(dImage.Width) * 2;
        int newHeight = Tools.ToNextPowerOfTwo(dImage.Height) * 2;

        double n = 6;
        double f0 = 0.5;
        double theta = 60;
        double a = 0.4;

        Array2d<Complex> kernel2d = CustomFft(newWidth, newHeight, n, f0, theta, a);

        dImage.PadTo(newWidth, newHeight);
        Array2d<Complex> cImage = DataConverter2d.ToComplex(dImage);
        Array2d<Complex> fImage = FourierTransform.ForwardFft(cImage);

        // FFT convolution .................................................
        Array2d<Complex> fOutput = new Array2d<Complex>(newWidth, newHeight);
        for (int x = 0; x < newWidth; x++)
        {
            for (int y = 0; y < newHeight; y++)
            {
                fOutput[x, y] = fImage[x, y] * kernel2d[x, y];
            }
        }

        Array2d<Complex> cOutput = FourierTransform.InverseFft(fOutput);
        Array2d<double> dOutput = Rescale2d.Rescale(DataConverter2d.ToDouble(cOutput));

        dOutput.CropBy((newWidth - image.Width) / 2, (newHeight - image.Height) / 2);

        Bitmap output = DataConverter2d.ToBitmap(dOutput, image.PixelFormat);

        Array2d<Complex> cKernel = FourierTransform.InverseFft(kernel2d);
        cKernel = FourierTransform.RemoveFFTShift(cKernel);
        Array2d<double> dKernel = DataConverter2d.ToDouble(cKernel);
        Array2d<double> dLimitedKernel = Rescale2d.Limit(dKernel);

        Bitmap kernel = DataConverter2d.ToBitmap(dLimitedKernel, image.PixelFormat);

        pictureBox1.Image = image;
        pictureBox2.Image = kernel;
        pictureBox3.Image = output;
    }

    private double Basic(double u, double v, double n, double f0, double rad, double a, double b)
    {
        double ua = u + f0 * Math.Cos(rad);
        double va = v + f0 * Math.Sin(rad);

        double ut = ua * Math.Cos(rad) + va * Math.Sin(rad);
        double vt = (-1) * ua * Math.Sin(rad) + va * Math.Cos(rad);

        double p = ut/a;
        double q = vt/b;

        double sqrt = Math.Sqrt(p*p + q*q);

        return 1.0 / (1.0+ 0.414 * Math.Pow(sqrt, 2*n));
    }

    private double Custom(double u, double v, double n, double f0, double theta, double a)
    {
        double rad1 = (Math.PI / 180) * (90-theta);
        double rad2 = rad1 + Math.PI;
        double b = (a / 5.0) / (2*n);

        double ka = Basic(u, v, n, f0, rad1, a, b);
        double kb = Basic(u, v, n, f0, rad2, a, b);

        return Math.Max(ka, kb);
    }

    private Array2d<Complex> CustomFft(double sizeX, double sizeY, double n, double f0, double theta, double a)
    {
        double halfX = sizeX / 2;
        double halfY = sizeY / 2;

        Array2d<Complex> kernel = new Array2d<Complex>((int)sizeX, (int)sizeY);

        for (double y = 0; y < sizeY; y++)
        {
            double v = y / sizeY;

            for (double x = 0; x < sizeX; x++)
            {
                double u = x / sizeX;

                double kw = Custom(u, v, n, f0, theta, a);

                kernel[(int)x, (int)y] = new Complex(kw, 0);
            }
        }

        return kernel;
    }
}

I created a 40-kernel filter bank where each kernel would have same size but different parameters.

This bank of filter is taking around 55000 milliseconds for the kernels to be generated.

How can I reduce this time?

Note: the following line is taking the most time:

Array2d<Complex> kernel2d = CustomFft(newWidth, newHeight, n, f0, theta, a);

.

Example source code for creating a kernel bank:

public partial class Form1 : Form
{
    public Form1()
    {
        InitializeComponent();

        Stopwatch sw = new Stopwatch();
        sw.Start();

        Bitmap image = DataConverter2d.Read(StandardImage.Scratch1);

        int newWidth = Tools.ToNextPowerOfTwo(image.Width);
        int newHeight = Tools.ToNextPowerOfTwo(image.Height);

        List<FrequencyDomainKernel> list = new List<FrequencyDomainKernel>();
        for (int i = 0; i < 40; i++)
        {          
            FrequencyDomainKernel k= new FrequencyDomainKernel();
            k.Size = new System.Drawing.Size(newWidth, newHeight);
            k.theta = i;
            k.Compute();

            list.Add(k);
        }

        sw.Stop();

        MessageBox.Show(sw.ElapsedMilliseconds.ToString());
    }
}

^{Note: this is not an executable code. This code is just to give an idea regarding my method of generating a kernel.}

Answer 1

As a starting point, try this:

private Array2d<Complex> CustomFft( double sizeX, double sizeY, double n, double f0, double theta,
                                    double a ) {
    double halfX = sizeX / 2;
    double halfY = sizeY / 2;

    //pre-calculated values
    double rad1, rad2, b, f1cos, f1sin, f2cos, f2sin, cosrad1, sinrad1, cosrad2, sinrad2;
    rad1 = ( Math.PI / 180 ) * ( 90 - theta );
    rad2 = rad1 + Math.PI;
    cosrad1 = Math.Cos( rad1 );
    sinrad1 = Math.Sin( rad1 );
    cosrad2 = Math.Cos( rad2 );
    sinrad2 = Math.Sin( rad2 );
    f1cos = f0 * Math.Cos( rad1 );
    f1sin = f0 * Math.Sin( rad1 );
    f2cos = f0 * Math.Cos( rad2 );
    f2sin = f0 * Math.Sin( rad2 );
    b = ( a / 5.0 ) / ( 2 * n );

    Array2d<Complex> kernel = new Array2d<Complex>((int)sizeX, (int)sizeY);

    /*
    for( double y = 0; y < sizeY; y++ ) {
        double v = y / sizeY;

        for( double x = 0; x < sizeX; x++ ) {
            double u = x / sizeX;

            double kw = Custom2( u, v, n, f0, theta, a, rad1, rad2, b, f1cos, f1sin, f2cos, f2sin,
                                 cosrad1, sinrad1 , cosrad2, sinrad2 );

            kernel[ (int)x, (int)y ] = new Complex( kw, 0 );
        }
    }*/

    //reverse the loops
    for( double x = 0; x < sizeX; x++ ) {    
        double u = x / sizeX;

        for( double y = 0; y < sizeY; y++ ) {
            double v = y / sizeY;

            double kw = Custom2( u, v, n, f0, theta, a, rad1, rad2, b, f1cos, f1sin, f2cos, f2sin,
                                  cosrad1, sinrad1, cosrad2, sinrad2 );

            kernel[ (int)x, (int)y ] = new Complex( kw, 0 );
        }
    }

    return kernel;
}

private double Custom2( double u, double v, double n, double f0, double theta, double a,
                        double rad1, double rad2, double b, double f1cos, double f1sin, 
                        double f2cos, double f2sin, double cosrad1, double sinrad1,
                        double cosrad2, double sinrad2 ) {

    double ka = Basic2( u, v, n, f0, rad1, a, b, f1cos, f1sin, cosrad1, sinrad1 );
    double kb = Basic2( u, v, n, f0, rad2, a, b, f2cos, f2sin, cosrad2, sinrad2 );

    return ( ka >= kb ? ka : kb); //Math.Max( ka, kb );
}

private double Basic2( double u, double v, double n, double f0, double rad, double a, double b,
                       double fcos, double fsin, double cosrad, double sinrad ) {
    double ua = u + fcos;
    double va = v + fsin;

    double ut = ua * cosrad + va * sinrad;
    double vt = ( -1 ) * ua * sinrad + va * cosrad;

    double p = ut / a;
    double q = vt / b;

    double sqrt = Math.Sqrt( p * p + q * q );

    return 1.0 / ( 1.0 + 0.414 * Math.Pow( sqrt, 2 * n ) );
}

I am going to work on it a little more.

EDIT 1

My first observation is that you are using too many classes and in these classes too many methods. The problem with OOP is that sometimes the code is divided into too many classes with very tiny methods which in turn can make the code hard to follow (going through all the classes etc..) and unnecessarily slow for simple things. Let use some numbers:

A simple:

Array2d<Complex> cImage = DataConverter2d.ToComplex( dImage );

takes arround 65msec! It just subtitutes the real part of cImage with the value of dImage and the imaginary is set to zero. This trully shouldn't take that much but as I said OOP and classes. Lets do a simple re-configuration and instead of Array2d<Complex> use a simple double array:

double[] cImagefast = new double[ newWidth * newHeight * 2 ];

Thats it. It stores real-imaginary parts linearly eg

cImagefast[0] -> real part of cImage[0,0]
cImagefast[1] -> imag part of cImage[0,0]
cImagefast[2] -> real part of cImage[0,1]
cImagefast[3] -> imag part of cImage[0,1]
....

Lets re-write the code and test it:

Stopwatch sw = new Stopwatch();
Bitmap image = DataConverter2d.ReadGray( "Scratch1.jpg" );
Array2d<double> dImage = DataConverter2d.ToDouble( image );

int newWidth = Tools.ToNextPowerOfTwo( dImage.Width );
int newHeight = Tools.ToNextPowerOfTwo( dImage.Height );

dImage.PadTo( newWidth, newHeight ); //~60msec

//Array2d<Complex> cImage = DataConverter2d.ToComplex( dImage ); //~65msec

double[] dImagefast = new double[ newWidth * newHeight ];
double[] cImagefast = new double[ newWidth * newHeight * 2 ];
int count = newWidth * newHeight * 2, pos = 0, pos_ = 0, i;

// fill it with true values before testing
for( int y = 0; y < newHeight; y++ ) {
    for( int x = 0; x < newWidth; x++ ) {
        dImagefast[ pos ] = dImage[ x, y ];

        pos++;
    }
}

//Start the testing!
sw.Start();

for( pos = 0; pos < count; pos += 2 ){
    cImagefast[ pos ] = dImagefast[ pos_ ];

    pos_++;
}

sw.Stop();

MessageBox.Show( sw.ElapsedMilliseconds.ToString() );

Time needed ~10msec. More than 6 times faster!!!!

Lets do another one. Lets test FFT convolution which takes ~105msec

// FFT convolution ................................................. ~105msec
Array2d<Complex> fOutput = new Array2d<Complex>( newWidth, newHeight );
for( int x = 0; x < newWidth; x++ ) {
    for( int y = 0; y < newHeight; y++ ) {
        fOutput[ x, y ] = fImage[ x, y ] * kernel2d[ x, y ];
    }
}

Re-write it with simple 1 dimension double arrays:

Stopwatch sw = new Stopwatch();
Bitmap image = DataConverter2d.ReadGray( "Scratch1.jpg" );
Array2d<double> dImage = DataConverter2d.ToDouble( image );
int newWidth = Tools.ToNextPowerOfTwo( dImage.Width );
int newHeight = Tools.ToNextPowerOfTwo( dImage.Height );
double n = 6, f0 = 0.5, theta = 60, a = 0.4;

Array2d<Complex> kernel2d = CustomFft( newWidth, newHeight, n, f0, theta, a );

dImage.PadTo( newWidth, newHeight ); //~60msec
Array2d<Complex> cImage = DataConverter2d.ToComplex( dImage ); //~65msec
Array2d<Complex> fImage = FourierTransform.ForwardFft( cImage );

/*
// FFT convolution ............................................. ~105msec
Array2d<Complex> fOutput = new Array2d<Complex>( newWidth, newHeight );
for( int x = 0; x < newWidth; x++ ) {
    for( int y = 0; y < newHeight; y++ ) {
        fOutput[ x, y ] = fImage[ x, y ] * kernel2d[ x, y ];
    }
}
*/

// Fast FFT convolution ..................................... ~31msec!!!
int pos = 0, count = newWidth * newHeight * 2;
double[] fImagefast = new double[ newWidth * newHeight * 2 ];
double[] kernel2dfast = new double[ newWidth * newHeight * 2 ];

//lets fill the arrays with values
for( int y = 0; y < newHeight; y++ ) {
    for( int x = 0; x < newWidth; x++ ) {
        kernel2dfast[ pos ] = kernel2d[ x, y ].Real;
        kernel2dfast[ pos + 1 ] = kernel2d[ x, y ].Imaginary;

        pos += 2;
    }
}

pos = 0;
for( int y = 0; y < newHeight; y++ ) {
    for( int x = 0; x < newWidth; x++ ) {
        fImagefast[ pos ] = fImage[ x, y ].Real;
        fImagefast[ pos + 1 ] = fImage[ x, y ].Imaginary;

        pos += 2;
    }
}

//Test!
sw.Start();
double[] fOutputfast = new double[ newWidth * newHeight * 2 ];

for( pos = 0; pos < count; pos += 2 ) {
    //( left.m_real * right.m_real ) - ( left.m_imaginary * right.m_imaginary );
    fOutputfast[ pos ] = fImagefast[pos] * kernel2dfast[ pos ] - fImagefast[ pos + 1 ] * kernel2dfast[ pos + 1 ];
    //( left.m_imaginary * right.m_real ) + ( left.m_real * right.m_imaginary );
    fOutputfast[ pos + 1 ] = fImagefast[ pos + 1 ] * kernel2dfast[ pos ] + fImagefast[ pos ] * kernel2dfast[ pos + 1 ];
}

sw.Stop();

MessageBox.Show( sw.ElapsedMilliseconds.ToString() );

The results, ~31msec!!! more than 3.5 times faster. Imagine how faster the whole code could be if just simple 1 dimension double arrays were used.

These are my first observations. New edits will come later.

EDIT2

In this edit we will optimize CustomFFt even more. Some observations. After some manipulation of the equations i come at this:

From the aboveCustomFft can be written as:

private void CustomFftFast( ref float[] kernel, int sizeX, int sizeY, double n, double f0, double theta, double a ) {
    int pos = 0;

    //pre-calculated values
    float rad1, b, cosrad1, sinrad1; //rad2, cosrad2, sinrad2;
    float Ax, Ay, Bx, By, Fo;
    float p1, p2, q1, q2, ka, kb, kw;

    rad1 = ( (float)Math.PI / 180.0f ) * ( 90.0f - (float)theta );
    //rad2 = rad1 + (float)Math.PI;
    cosrad1 = (float)Math.Cos( rad1 );
    sinrad1 = (float)Math.Sin( rad1 );
    //cosrad2 = (float)Math.Cos( rad2 );
    //sinrad2 = (float)Math.Sin( rad2 );
    b = ( (float)a / 5.0f ) / ( 2.0f * (float)n );

    Ax = cosrad1 / ( (float)sizeX * (float)a );
    Bx = sinrad1 / ( (float)sizeX * (float)b );

    Ay = sinrad1 / ( (float)sizeY * (float)a );
    By = cosrad1 / ( (float)sizeY * (float)b );

    Fo = (float)f0 / (float)a;

    kernel = new float[ sizeX * sizeY * 2 ];

    for( int x = 0; x < sizeX; x++ ) {
        for( int y = 0; y < sizeY; y++ ) {
            p1 = Ax * (float)x + Ay * (float)y + Fo;
            q1 = By * (float)y - Bx * (float)x;

            p2 = -Ax * (float)x - Ay * (float)y + Fo;
            q2 = -By * (float)y + Bx * (float)x;

            ka = p1 * p1 + q1 * q1;
            kb = p2 * p2 + q2 * q2;

            //Normal Pow
            ka = 1.0f / ( 1.0f + 0.414f * (float)Math.Pow( ka, n ) );
            kb = 1.0f / ( 1.0f + 0.414f * (float)Math.Pow( kb, n ) );

            /*//Faster Pow, loses in percision though. Ther is another one
              //with better percision but a little slower.
            unsafe
            {
                uint *c = (uint *)&ka;
                float y_ = (float)*c;
                y_ *= 1.0f / (float)( 1 << 23 );
                y_ -= 126.94269504f;

                uint a_ = (uint)( ( 1 << 23 ) * ( y_ * (float)n + 126.94269504f ) );
                float *b_ = (float *)&a_;
                ka = 1.0f / ( 1.0f + 0.414f * ( *b_ ) );

                c = (uint *)&kb;
                y_ = (float)*c;
                y_ *= 1.0f / (float)( 1 << 23 );
                y_ -= 126.94269504f;

                a_ = (uint)( ( 1 << 23 ) * ( y_ * (float)n + 126.94269504f ) );
                b_ = (float *)&a_;
                kb = 1.0f / ( 1.0f + 0.414f * ( *b_ ) );

            }*/

            kw = ka >= kb ? ka : kb;

            kernel[ pos ] = kw;;
            //kernel[ pos + 1 ] = 0.0f;

            pos += 2;
        }
    }
}

The bottleneck of CustomFFt is the Pow call. Times:

With Math.Pow(),    time = ~110msec
With faster Pow,    time = ~90msec, strange for such a small improvement
Without Pow at all, time = ~10msec! <--