Optimizing C++ code for performance

Question

Can you think of some way to optimize this piece of code? It's meant to execute in an ARMv7 processor (Iphone 3GS):

4.0%  inline float BoxIntegral(IplImage *img, int row, int col, int rows, int cols) 
      {
0.7%    float *data = (float *) img->imageData;
1.4%    int step = img->widthStep/sizeof(float);

        // The subtraction by one for row/col is because row/col is inclusive.
1.1%    int r1 = std::min(row,          img->height) - 1;
1.0%    int c1 = std::min(col,          img->width)  - 1;
2.7%    int r2 = std::min(row + rows,   img->height) - 1;
3.7%    int c2 = std::min(col + cols,   img->width)  - 1;

        float A(0.0f), B(0.0f), C(0.0f), D(0.0f);
8.5%    if (r1 >= 0 && c1 >= 0) A = data[r1 * step + c1];
11.7%   if (r1 >= 0 && c2 >= 0) B = data[r1 * step + c2];
7.6%    if (r2 >= 0 && c1 >= 0) C = data[r2 * step + c1];
9.2%    if (r2 >= 0 && c2 >= 0) D = data[r2 * step + c2];

21.9%   return std::max(0.f, A - B - C + D);
3.8%  }

All this code is taken from the OpenSURF library. Here's the context of the function (some people were asking for the context):

//! Calculate DoH responses for supplied layer
void FastHessian::buildResponseLayer(ResponseLayer *rl)
{
  float *responses = rl->responses;         // response storage
  unsigned char *laplacian = rl->laplacian; // laplacian sign storage
  int step = rl->step;                      // step size for this filter
  int b = (rl->filter - 1) * 0.5 + 1;         // border for this filter
  int l = rl->filter / 3;                   // lobe for this filter (filter size / 3)
  int w = rl->filter;                       // filter size
  float inverse_area = 1.f/(w*w);           // normalisation factor
  float Dxx, Dyy, Dxy;

  for(int r, c, ar = 0, index = 0; ar < rl->height; ++ar) 
  {
    for(int ac = 0; ac < rl->width; ++ac, index++) 
    {
      // get the image coordinates
      r = ar * step;
      c = ac * step; 

      // Compute response components
      Dxx = BoxIntegral(img, r - l + 1, c - b, 2*l - 1, w)
          - BoxIntegral(img, r - l + 1, c - l * 0.5, 2*l - 1, l)*3;
      Dyy = BoxIntegral(img, r - b, c - l + 1, w, 2*l - 1)
          - BoxIntegral(img, r - l * 0.5, c - l + 1, l, 2*l - 1)*3;
      Dxy = + BoxIntegral(img, r - l, c + 1, l, l)
            + BoxIntegral(img, r + 1, c - l, l, l)
            - BoxIntegral(img, r - l, c - l, l, l)
            - BoxIntegral(img, r + 1, c + 1, l, l);

      // Normalise the filter responses with respect to their size
      Dxx *= inverse_area;
      Dyy *= inverse_area;
      Dxy *= inverse_area;

      // Get the determinant of hessian response & laplacian sign
      responses[index] = (Dxx * Dyy - 0.81f * Dxy * Dxy);
      laplacian[index] = (Dxx + Dyy >= 0 ? 1 : 0);

#ifdef RL_DEBUG
      // create list of the image coords for each response
      rl->coords.push_back(std::make_pair<int,int>(r,c));
#endif
    }
  }
}

Some questions:
Is it a good idea that the function is inline? Would using inline assembly provide a significant speedup?

Answer 1

Specialize for the edges so that you don't need to check for them in every row and column. I assume that this call is in a nested loop and is called a lot. This function would become:

inline float BoxIntegralNonEdge(IplImage *img, int row, int col, int rows, int cols) 
{
  float *data = (float *) img->imageData;
  int step = img->widthStep/sizeof(float);

  // The subtraction by one for row/col is because row/col is inclusive.
  int r1 = row - 1;
  int c1 = col - 1;
  int r2 = row + rows - 1;
  int c2 = col + cols - 1;

  float A(data[r1 * step + c1]), B(data[r1 * step + c2]), C(data[r2 * step + c1]), D(data[r2 * step + c2]);

  return std::max(0.f, A - B - C + D);
}

You get rid of a conditional and branch for each min and two conditionals and a branch for each if. You can only call this function if you already meet the conditions -- check that in the caller for the whole row once instead of each pixel.

I wrote up some tips for optimizing image processing when you have to do work on each pixel:

http://www.atalasoft.com/cs/blogs/loufranco/archive/2006/04/28/9985.aspx

Other things from the blog:

1. You are recalculating a position in the image data with 2 multiplies (indexing is multiplication) -- you should be incrementing a pointer.

2. Instead of passing in img, row, row, col and cols, pass in pointers to the exact pixels to process -- which you get from incrementing pointers, not indexing.

3. If you don't do the above, step is the same for all pixels, calculate it in the caller and pass it in. If you do 1 and 2, you won't need step at all.