C++实现图像局部区域的亚像素极值定位与曲面拟合

引言

我们在做定位的时候，比如，计算模板匹配的亚像素位置。传统的像素级定位往往无法满足高精度需求，而亚像素定位技术通过数学插值和曲面拟合，可以实现超越单个像素的定位精度。

本文将详细介绍如何使用C++对3×3图像区域进行二次曲面拟合，求解亚像素极值位置及其像素值。本文主要针对这个模板匹配位置进行亚像素位置拟合来展开说明。

原理介绍：

对于计算亚像素位置，一般有两种方法；线性插值方法以及拟合的方法。

这篇文章主要介绍拟合的方法。

1.一维数据拟合（单行像素拟合）

我们都知道，像素是离散分布的，因此对于其位置，只能获取到其像素级别的位置。

若要获取真正的极值点，需要根据几个离散点进行函数拟合，然后计算极值点。

2.二维数据拟合

对于一维数据，可以将其拟合成线进行亚像素极值点计算；但是对于二维数据，可能就需要将其拟合成面了。但是怎么去计算面上的极值点，可能就需要求解特殊函数了。

而在模板匹配计算处理之后，会得到一个匹配得分图像，可通过取 极大值点3*3区域内的点进行曲面拟合，进而可计算出亚像素的极大值位置，而曲面拟合之后，怎么去计算这个极大值呢？

算法原理：

1.二次曲面模型

拟合3*3像素区域内的二次曲面函数为：

其中：

• (x,y) 是相对于中心点(0,0)的坐标

• a₀ 到 a₅ 是需要求解的系数

• 3×3网格坐标范围为 x ∈ [-1, 1], y ∈ [-1, 1]

2.最小二乘拟合

这里借鉴了这位博主：最小二乘法曲面拟合及C++代码，用于求取图像定位的亚像素级精度_曲面拟合求亚像素点位置-CSDN博客

并使用最小二乘法进行求解：是下面式子的误差最小：
则应该满足：







写为矩阵的表达式：

下一步就是要根据拟合点的坐标信息，解出[a,b,c,d,e,f]。在本例子中，共选取9个点，对曲面进行拟合（具体可以根据实际情况进行选择）。9点的坐标可以表示为。

代码：

#include <iostream>
#include <cmath>
#include <vector>

// 求解线性方程组 Ax = b
bool solveLinearSystem(const std::vector<std::vector<double>>& A, 
                       const std::vector<double>& b, 
                       std::vector<double>& x) {
    int n = A.size();
    if (n == 0) return false;
    
    // 创建增广矩阵
    std::vector<std::vector<double>> Ab(n, std::vector<double>(n + 1));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            Ab[i][j] = A[i][j];
        }
        Ab[i][n] = b[i];
    }
    
    // 高斯消元（部分主元）
    for (int col = 0; col < n; ++col) {
        // 寻找主元
        int max_row = col;
        for (int row = col + 1; row < n; ++row) {
            if (std::abs(Ab[row][col]) > std::abs(Ab[max_row][col])) {
                max_row = row;
            }
        }
        
        // 交换行
        if (max_row != col) {
            std::swap(Ab[col], Ab[max_row]);
        }
        
        // 检查主元是否为零
        if (std::abs(Ab[col][col]) < 1e-10) {
            return false;
        }
        
        // 消元
        for (int row = col + 1; row < n; ++row) {
            double factor = Ab[row][col] / Ab[col][col];
            for (int k = col; k <= n; ++k) {
                Ab[row][k] -= factor * Ab[col][k];
            }
        }
    }
    
    // 回代
    x.resize(n);
    for (int i = n - 1; i >= 0; --i) {
        x[i] = Ab[i][n];
        for (int j = i + 1; j < n; ++j) {
            x[i] -= Ab[i][j] * x[j];
        }
        x[i] /= Ab[i][i];
    }
    return true;
}

// 计算亚像素极值位置和值
bool computeSubpixelExtremum(const double pixels[3][3], 
                            double& dx, double& dy, double& extremumValue) {
    // 定义3x3网格坐标 (x,y) ∈ [-1,1]×[-1,1]
    const double positions[9][2] = {
        {-1, -1}, {0, -1}, {1, -1},
        {-1,  0}, {0,  0}, {1,  0},
        {-1,  1}, {0,  1}, {1,  1}
    };
    
    // 构建设计矩阵A和观测向量b
    std::vector<std::vector<double>> A(9, std::vector<double>(6));
    std::vector<double> b(9);
    
    for (int i = 0; i < 9; ++i) {
        double x = positions[i][0];
        double y = positions[i][1];
        int row = static_cast<int>(y + 1);  // 从y坐标到行索引
        int col = static_cast<int>(x + 1);  // 从x坐标到列索引
        
        A[i][0] = 1.0;
        A[i][1] = x;
        A[i][2] = y;
        A[i][3] = x * y;
        A[i][4] = x * x;
        A[i][5] = y * y;
        
        b[i] = pixels[row][col];
    }
    
    // 构建正规方程 A^T A x = A^T b
    std::vector<std::vector<double>> ATA(6, std::vector<double>(6, 0.0));
    std::vector<double> ATb(6, 0.0);
    
    for (int k = 0; k < 6; ++k) {
        for (int j = 0; j < 6; ++j) {
            for (int i = 0; i < 9; ++i) {
                ATA[k][j] += A[i][k] * A[i][j];
            }
        }
        for (int i = 0; i < 9; ++i) {
            ATb[k] += A[i][k] * b[i];
        }
    }
    
    // 求解系数向量 [a0, a1, a2, a3, a4, a5]
    std::vector<double> coeffs;
    if (!solveLinearSystem(ATA, ATb, coeffs) || coeffs.size() != 6) {
        return false;
    }
    
    double a1 = coeffs[1];
    double a2 = coeffs[2];
    double a3 = coeffs[3];
    double a4 = coeffs[4];
    double a5 = coeffs[5];
    
    // 计算分母 (4*a4*a5 - a3^2)
    double denominator = 4 * a4 * a5 - a3 * a3;
    
    // 检查分母是否接近零
    if (std::abs(denominator) < 1e-10) {
        return false;
    }
    
    // 计算亚像素偏移
    dx = (a3 * a2 - 2 * a5 * a1) / denominator;
    dy = (a3 * a1 - 2 * a4 * a2) / denominator;
    
    // 检查偏移是否在合理范围内
    if (std::abs(dx) > 1.0 || std::abs(dy) > 1.0) {
        return false;
    }
    
    // 计算亚像素极值
    extremumValue = coeffs[0] + 
                   coeffs[1] * dx + 
                   coeffs[2] * dy + 
                   coeffs[3] * dx * dy + 
                   coeffs[4] * dx * dx + 
                   coeffs[5] * dy * dy;
    
    return true;
}

int main() {
    // 示例输入：3x3像素矩阵
    double pixels[3][3] = {
        {0.2, 0.3, 0.1},
        {0.4, 0.9, 0.5},
        {0.1, 0.4, 0.3}
    };
    
    double dx, dy, extremumValue;
    if (computeSubpixelExtremum(pixels, dx, dy, extremumValue)) {
        std::cout << "亚像素极值位置偏移: (" << dx << ", " << dy << ")\n";
        std::cout << "相对于中心点的位置: (" << 1.0 + dx << ", " << 1.0 + dy << ")\n";
        std::cout << "亚像素极值: " << extremumValue << std::endl;
    } else {
        std::cout << "无法计算亚像素极值" << std::endl;
    }
    
    return 0;
}