【NCC】之二：积分图加速均值计算

积分图 integral image

1. 原理：

Summed Area Table是一种数据结构和算法，用于快速有效地生成网格矩形子集中的值总和。在图像处理领域，它也被称为积分图像。

记源图像为$I_{m\times n}$，积分图为SAT（Summed Area Table）则
$$SAT(x,y)=\underset{i\le x,j\le y}{\Sigma }I(i,j)$$

$$SQAT(x,y)=\underset{i\le x,j\le y}{\Sigma }I(i,j{)}^2$$

• 计算积分图的递推式
  $$SAT(x,y)=SAT(x-1,y)+SAT(x,y-1)-SAT(x-1,y-1)+I(x,y)$$
  上面的公式中起始项必须从第二行第二列开始，因为递推式中包含-1项，也就是需要单独把第一行和第一列先单独计算出来，后面的所有点都可以用递推式写出来。

• 一点漂亮的改进
  让积分图比原图多一行一列：上面一行，左边一列，且行列初始值为0！这样的话就可以不需要考虑-1项完全采用递推式得到积分图！opencv就是这么干的

2. 示例

图片摘自Wikipedia。1为源图，2为积分图每一个像素点的值都是它在源图上对应位置的上面、左边所有像素值的和。

3. 计算区域均值

当要计算源图中的某个矩形区域的像素和时可以用查表的方式快速计算，构建好积分图后求区域和的复杂度将是$O(1)$。

设在源图中的矩形$R$左上角坐标为$(tpx,tpy)$，右下角坐标为$(btx,bty)$,则计算该矩形区域内所有像素值和可以写为下式。

• 区域像素值之和
  $$sum=SAT(btx,bty)-SAT(tpx-1,bty)-SAT(btx,tpy-1)+SAT(tpx-1,tpy-1)$$
• 区域均值
  $$mean=\frac{sum}{(bty-tpy+1)(btx-tpx+1)}$$

4. 计算区域方差

以【NCC】之一：从Pearson相关系数到模板匹配的NCC方法中的$S_{x,y}$的方差计算为例
计算区域方差
$$\sigma (S_{x,y})=\sqrt{var(S_{x,y})}=\sqrt{\frac{{\Sigma }_{i=1}^m{\Sigma }_{j=1}^n(S_{x,y}(i,j)-\overline{S_{x,y}}{)}^2}{mn}}$$
let $A={\Sigma }_{i=1}^m{\Sigma }_{j=1}^n(S_{x,y}(i,j)-\overline{S_{x,y}}{)}^2$ 则有
$$\begin{array}{r}\begin{array}{rl}A& ={\Sigma }_{i=1}^m{\Sigma }_{j=1}^n(S_{x,y}(i,j{)}^2-2\overline{S_{x,y}}{S}_{x,y}(i,j)+{\overline{S_{x,y}}}^2)\\ & =\underset{S_{x,y}}{SQAT(m,n)}-\underset{S_{x,y}}{2\overline{S_{x,y}}\times SAT(m,n)}+mn\times {\overline{S_{x,y}}}^2\\ & =\underset{S_{x,y}}{SQAT(m,n)}-\underset{S_{x,y}}{2\overline{S_{x,y}}\times SAT(m,n)}+\overline{S_{x,y}}\times SAT(m,n)\\ & =\underset{S_{x,y}}{SQAT(m,n)}-\underset{S_{x,y}}{\overline{S_{x,y}}\times SAT(m,n)}\end{array}\end{array}$$

5. 积分图示例

• 原图
  

• 积分图
  

• 平方积分图
  

6. 计算积分图的源码

/**
 * @brief 计算输入图的积分图,为了提高计算效率，可以让积分图比输入图多一行一列，
 * 具体的就是在原图左边插入一列0，上面插入一行0，设原图为I，积分图为SAT(sum area of table)
 * 则：SAT(i,j) = SAT(i,j-1) + SAT(i-1,j) - SAT(i-1,j-1) + I(i,j)
 * 这样就不用考虑下边界的情况，省掉很多判断条件
 * 
 * @param image  : 输入图CV_8UC1，MxN
 * @param integral_image  : 积分图CV_32FC1,(M+1)x(N+1)
 * @return int 
 */
int integral(const cv::Mat &image,cv::Mat &integral_image)
{
    if(image.empty())
    {
        MYCV_ERROR(kImageEmpty,"image empty");
        return kImageEmpty;
    }

    int h = image.rows;
    int w = image.cols;
    integral_image = cv::Mat::zeros(cv::Size(w+1,h+1),CV_64FC1);

    //SAT(i,j) = SAT(i,j-1) + SAT(i-1,j) - SAT(i-1,j-1) + I(i,j)
    for(int i = 0; i < h ; i++)
    {
        const uchar *ps = image.ptr<uchar>(i);
        double *pd_m1 = integral_image.ptr<double>(i);//integral 的"上一行"
        double *pd = integral_image.ptr<double>(i+1); //integral 的"当前行"
        for(int j = 0; j < w; j++)
        {
            pd[j+1] = pd[j] + pd_m1[j+1] - pd_m1[j] + ps[j];
        }
    }

    return kSuccess;
}


/**
 * @brief 计算输入图的积分图,为了提高计算效率，可以让积分图比输入图多一行一列，
 * 具体的就是在原图左边插入一列0，上面插入一行0，设原图为I，积分图为SAT(summed area table)
 * 则：SAT(i,j) = SAT(i,j-1) + SAT(i-1,j) - SAT(i-1,j-1) + I(i,j)
 * SQAT(i,j) = SQAT(i,j-1) + SQAT(i-1,j) - SQAT(i-1,j-1) + I(i,j) * I(i,j)
 * 这样就不用考虑下边界的情况，省掉很多判断条件
 * 
 * @param image  : 输入图CV_8UC1，MxN
 * @param integral_image  : 积分图CV_32FC1,(M+1)x(N+1)
 * @param integral_sq : 平方的积分图CV_32FC1,(M+1)x(N+1)
 * @return int 
 */
int integral(const cv::Mat &image,cv::Mat &integral_image,cv::Mat &integral_sq)
{
     if(image.empty())
    {
        MYCV_ERROR(kImageEmpty,"image empty");
        return kImageEmpty;
    }

    int h = image.rows;
    int w = image.cols;
    integral_image = cv::Mat::zeros(cv::Size(w+1,h+1),CV_64FC1);
    integral_sq = cv::Mat::zeros(cv::Size(w+1,h+1),CV_64FC1);

    //SAT(i,j) = SAT(i,j-1) + SAT(i-1,j) - SAT(i-1,j-1) + I(i,j)
    //SQAT(i,j) = SQAT(i,j-1) + SQAT(i-1,j) - SQAT(i-1,j-1) + I(i,j) * I(i,j)
    for(int i = 0; i < h ; i++)
    {
        const uchar *ps = image.ptr<uchar>(i);
        double *pd_m1 = integral_image.ptr<double>(i);//integral 的"上一行"
        double *pd = integral_image.ptr<double>(i+1); //integral 的"当前行"
        double *pqd_m1 = integral_sq.ptr<double>(i);
        double *pqd = integral_sq.ptr<double>(i+1);
        for(int j = 0; j < w; j++)
        {
            pd[j+1] = pd[j] + pd_m1[j+1] - pd_m1[j] + (double)ps[j];
            pqd[j+1] = pqd[j] + pqd_m1[j+1] - pqd_m1[j] + (double)ps[j] * (double)ps[j];
        }
    }


    return kSuccess;
}

7. 用积分图加速NCC

采用积分图计算均值和方差

source image size w,h = (1095,680)
target image size w,h = (89,91)
my NCC run 5 times,average use 5699.000000 ms
min_value=-0.431418 , min_loc(x,y)=(799,210),    max_value=0.998323,max_loc(x,y)=(393,286)
opencv NCC run 5 times,average use 20.000000 ms
min_value=-0.431418 , min_loc(x,y)=(799,210),    max_value=0.998322,max_loc(x,y)=(393,286)

opencv的速度是这一版本的285倍，相比上一版882.78倍还是有点进步。

8. 补充

测试自己的积分图和opencv的积分图在数值上完全一致了，但是速度上opencv的是我的3-5倍，今天（2023年6月5日）采用IPP中的ippiSqrIntegral_8u32f64f_C1R()方法后速度就和opencv中的积分图一致了。IPP是不开源的，IPPICV是一个IPP提取出来的一个子库专门用在图像处理里面，现在已经打包在opencv中，在IPPICV中对应的函数为 ippicviSqrIntegral_8u32f64f_C1R

参考

[1] https://en.wikipedia.org/wiki/Summed-area_table
[2] 积分图(一) - 原理及应用