数字图像处理之积分图计算优化

上一篇文章讲解了积分图的计算原理，并使用C++实现了基本的积分图算法。在WIN10 i5-7500 CPU的运行环境下，计算一帧1024*1024的图像，耗时7 ms左右。本文将在此基础上做进一步计算优化以减少计算耗时。

如上图所示，在原算法中计算I(x,y)需要I(x-1,y-1)，I(x,y-1)和I(x-1,y)，而在优化算法中要计算I(x,y)，只需要计算I(x,y-1)和点(x,y)所在行的前缀和即可，表示为如下公式。从公式看运算量很大，因为每次都要计算前缀和(x+1次加法)，实际上在代码实现过程中，可以使用一个变量来累加前缀和，每次在前面累加的基础上再加上点(x,y)的像素值即可(1次加法)。

计算过程中分为两步，第一步计算第一行的积分值，也即前缀和，第二步计算图像所有剩余像素点的积分值，C++实现代码如下。该算法计算1024*1024图像的积分图，耗时2 ms左右。

void Integal_row(Mat src, Mat &integal_out)
{
    src.convertTo(integal_out, CV_32F);
    float *p = integal_out.ptr<float>(0);
    for(int i = 1; i < src.cols; i++)   //计算第一行像素点的积分值
    {
        p[i] += p[i-1];
    }
    float *p1;
    for(int i = 1; i < src.rows; i++)   //从第二行开始，第i行
    {
        float sum = 0.0;   //累加和变量
        p = integal_out.ptr<float>(i);
        p1 = integal_out.ptr<float>(i-1);
        for(int j = 0; j < src.cols; j++)   //第j列
        {
            sum += p[j];      //累加当前像素值
            p[j] = p1[j] + sum;    //计算当前像素点的积分值
        }
    }


}

在以上基础上，再进一步使用SSE指令进行优化。数据流单指令序列扩展指令，简称为SSE指令，是Intel公司最早在Pentium III芯片中引入的指令集，SSE指令集是MMX指令集的扩展。SSE指令是一种“单指令多数据”的指令集，即在一个指令周期内同时对多个数据进行相同运算的指令集，因此可以有效提升CPU的数据运算效率。自从Visual Studio.NET 2003增加了对SSE指令集的编程支持，在VisualStudio开发环境中，开发人员在C/C++代码中不再需要嵌入汇编代码，只需调用已封装好的SSE指令函数库就可以直接使用SSE指令功能。

SSE指令的基本数据类型有_m128（单精度浮点数）、_m128d（双精度浮点数）和_m128i（整型）等，这些数据类型的二进制位数都是128位，区别在于包含的子数据类型不相同：

1. _m128包含四个单精度浮点数，如图4.5所示。                        

2. _m128d包含两个双精度浮点数，如图4.6所示：

3. _m128i数据类型的定义是一个C/C++共用体，其可以被当作不同位数的整型数进行读写，比如可以分别当作16个8位、8个16位、4个32位、2个64位的整型数，其中不同位数的整型数又可以分为有符号数和无符号数。

在C/C++中进行SSE指令运算的基本流程是先把C/C++的基本数据类型数据载入SSE的基本数据类型变量，然后调用SSE库函数对变量进行运算，运算完毕后再把SSE变量中的运算结果导出到C/C++的基本数据类型变量。SSE库函数的调用方式非常简单，只需要包含相应的头文件即可调用：

#include <xmmintrin.h>
#include <emmintrin.h>
#include <immintrin.h>

首先，把以上代码的内循环展开，以便进行SSE优化：

void Integal_row(Mat src, Mat &integal_out)
{
    src.convertTo(integal_out, CV_32F);    //将uchar型转换为float型
    float *p = integal_out.ptr<float>(0);
    for(int i = 1; i < src.cols; i++)   //计算第一行像素点的积分值
    {
        p[i] += p[i-1];
    }
    
    float *p1;


    for(int i = 1; i < src.rows; i++)   //第i行


    {


        float sum = 0.0;


        p = integal_out.ptr<float>(i);


        p1 = integal_out.ptr<float>(i-1);


        int j = 0;


        for(; j < src.cols-4; j+=4)   //每次循环计算4个像素点的积分值


        {
            //展开循环部分
            sum += p[j];


            p[j] = p1[j] + sum;


            sum += p[j+1];


            p[j+1] = p1[j+1] + sum;


            sum += p[j+2];


            p[j+2] = p1[j+2] + sum;


            sum += p[j+3];


            p[j+3] = p1[j+3] + sum;


        }


        for(; j < src.cols; j++)   //计算该行剩余点的积分值


        {


            sum += p[j];


            p[j] = p1[j] + sum;


        }


    }


}

接着，使用SSE指令优化内循环的展开部分：

void Integal_row(Mat src, Mat &integal_out)
{
    src.convertTo(integal_out, CV_32F);    //将uchar型转换为float型
    float *p = integal_out.ptr<float>(0);
    for(int i = 1; i < src.cols; i++)   //计算第一行像素点的积分值
    {
        p[i] += p[i-1];
    }
    
    float *p1;


    for(int i = 1; i < src.rows; i++)   //第i行


    {


        float sum = 0.0;


        p = integal_out.ptr<float>(i);


        p1 = integal_out.ptr<float>(i-1);


        int j = 0;


        for(; j < src.cols-4; j+=4)   //每次循环计算4个像素点的积分值


        {


                float a1 = p[j] + p[j+1];   
                float a2 = a1 + p[j+2];
                float a3 = a2 + p[j+3];
                //_mm_set1_ps为将同一个float数据x载入__m128变量中，即x x x x
               //_mm_set_ps为将4个float数据载入__m128变量中，即x3 x2 x1 x0
               //_mm_add_ps为相加两个__m128变量，即x+x3 x+x2 x+x1 x+x0
                __m128 sum_rst = _mm_add_ps(_mm_set1_ps(sum), _mm_set_ps(a3, a2, a1, p[j]));
                 __m128 pi = _mm_add_ps(_mm_set_ps(p1[j+3], p1[j+2], p1[j+1], p1[j]), sum_rst);
                _mm_storeu_ps(&p[j], pi);      //将计算结果从__m128变量下载到数组中
                sum += a3;        //将本次循环的像素点的像素值累加到sum中
        }


        for(; j < src.cols; j++)   //计算该行剩余点的积分值


        {


            sum += p[j];


            p[j] = p1[j] + sum;


        }


    }


}

使用SSE指令优化之后，耗时由2 ms左右减小到1.3 ms左右。积分图结果如下：