The CUDA Handbook ▸ Chapter 15. Image Processing: Normalized Correlation

15.1 Overview

Two 2D images, the image and the template, are compared by computing a correlation coefficient as follows:

\[\gamma(s,t) = \frac{\sum_{x}^{}{\sum_{y}^{}{\left\lbrack I(x,y) - \overline{I}(x,y) \right\rbrack\left\lbrack T(x - s,y - t) - \overline{T} \right\rbrack}}}{\sqrt{\sum_{x}^{}{\sum_{y}^{}{\left\lbrack I(x,y) - \overline{I}(x,y) \right\rbrack^{2}\sum_{x}^{}{\sum_{y}^{}\left\lbrack T(x - s,y - t) - \overline{T} \right\rbrack^{2}}}}}}\]

Where I and T are the image and template, respectively; \(\overline{T}\) is the average value of the template, while \(\overline{I}\) is the average value of the image pixels corresponding to the template.

The value of this coefficient falls into the range [-1.0, 1.0]; a value of 1.0 corresponds to a perfect match.

An optimized implementation of normalized correlation factors out the statistics that may be precomputed, and computes sums instead of averages to avoid a separate pass over the input data. If N pixels are being compared, replacing \(\overline{I}\) by \(\frac{\sum_{x}^{}{\sum_{y}^{}\left\lbrack I(x,y) \right\rbrack}}{N}\) and multiplying numerator and denominator by N yields a coefficient that can be expressed entirely in terms of sums. Rewriting without the coordinate notation:

\[\frac{N\sum_{}^{}{IT} - \sum_{}^{}I\sum_{}^{}T}{\sqrt{\left( N\sum_{}^{}{I^{2} - \left( \sum_{}^{}I \right)^{2}} \right)\left( N\sum_{}^{}{T^{2} - \left( \sum_{}^{}T \right)^{2}} \right)}}\]

Assuming the template will be the same for many correlation computations, the statistics on the template \(\sum_{}^{}T\) and \(\sum_{}^{}T^{2}\) can be precomputed, as can the subexpression \(\left( N\sum_{}^{}{T^{2} - \left( \sum_{}^{}T \right)^{2}} \right)\) in the denominator.

Translating this notation to C variable names:

Then a normalized correlation value may be computed using this function:

float
CorrelationValue( float SumI, float SumISq,
  float SumT, float SumTSq, float SumIT, float N )
{
    float Numerator = N*SumIT - SumI*SumT;
    float Denominator = (N*SumISq - SumI*SumI)*
                        (N*SumTSq - SumT*SumT);
    return Numerator / sqrtf(Denominator);
}

In practical applications for this algorithm, the template is kept fixed across many invocations, matching against different offsets into an image; then it makes sense to precompute the template statistics and the denominator subexpression:

float fDenomExp = N*SumSqT - SumT*SumT;

In practice, it’s best to use double precision to compute fDenomExp:

float fDenomExp = (float) ((double) N*SumSqT – (double) SumT*SumT);

and, on GPUs at least, it is faster to multiply by the reciprocal square root than to divide by the square root. The resulting CorrelationValue() function is as follows:

float
CorrelationValue( float SumI, float SumISq, float SumIT,
                  float N, float fDenomExp )
{
    float Numerator = cPixels*SumIT - SumI*SumT;
    float Denominator = (cPixels*SumISq - SumI*SumI)*fDenomExp;
    return Numerator * rsqrtf(Denominator);
}

Hence, an optimized implementation of this algorithm need only compute three sums over the pixels to compute a given correlation coefficient: \(\sum_{}^{}I\), \(\sum_{}^{}I^{2}\) and \(\sum_{}^{}{IT}\). Since the SMs include hardware support for integer multiply-add, NVIDIA GPUs are able to perform this computation extremely fast.

Of those three sums, \(\sum_{}^{}{IT}\) is the most expensive. \(\sum_{}^{}I\) and \(\sum_{}^{}I^{2}\) – the first and second moment of the pixels under the window, their sum and the sum of their squares – depend only on the image and the window’s position, so they can be tabulated once and read back in constant time (Section 15.8); but \(\sum_{}^{}{IT}\) combines the image with the template and changes with every shift, so the number and character of these sums determines the choice of algorithm. This chapter computes \(\sum_{}^{}{IT}\) three ways, each matched to a different type of search:

• a spatial kernel that slides one template across the image, the direct approach optimized through Sections 15.2–15.5 – best for a single, small template;

• an INT8 Tensor-Core matrix multiply that matches a whole bank of templates at once (Section 15.6) – best when many templates are searched together;

• an FFT that evaluates every shift in one transform, at a cost independent of the template size (Section 15.7) – best for a single large template searched exhaustively.

All three share one denominator, built once as a summed-area table (Section 15.8).

CUDA offers a number of paths that could be used to deliver the data to the hardware:

• global memory or texture memory for the image, the template, or both;

• constant memory for the template and possibly other template-specific parameters (up to 64K); and

• shared memory to hold image and/or template values for reuse.

This chapter will assume the pixels are 8-bit grayscale. The hardware works very well on images with higher precision, but if anything, that simplifies the problem by making it easier to efficiently address global and shared memory.

All of the CUDA implementations in this chapter use texture to read the image that is being compared with the template. There are several reasons for this:

• the texture units deal with boundary conditions gracefully and efficiently,

• the texture cache aggregates external bandwidth on reuse, which will occur as nearby correlation values are computed, and

• the 2D locality of the texture cache is a good fit with the access patterns exhibited by correlation search algorithms.

We’ll explore the tradeoffs of using texture versus constant memory for the template.