Not sure if I missed some obvious CPU brakes

Question

I have three functions (out of a medium 2 digit number of functions) which use up 80% of CPU time, so I kind of wonder if there are points for optimization that I am missing.

Function 1, Extracting a bilinear interpolated point from a width * height float value image:

float ExtractBilinear(float* image, int w, int h, float x, float y)
{
    int x0 = (int)floor(x);
    int y0 = (int)floor(y);
    int x1 = x0 + 1;
    int y1 = y0 + 1;
    float c00, c01, c11, c10;
    c00 = c01 = c11 = c10 = 0.0f;

    if(x0 < -1 || x1 > w || y0 < -1 || y1 > h)
    {
        return 0.0f;
    }

    if(x0<0)
    {
        c00 = 0.0f;
        c01 = 0.0f;
    }
    else
    {
        if(y0<0)
        {
            c00 = 0.0f;
        }
        else
        {
            c00 = image[y0*w + x0];
        }
        if(y1 < h)
        {
            c01 = image[y1*w + x0];
        }
        else
        {
            c01 = 0.0f;
        }
    }

    if(x1 < w)
    {
        if(y0<0)
        {
            c10 = 0.0f;
        }
        else
        {
            c10 = image[y0*w + x1];
        }
        if(y1 < h)
        {
            c11 = image[y1*w + x1];
        }
        else
        { 
            c11 = 0.0f;
        }
    }
    else
    {
        c10 = 0.0f;
        c11 = 0.0f;
    }

    float c0 = c10 * (x - x0) + c00 * (x1 - x);
    float c1 = c11 * (x - x0) + c01 * (x1 - x);
    return c1 * (y - y0) + c0 * (y1 - y);
}

Function 2: Taking 2 float samples of data (one is image data, one is simulated data), and calculating the sum of squared differences error, with a scalar 'a' to minimize the error as much as possible. SampleX is the size of the input, compX the size of the output, and offX is deprecated, and always 0 (due to legacy code I am keeping it in there).

float PatternMatcher::GetSADFloatRel(float* sample, float* compared, int sampleX, int compX, int offX)
{
    if (sampleX != compX)
    {
        return 50000.0f;
    }
    float result = 0;

    float* pTemp1 = sample;
    float* pTemp2 = compared + offX;

    float w1 = 0.0f;
    float w2 = 0.0f;
    float w3 = 0.0f;

    for(int j = 0; j < sampleX; j ++)
    {
        w1 += pTemp1[j] * pTemp1[j];
        w2 += pTemp1[j] * pTemp2[j];
        w3 += pTemp2[j] * pTemp2[j];
    }               
    float a = w2 / w3;
    result = w3 * a * a - 2 * w2 * a + w1;
    return result / sampleX;
} 

Function 3: a 2 dimensional convolution over the image with rotations being applied. This is working multi-threaded, and is, actually, the most CPU intensive function in our whole algorithm.

Input: image is the original float image.

convolution_image holds the highest result for that particular thread (I start #CPU Cores - 2 threads, and each thread has its own copy of convolution_image, r_map_image and orientation_map_image).

r_map_image saves the currently found radius at that position (with highest convolution value). orientation_map_image saves the angle, and goes from 0..180.

middleReflex is a pre-calculated map of where the image has to be invalidated due to unavoidable reflexes happening in the optical system.

kernel1DAll holds the pre-calculated sample of what an area has to look like to be used later in the algorithm.

startMin, startMax hold the angles in which the calculation for the current thread should run (for example, with 3 threads, 0..60, 61..120, 121..180)

outer is a Circle, containing a position (2D coordinates), and a radius, and every point outside that circle is completely unimportant and should be ignored.

void rotateImageConvolution(float* image, float* convolution_image, unsigned char* r_map_image, unsigned char* orientation_map_image, unsigned char* middleReflex, float* kernel1DAll, int startMin, int startMax, Circle outer)
{
    const int imgsize = Width*Height;
    printf("Started thread from %d to %d\n", startMin, startMax);

    float half_w = Width/2.0f;
    float half_h = Height/2.0f;
    const int radius =  int(sqrt(half_w*half_w + half_h*half_h) + 0.5f);
    const int RotateW = 2*radius + 1;
    const int s = RotateW * RotateW; //* nAngles;

    float* rotate_image = (float*)malloc(s*sizeof(float));
    float* summingup_image = (float*)malloc(s*sizeof(float));

    for(int a = startMin ; a<startMax ; a++)
    {
        float theta = (float)(PI*a/nAngles);
        const float costheta = cos(theta);
        const float sintheta = sin(theta);
        Cos[a] = costheta;
        Sin[a] = sintheta;

        for(int j = 0; j<RotateW; j++)
        {
            int y = j - radius;
            float ysina = y*sintheta;
            float ycosa = y*costheta;
            for(int i = 0; i<RotateW; i++)
            {
                int x = i - radius;

                float xf = x*costheta - ysina + Width/2;
                float yf = x*sintheta + ycosa + Height/2;

                rotate_image[j*RotateW + i] = ExtractBilinear(image, Width, Height, xf, yf);
            }
        }

        for(int j = 0; j<RotateW; j++)
        {
            int yoff = j*RotateW;
            float subSum = 0.0f;
            int num = 0;
            for(int kx=0; kx<=half_Integ_L && kx < RotateW; kx++)
            {               
                subSum += rotate_image[yoff + kx];
                num++;
            }

            summingup_image[yoff] = subSum/num;

            for(int i=1; i<RotateW; i++)
            {
                float sum = subSum;
                int istart = i - half_Integ_L;
                int iend = i + half_Integ_L;

                if(istart>=0)
                {
                    num--;
                    subSum -= rotate_image[yoff + istart];
                }
                if(iend<RotateW)
                {
                    num++;
                    float endv = rotate_image[yoff + iend];
                    sum += endv;
                    subSum += endv;
                }

                summingup_image[yoff+i] = sum/num;
            }
        }

        for(int i=0;i<imgsize;i++)
        {
            if(middleReflex[i] == 0xFF)
            {
                convolution_image[i] = 0.0f;
                r_map_image[i] = 10;
                orientation_map_image[i] = 0;
                continue;
            }

            int x = i%Width - Width/2;
            int y = i/Width - Height/2;

            float d = sqrt((i/Width - outer.y) * (i/Width - outer.y) + (i%Width - outer.x) * (i%Width - outer.x));
            if(d >= outer.r-16)
            {
                convolution_image[i] = 0.0f;
                r_map_image[i] = 10;
                orientation_map_image[i] = 0;
                continue;
            }

            int xx = (int)floor(x*costheta + y*sintheta + radius + 0.5f);
            int yy = (int)floor(-x*sintheta + y*costheta + radius + 0.5f);

            float I[KWIDTH];
            for(int ky=-half_MaxK; ky<=half_MaxK; ky++)
            {
                int dy = yy + ky;
                int ii = ky+half_MaxK;

                if(dy >= 0 && xx >=0 && dy < RotateW && xx < RotateW)
                {
                    I[ii] = summingup_image[dy*RotateW + xx];
                }
                else I[ii] = 0.0f;
            }

            for(int ik = NKernel-1; ik>=0; ik--)
            {
                float *pKernel = kernel1DAll + ik*KWIDTH;

                float convrst = 0.0f;

                int r0 = ik + half_MinK;
                for(int x0 = -r0; x0<=r0; x0++ )
                {
                    int ii = x0 + half_MaxK;
                    convrst += (I[ii]*pKernel[KWIDTH-ii-1]);                
                }

                if(convolution_image[i] <= convrst)
                {
                    convolution_image[i] = convrst;
                    r_map_image[i] = (unsigned char)r0;
                    orientation_map_image[i] = (unsigned char)a;
                }
            }
        }
    }

    free(summingup_image);
    free(rotate_image);
}

Am I doing something horribly wrong with this code with regards to performance? The code does what it is supposed to do, but unless running on a high-end system, it can take really long to finish.

Answer 1

Pre-computing things is useful when it is possible. As some one has already said in your case may be advisable to precomunte things as Width/2. But, if Width is defined through a "Define" your compiler shall be able to do this optimization by itself.

But, I believe that the problem in your code is mainly another one. You have a first loop between startMin and StartMax and than another nested loop from o to imageSize (which i expect being a large number). And this ends up to be the most time expensive part of your code.

There you have other two loops such as: for(int ik = NKernel-1; ik>=0; ik--) or for(int ky=-half_MaxK; ky<=half_MaxK; ky++)

This you may think that is an "optimized" loop, but actually it isn't, because the compiler is (generally) designed to be really good at optimizing loops from 0 to N.

Thus you can have significant improvement. Write your two inner loops in such way and then just compute your index (ik and ky) by subtracting an offset.

Answer 2

How often is the third function called?

The allocation of the images could be a bottleneck. You could try to reuse the memory for the images by passing them to the function instead of reallocating and freeing them every time.

Also things like Width/2.0f, Height/2.0f and PI*a/nAngles should not be evaluated every time they are needed if they are not changed.