数字图像处理之积分图计算优化
上一篇文章讲解了积分图的计算原理,并使用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左右。积分图结果如下: