图像处理 -- 图像模版匹配算法之NCC

图像模版匹配算法NCC的数学原理及实现

背景

NCC（Normalized Cross-Correlation，归一化互相关）是一种常用于图像处理中的模板匹配算法。其核心思想是计算模板图像与目标图像局部区域之间的相似度，并通过最大化NCC值找到最佳匹配位置。

数学原理

NCC的计算公式如下：

$$\text{NCC}(u,v)=\frac{{\sum }_{x,y}(T(x,y)-\bar{T})(I(x+u,y+v)-{\bar{I}}_{u,v})}{\sqrt{{\sum }_{x,y}(T(x,y)-\bar{T}{)}^2{\sum }_{x,y}(I(x+u,y+v)-{\bar{I}}_{u,v}{)}^2}}$$

其中：

• $T(x,y)$ 是模板图像的像素值。
• $I(x+u,y+v)$ 是目标图像在位置 $(u,v)$ 处的像素值。
• $\bar{T}$ 是模板图像的均值：
  $$\bar{T}=\frac{1}{N}\sum _{x,y}T(x,y)$$
• ${\bar{I}}_{u,v}$ 是目标图像在位置 $(u,v)$ 处对应窗口的均值：
  $${\bar{I}}_{u,v}=\frac{1}{N}\sum _{x,y}I(x+u,y+v)$$

1. 归一化处理

为了确保匹配计算不受亮度变化影响，我们需要对图像进行归一化处理：
$$T^{\prime }(x,y)=T(x,y)-\bar{T},I_{u,v}^{\prime }(x,y)=I(x+u,y+v)-{\bar{I}}_{u,v}$$

2. 计算标准差

标准差用于衡量图像局部区域的对比度：
$${\sigma }_T=\sqrt{\frac{1}{N-1}\sum _{x,y}(T(x,y)-\bar{T}{)}^2}$$
$${\sigma }_{I_{u,v}}=\sqrt{\frac{1}{N-1}\sum _{x,y}(I(x+u,y+v)-{\bar{I}}_{u,v}{)}^2}$$

3. 计算NCC匹配值

最终的NCC公式为：
$$\text{NCC}(u,v)=\frac{{\sum }_{x,y}(T(x,y)-\bar{T})(I(x+u,y+v)-{\bar{I}}_{u,v})}{{\sigma }_T\cdot {\sigma }_{I_{u,v}}}$$

代码实现

下面是一个优化后的 C 代码示例，进行NCC匹配，并修正了标准差计算和边界处理问题。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

#define IMAGE_WIDTH  640
#define IMAGE_HEIGHT 480
#define TEMPLATE_WIDTH 16
#define TEMPLATE_HEIGHT 16

// 读取灰度图像数据（假设数据为0-255的8位像素）
void load_image(const char *filename, uint8_t image[IMAGE_HEIGHT][IMAGE_WIDTH]) {
    FILE *file = fopen(filename, "rb");
    if (!file) {
        perror("Failed to open image file");
        exit(EXIT_FAILURE);
    }
    fread(image, sizeof(uint8_t), IMAGE_WIDTH * IMAGE_HEIGHT, file);
    fclose(file);
}

// 计算均值
double compute_mean(uint8_t *data, int width, int height) {
    double sum = 0.0;
    for (int i = 0; i < width * height; i++) {
        sum += data[i];
    }
    return sum / (width * height);
}

// 计算标准差
double compute_stddev(uint8_t *data, int width, int height, double mean) {
    double sum_sq = 0.0;
    for (int i = 0; i < width * height; i++) {
        sum_sq += (data[i] - mean) * (data[i] - mean);
    }
    return sqrt(sum_sq / (width * height - 1));
}

// 计算NCC
double compute_ncc(uint8_t template[TEMPLATE_HEIGHT][TEMPLATE_WIDTH], uint8_t *region, int width, int height) {
    // 计算模板均值和标准差
    double template_mean = compute_mean((uint8_t *)template, TEMPLATE_WIDTH, TEMPLATE_HEIGHT);
    double template_std = compute_stddev((uint8_t *)template, TEMPLATE_WIDTH, TEMPLATE_HEIGHT, template_mean);

    // 计算当前区域均值和标准差
    double region_mean = compute_mean(region, width, height);
    double region_std = compute_stddev(region, width, height, region_mean);

    // 计算NCC分子部分
    double numerator = 0.0;
    for (int y = 0; y < TEMPLATE_HEIGHT; y++) {
        for (int x = 0; x < TEMPLATE_WIDTH; x++) {
            numerator += (template[y][x] - template_mean) * (region[y * width + x] - region_mean);
        }
    }

    return numerator / (template_std * region_std);
}

// 进行NCC匹配
void ncc_matching(uint8_t image[IMAGE_HEIGHT][IMAGE_WIDTH], uint8_t template[TEMPLATE_HEIGHT][TEMPLATE_WIDTH], int *best_x, int *best_y, double *best_ncc) {
    *best_ncc = -1.0;

    for (int y = 0; y <= IMAGE_HEIGHT - TEMPLATE_HEIGHT; y++) {
        for (int x = 0; x <= IMAGE_WIDTH - TEMPLATE_WIDTH; x++) {
            // 取出子区域
            uint8_t region[TEMPLATE_HEIGHT * TEMPLATE_WIDTH];
            for (int ty = 0; ty < TEMPLATE_HEIGHT; ty++) {
                for (int tx = 0; tx < TEMPLATE_WIDTH; tx++) {
                    region[ty * TEMPLATE_WIDTH + tx] = image[y + ty][x + tx];
                }
            }

            // 计算NCC值
            double ncc = compute_ncc(template, region, TEMPLATE_WIDTH, TEMPLATE_HEIGHT);
            if (ncc > *best_ncc) {
                *best_ncc = ncc;
                *best_x = x;
                *best_y = y;
            }
        }
    }
}

int main() {
    uint8_t image[IMAGE_HEIGHT][IMAGE_WIDTH];
    uint8_t template[TEMPLATE_HEIGHT][TEMPLATE_WIDTH];

    // 加载图像数据
    load_image("image.raw", image);

    // 提取模板（假设固定位置）
    for (int y = 0; y < TEMPLATE_HEIGHT; y++) {
        for (int x = 0; x < TEMPLATE_WIDTH; x++) {
            template[y][x] = image[200 + y][300 + x];
        }
    }

    int best_x, best_y;
    double best_ncc;
    ncc_matching(image, template, &best_x, &best_y, &best_ncc);

    printf("Best match location: (%d, %d), NCC value: %lf\n", best_x, best_y, best_ncc);

    return 0;
}