// Computes the integral image of image and the integral image of
// the squares of the elements in image.  Assumes image to be a 32-bit floating
// point.  Resultant intimage and sqintimage are also 32-bit floating point and
// are created automatically.
void integralImage(IplImage *image, IplImage *&intimage, IplImage *&sqintimage)
{
    int W = image->width, H = image->height;

    // Create resulting images
    intimage = cvCreateImage( cvSize(W, H), IPL_DEPTH_32F, 1 );
    sqintimage = cvCreateImage ( cvSize(W, H), IPL_DEPTH_32F, 1);

    float *i_m = (float *) image->imageData;
    float *ii_m = (float *) intimage->imageData;
    float *sii_m = (float *) sqintimage->imageData;
    
    ii_m[0] = i_m[0];
    sii_m[0] = i_m[0] * i_m[0];
    
    // Create the first row of the integral image
    for (int x = 1; x < W; x++)
    {
        ii_m[x] = ii_m[x-1] + i_m[x];
        sii_m[x] = sii_m[x-1] + i_m[x]*i_m[x];
    }
    
    // Compute each other row/column
    for (int y = 1, Y = W, YY=0; y < H; y++, Y+=W, YY+=W)
    {
        // Keep track of the row sum
        float r = 0, rs = 0;
        
        for (int x = 0; x < W; x++)
        {
            r += i_m[Y + x];
            rs += i_m[Y + x]*i_m[Y + x];
            ii_m[Y + x] = ii_m[YY + x] + r;
            sii_m[Y + x] = sii_m[YY + x] + rs;
        }
    }
}
//---------------------------------------------------------------------------