Correction for barrel and pincushion distortion

Question

I have a c++ dll I am calling (PInvoking) from C# which corrects barrel and pincushion distortion. Basically I call a method which generates a set of look up arrays and I adjust the pixels of the input image accordingly. However this is causing some serious performance issues with 50% of my applications CPU usage being in this method. Could I improve the speed and or CPU usage any more? I am not a native C++ coder and I have done my best to speed this up but I cant seem to squeeze anything else out of the method.

//Global arrays for the image distortion correction
unsigned short MapBuffer1_x[5067304], MapBuffer2_x[5067304];
unsigned long MapBuffer1_y[5067304], MapBuffer2_y[5067304];
unsigned char MapBuffer_ratio_x[5067304], MapBuffer_ratio_y[5067304];

////////////////////////////////////////////////////////////////////////////////
// Corrects image distortion, GenerateLensMap must be called first
//
// srcBuffer = pointer to the source buffer
// dstBuffer = pointer to destination buffer
// nWidth = width of image
// nHeight = height of image
//
////////////////////////////////////////////////////////////////////////////////
BOOL __stdcall CorrectDistortion(unsigned char *srcBuffer, unsigned char *dstBuffer, DWORD nWidth, DWORD nHeight)
{
    unsigned long x1, y1, r1, r2, r, end, dstPix, srcPix, srcPix_x, srcPix_y, pix, v1, v2;
    //unsigned long srcPixB;
    signed long x2, y2, width3;

    width3 = nWidth * 3;
    end = nHeight * width3;

    pix = 0;

    //scan dstPix through entire destination image
    for (dstPix = 0; dstPix < end; dstPix += 3)
    {
        //Required source pixel must be interpolated between three adjacent pixels
        //Get 1st (main) required source pixel from mapping buffers MapBuffer1
        x1 = MapBuffer1_x[pix];
        y1 = MapBuffer1_y[pix];

        //Get 2nd required source pixel in x direction from mapping buffer MapBuffer2_x. This will share same y value as 1st pixel
        x2 = MapBuffer2_x[pix];

        //Get 3rd required source pixel in y direction from mapping buffer MapBuffer2_y. This will share same x value as 1st pixel
        y2 = MapBuffer2_y[pix];

        v1 = MapBuffer_ratio_x[pix];
        v2 = MapBuffer_ratio_y[pix];

        //The mix ratio between the first pixel and these two other pixels is next obtained
        //(Ratio value in buffer is 128 times too big. This avoids using slow floating point numbers.)
        r = v1 + v2;//proportion of the first pixel in both x and y direction to use

        r1 = 128 - v1;//This is effectively 1 - ratio being the proportion of the second pixel (x direction) to use
        r2 = 128 - v2;//This is effectively 1 - ratio being the proportion of the third pixel (y direction) to use

        srcPix = y1 + x1;//main source pixel
        srcPix_x = y1 + x2;//second source pixel (x direction)
        srcPix_y = y2 + x1;//third source pixel (y direction)

                           //destination pixel is average of (first pixel interpolated with (second) pixel offset in x direction) AND
                           //(first pixel interpolated with (third) pixel offest in y direction)

        dstBuffer[dstPix] = (r* srcBuffer[srcPix] + r1 * srcBuffer[srcPix_x] + r2 * srcBuffer[srcPix_y]) >> 8;   //B
        dstBuffer[dstPix + 1] = (r*srcBuffer[srcPix + 1] + r1*srcBuffer[srcPix_x + 1] + r2*srcBuffer[srcPix_y + 1]) >> 8; //G
        dstBuffer[dstPix + 2] = (r*srcBuffer[srcPix + 2] + r1*srcBuffer[srcPix_x + 2] + r2*srcBuffer[srcPix_y + 2]) >> 8; //R

        pix++;
    }
    return 0;
}

This is the method that generates the lookup tables (Map Buffers)

////////////////////////////////////////////////////////////////////////////////
//
// nWidth = width of image
// nHeight = height of image
// distortion = amount of distortion
//
LONG __stdcall GenerateLensMap(DWORD nWidth, DWORD nHeight, LONG Distortion)
//Generates the re-mapping look-up tables for correcting lens distortion.
//This routine is called whenever the zoom lens is zoomed.
//Distortion is reciprocal of required value
//Corrected image pixel is average of (first pixel interpolated with (second) pixel offest in x direction) AND
//(first pixel interpolated with (third) pixel offest in y direction)
{
    long       py, px, pix, x1, y1;
    short width3;
    unsigned long x, y, halfwidth, halfheight;
    float       k2, m, adj, x_sub, y_sub, remainder;

    width3 = short(nWidth * 3);
    halfwidth = nWidth / 2;
    halfheight = nHeight / 2;
    k2 = (float)1 / Distortion;

    adj = (float)(1 / (1 + fabs(k2) * ((halfwidth * halfwidth) + (halfheight*halfheight)))); //reciprocal of maximum correction based on corner pixel

    for (y = 0, pix = 0; y < nHeight; y++) {
        py = y - halfheight; //will scan between -halfheight to +halfheight
        for (x = 0; x < nWidth; x++, pix++) {
            px = x - halfwidth; //will scan between -halfwidth to +halfwidth
            m = (1 + k2 * ((py * py) + (px * px))) * adj; //correction to be applied to current pixel based on distance from centre

                                                          //calculate corrected pixel x position to sub pixel accuracy
            x_sub = (px * m) + halfwidth;
            if (x_sub < 0) x_sub = 0;//prevent overflowing buffer
            if (x_sub > (signed)nWidth) x_sub = (float)nWidth;
            remainder = x_sub - (int)x_sub;
            MapBuffer1_x[pix] = (int)x_sub * 3; //*3 required to step thru colour channels
            if (remainder < 0.5)//if subpixel to left of centre, interpolate between left pixel else interpolate between right pixel
                x1 = (int)x_sub - 1;
            else
                x1 = (int)x_sub + 1;
            if (x1 < 0) x1 = 0;//prevent overflowing buffer
            if (x1 > (signed)nWidth) x1 = nWidth;
            MapBuffer2_x[pix] = (int)x1 * 3;  //*3 required to step thru colour channels

                                              //Calculate mix ratio as 2*(0.5 - |remainder - 0.5|), i.e. ranges from 0 to 1 to 0 as passing through pixel
            remainder = remainder - 0.5F;
            if (remainder < 0) remainder = -remainder;//make positive
            MapBuffer_ratio_x[pix] = (unsigned char)(128 * (0.5F - remainder));//scale to 128

                                                                               //calculate corrected pixel y position to sub pixel accuracy
            y_sub = (py * m) + halfheight;
            if (y_sub < 0) y_sub = 0;//prevent overflowing buffer
            if (y_sub > (signed)nHeight) y_sub = (float)nHeight;
            remainder = y_sub - (int)y_sub;
            MapBuffer1_y[pix] = (int)y_sub * width3; //* width3 required to step thru lines in VSCDistortionCorrection
            if (remainder < 0.5)//if subpixel to above centre, interpolate between upper pixel else interpolate between lower pixel
                y1 = (int)y_sub - 1;
            else
                y1 = (int)y_sub + 1;
            if (y1 < 0) y1 = 0;//prevent overflowing buffer
            if (y1 > (signed)nHeight) y1 = nHeight;
            MapBuffer2_y[pix] = (int)y1 * width3; //* width3 required to step thru lines in VSCDistortionCorrection;

                                                  //Calculate mix ratio as 2*(0.5 - |remainder - 0.5|), i.e. ranges from 0 to 1 to 0 as passing through pixel
            remainder = remainder - 0.5F;
            if (remainder < 0) remainder = -remainder;//make positive
            MapBuffer_ratio_y[pix] = (unsigned char)(128 * (0.5F - remainder));//scale to 128
        }
    }

    return Distortion;
}

Answer 1

Globals

There's no reason to use global variables here. You can pass the tables into both the GenerateLensMap() and CorrectDistortion() functions. If you need to wrap them in a struct or class so you're only passing a single parameter around, that's fine. Just don't use global variables. You never know when some other function is going to change them.

Magic Numbers

Why are your buffers all 5,067,304 elements long? What is the significance of that number? Whatever it is, you should name it. I recommend using a named constant like:

const size_t kMaxBufferSize = 5067304;

Even better if you can name it something more meaningful like kMaxImageSizeInPixels or kMaxPixelsPerImage or whatever it represents. (Also, why make it a fixed value? Why not allocate it with the exact size you need?)

Structures

You have data that has a definite structure, yet you've decided to ignore that structure and treat everything like a big 1 dimensional array of unsigned characters. It makes your code much harder to read and understand. Define a pixel type like this:

typedef struct BGRPixel {
    unsigned char blue;
    unsigned char green;
    unsigned char red;
} BGRPixel;

Then instead of writing this:

    dstBuffer[dstPix] = (r* srcBuffer[srcPix] + r1 * srcBuffer[srcPix_x] + r2 * srcBuffer[srcPix_y]) >> 8;   //B
    dstBuffer[dstPix + 1] = (r*srcBuffer[srcPix + 1] + r1*srcBuffer[srcPix_x + 1] + r2*srcBuffer[srcPix_y + 1]) >> 8; //G
    dstBuffer[dstPix + 2] = (r*srcBuffer[srcPix + 2] + r1*srcBuffer[srcPix_x + 2] + r2*srcBuffer[srcPix_y + 2]) >> 8; //R

you can write this:

    dstBuffer[dstPix].blue = (r* srcBuffer[srcPix].blue + r1 * srcBuffer[srcPix_x].blue + r2 * srcBuffer[srcPix_y].blue) >> 8;
    dstBuffer[dstPix].green = (r*srcBuffer[srcPix].green + r1*srcBuffer[srcPix_x].green + r2*srcBuffer[srcPix_y].green) >> 8;
    dstBuffer[dstPix].red = (r*srcBuffer[srcPix].red + r1*srcBuffer[srcPix_x].red + r2*srcBuffer[srcPix_y].red) >> 8;

Now you don't have to explain that the first line handles blue, the second green, and the third red.

Furthermore, you can overload operator*() to multiply a scalar by a pixel and overload operator+() to add 2 pixels together. And you could define the operator>>() to do the final shift. That would leave you with code like this:

dstBuffer[dstPix] = (r * srcBuffer[srcPix] + r1 * srcBuffer[srcPix_x + r2 * srcBuffer[srcPix_y]) >> 8;

Much easier to read and understand!

Making the struct also removes the need to constantly add and multiply various values by 3. The for loop in CorrectDistortion() would just increment by 1 each time instead of by 3, and the various multiplies by 3 in GenerateLensMap() would go away.

Likewise, as mentioned in the comments by @Toby Speight, you can also make the various global arrays into a single array of structs something like this:

typedef struct LensDistortionMapEntry {
    unsigned short x1_offset;
    unsigned long y1_offset;
    unsigned short x2_offset;
    unsigned long y2_offset;
    unsigned char x_ratio;
    unsigned char y_ratio;
} LensDistortionMapEntry;

Interpolation

Your method of interpolation is a bit odd. Usually one would use 4 pixels to interpolate between using something like bilinear interpolation. This would do 2 linear interpolations in x (one on the row just below the actual location and one on the row just above), and then do another linear interpolation between those 2 results. I've never seen the type of interpolation you're doing before.

Performance

Here's where it gets fun. Overall what you have is pretty good and straightforward. There are a number of ways you can speed this up. The easiest would be to perform the work on multiple threads at the same time. For example, you could figure out how many cores the machine you're running on has and spin up that many threads. For each image send (height / numThreads) rows of pixels to each thread. You're likely to get something approaching a linear speedup out of this since it's so easily parallelizable.

It addition, if you instead did tiles that always fit within the CPU's cache lines you may be able to eke out even more performance from each of those threads.

Another way to do it would be to send it off to the GPU to do the work. I often do image processing on the GPU using OpenGL. There's a bit of a learning curve to get it up and running, but once you've done it once, it's pretty easy to make new filters that work on the GPU. The speedup for something like this with a GPU also scales well with the number of cores, and modern GPUs have hundreds to thousands of cores, so it's a pretty huge speedup. (Far better than anything you'll get on the CPU alone.) You'd want to generate the lens map as a texture, then send that as a texture along with the input image to image processing piece. Both can be done on the GPU.

Errors

I should also mention that it looks like you might have an error in GenerateLensMap(). In it, you have this line:

if (x_sub > (signed)nWidth) x_sub = (float)nWidth;

However, you don't want to ever access the nWidthth pixel on a line because you'll go past the end of the line. I think you want nWidth - 1 for both the conditional and the assignment. (Or the conditional could be >= instead of >.) Likewise with the checks against the height of the image.