Банк фильтров занимает слишком много времени для создания ядер

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;
    }
}

Я создал 40-ядерный банк фильтров , где каждое ядро ​​будет иметь одинаковый размер, но разные параметры.

Этот банк фильтров занимает около 55000 миллисекунд для генерируемых ядер.

Как мне сократить это время?

Примечание: следующая строка занимает больше всего времени:

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

.

Пример исходного кода для создания банка ядра:

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());
    }
}

^{Примечание: это не исполняемый код.Этот код просто для того, чтобы дать представление о моем методе генерации ядра.}

Answer 1

В качестве отправной точки попробуйте это:

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 ) );
}

Я собираюсь еще немного поработать над этим.

РЕДАКТИРОВАТЬ 1

Мое первое наблюдение заключается в том, что вы используете слишком много классов, а в этих классах слишком много методов. Проблема с ООП состоит в том, что иногда код делится на слишком много классов с очень крошечными методами, что, в свою очередь, может затруднить выполнение кода (прохождение всех классов и т. Д.) И излишне медленное выполнение простых задач. Позвольте использовать несколько чисел:

Простое:

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

занимает около 65 мсек ! Он просто заменяет действительную часть cImage значением dImage, а мнимое устанавливается на zero. Это не должно занимать так много, но, как я сказал, ООП и классы. Давайте сделаем простую переконфигурацию и вместо Array2d<Complex> используем простой двойной массив:

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

Вот и все. Он хранит real-imaginary частей линейно, например

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]
....

Позволяет переписать код и проверить его:

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() );

Время, необходимое ~ 10 мсек . Более чем в 6 раз быстрее !!!!

Давайте сделаем еще один. Давайте попробуем FFT convolution, который займет ~ 105 мсек

// 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 ];
    }
}

Перепишите его с простыми одномерными двойными массивами:

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() );

Результаты ~ 31 мсек. более чем в 3,5 раза быстрее. Представьте, насколько быстрее мог бы быть весь код, если бы использовались только простые одномерные двойные массивы.

Это мои первые наблюдения. Новые правки появятся позже.

EDIT2

В этом редакторе мы оптимизируем CustomFFt еще больше. Некоторые наблюдения. После некоторой манипуляции с уравнениями я прихожу к этому:

Из вышеприведенного CustomFft можно записать как:

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;
        }
    }
}

Узким местом CustomFFt является Pow вызов. Times:

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