feature: correlation coefficient template matching method

[sby1ce]

Hello.

imageproc currently has 2 template matching methods, sum of squared errors (SQDIFF) and cross correlation (CCORR), as opposed to OpenCV's 3 methods, which also has correlation coefficient (CCOEFF).

Saying it outright, this feature is mostly pointless.
From my own work, it's obvious that without using FFT and the work being done purely on the CPU this method is way too slow to be used for any practical purpose, especially when OpenCV already exists.
However I would like to share the math, because frankly speaking, OpenCV documentation and source code are unreadable.
I would also like to use this as a way to ask if anyone knows how is OpenCV able to make its OpenCL so fast despite seemingly suboptimal code.

So here's the math.

$$\begin{align*} R(x, y) = \sum_{x',y'} T'(x',y') \cdot I'(x+x',y+y') \\\ = \sum_{x',y'} \big((T(x',y') - \mu_T) \cdot (I(x+x',y+y') - \mu_{I_w})\big) \\\ = \sum_{x',y'} \big( (T(x',y') \cdot I(x+x',y+y') - (T(x',y') \cdot \mu_{I_w}) - (I(x+x',y+y') \cdot \mu_T) + (\mu_T \cdot \mu_{I_w}) \big) \\\ = \sum_{x',y'} \big( (T(x',y') \cdot I(x+x',y+y') - (T(x',y') \cdot \mu_{I_w}) - (I(x+x',y+y') \cdot \mu_T) \big) + \mu_T \mu_{I_w} \text{Area} \\\ = \sum_{x',y'} \big(T(x',y') \cdot I(x+x',y+y')\big) - \mu_{I_w} \sum_{t \in T} t - \mu_T \sum_{x',y'} I(x+x',y+y') + \mu_{I_w} \sum_{t \in T} t \\\ = \sum_{x',y'} \big(T(x',y') \cdot I(x+x',y+y')\big) - \mu_T \sum_{x',y'} I(x+x',y+y') \\ \text{where } T' \text{ and } I' \text{ are mean shifted} \\\ T'(x',y') = T(x', y') - 1 / (w \cdot h) \cdot \sum_{x'', y''} T(x'', y'') \\\ I'(x+x',y+y') = I(x+x',y+y') - 1 / (w \cdot h) \cdot \sum_{x'',y''} I(x+x'',y+y'') \\\ w \text{ - width of the template, } h \text{ - height of the template} \end{align*}$$

As mean of the template is constant and sum over image window can be computed with an integral image we're left with the pairwise multiplication of image window over template for every pixel in the output image.
This is obviously not very good and I don't see a way around it without FFT convolution because even on the GPU with GLSL, it's way too slow.

Here's what the Rust code would look like

/// R(x,y)
fn calculate_ccoeff_at(
    x: u32,
    y: u32,
    image: &ImageBuffer<Luma<u8>, Vec<u8>>,
    template: &ImageBuffer<Luma<u8>, Vec<u8>>,
    template_mean: f32,
    summed_area_table: &ImageBuffer<Luma<f32>, Vec<f32>>,
) -> f32 {
    let image_window_sum: f32 = calculate_window(x, y, template.dimensions(), summed_area_table);

    let product: f32 = template
        .enumerate_pixels()
        .map(|(x_prime, y_prime, template_pixel)| -> f32 {
            let template_pixel: f32 = f32::from(template_pixel[0]);
            let image_pixel: f32 = f32::from(image.get_pixel(x + x_prime, y + y_prime)[0]);

            template_pixel * image_pixel
        })
        .sum();

    let image_window_product: f32 = template_mean * image_window_sum;

    product - image_window_product
}