图像编程要点，如何加速对图像的处理

图像本身数据特点

图像常用矩阵形式进行储存；但图像本身的数据量是极其大的。以1080P视频为例，每秒60帧1920*1080的彩色图像，原始字节数高达460M每秒。但图像本身存在一些规律，因此形成了独特的处理方法。
存储特点：
凡谈到图像，一般都是指的是一个二维的矩阵(数组)，其在计算机内存的存放是一个连续的地址空间，该地址空间可以由第一象素和最末一个象素的存储地址决定，也可由第一象素和总的象素个数决定。二维数组的各元素的在内存中的存放顺序是按行存储的，即在此连续地址的内存空间中先存放第一行的所有数据，再存放第二行的内存数据，依次类推，直到图像的最后一行。数组在计算内存中是占用一块连续的地址空间，因此可以当做一维数组来访问。
图像表示：
原始图像一般用pImg命名，有时灰度图像、彩色图像、二值图像分
别用pGryImg、 pRGBImg、 pBinImg命名，处理结果的图像一般用pResImg命名， p代表指针数据类型；图像的宽度为width，图像的高度为height

图像矩阵访问形式

由于存储在连续内存空间中，应充分利用指针访问的便利性，无需自己算计在数组中的位置。
点运算访问：

BYTE *pCur,*pEnd,*pRes;
for(pCur=pImg,pEnd=pImg+width*height,pRes=pResImg;pCur<pEnd;pCur++)
{
*(pRes++)=f(*pCur); //f代表点运算的语句体
}

邻域访问：

int x,y;
BYTE *pCur,*pRes;
for(y=0,pCur=pImg,pRes=pResImg;y<height;y++)
{
for(x=0;x<width;x++,pCur++,pRes++)
{
*pRes=f(pCur,x,y); //f代表邻域运算的语句体(与位置x,y有关)
}
}

不建议访问方式：

for(row=0;row<height;row++)
{
for(col=0;col<width;col++)
{
index=row*height+col; //此处多了乘法和加法运算
pResImg[index]=f(pImg[index]); //f代表运算的语句体
}
}

for(y=0;y<height;y++)
{
for(x=0;x<width;x++)
{
pResImg[y][x]=f(pImg[y][x];
//pResImg[x][y]=f(pImg[x][y]; //会更加无效，为什么？
}
}
此种形式，若不考虑C编译器对程序的优化， pImg[y][x]实质上是
*(pImg+y*width+x)，其中包含了一个乘法和两个加法。

Cache访问

高速缓冲存储器Cache是位于CPU与内存之间的临时存储器，容量比内存小， 交换速度快。在Cache中的数据是内存中的一小部分，这一小部分是短时间内CPU即将访问的数据。 当CPU调用大量数据时，就可避开内存直接从Cache中调用，从而加快读取速度。
Cache对CPU的性能影响很大，主要是由CPU的数据交换顺序和CPU与Cache间的带宽引起的。

这样的读取机制使CPU读取Cache的命中率非常高（大多数CPU可达90%左右），也就是说CPU下一次要读取的数据90%都在Cache中，只有大约10%需要从内存读取。这大大节省了CPU直接读取内存的时间，也使CPU读取数据时基本无需等待。总的来说，CPU读取数据的顺序是先Cache后内存。
已经知道， 数据在Cache中时的速度最高，那么若能通过合理的机制， 增加访问
的数据在Cache中的概率，则能大大提高CPU访问数据的速度，从而提高程序的执
行速度。
由于计算机的Cache相对于数据量动辄几百万的图像而言还是太小了，不足以放
下一幅完整的图像，常规内存的数据块将会频繁出入Cache， 因此若能降低数据
块进出Cache的次数，则能提高速度。
按行访问图像数据:
结合图像数据的存放特点，显然，程序按行访问图像数据肯定比按列访问图像数
据的执行速度快；

注： 算法复杂度与执行速度并不完全对应，如果算法不便于缓冲区优化，即使算
法复杂度再低，也许执行速度并不快。 图像处理的速度受到缓冲区命中率和内存
速度的影响不可小视。
将图像分块或分多种分辨率处理:
将一个大图像分成宽度不等的N个小图像，对N个小图像分别处理所花费的总时间
肯定小于对大图像处理的时间。因为小图像中每行的数据更少，一次内存中的数
据块进入Cache将进入更多的行数，从而提高了数据在Cache中的命中率。 (当然
，要评估一下图像分成N个小图像所花费的代价，才能决定值得不值得)。
同样的原理， 如果一个图像处理的程序能够按照由粗到精的策略进行，则先把原
图像缩小M倍，进行粗略的处理(比如目标定位等)；然后在原图像中取出感兴趣
的区域(肯定远远小于原图像)，再进行精确的处理。这样，由于处理的都是小图
像，也就提高了数据在Cache中的命中率.

LUT查找表

在空间域的图像计算一般可以分为点运算和邻域运算两大类。点运算是指运算后的
结果仅仅与该单个像素的灰度值有关，而与该像素所处的位置无关,比如线性拉伸
、直方图均衡化等；邻域运算则是运算后的结果与该像素所处的位置有关，比如均
值滤波、中值滤波、边缘检测等。
所有点运算均可表示为一个函数F: G=F(g)； g是像素的灰度值(自变量)， G是变换
后的值(因变量)。
然而一幅图像中的像素数量非常可观，动辄几百万像素；因此F虽看起来简单，但
一幅图像就几百万次的函数调用，其实很浪费时间。
考虑到图像的特点：数据量特别大，灰度值范围非常小； 比如通常的灰度图像
0<=g<=255最多只有256种取值。因此，任意以图像灰度为自变量的函数，都可以使
用查找表（LUT , Look Up Table） 的方法实现。
转换规则与一般形式：

int LUT[256],g；
for(g=0;g<256;g++) LUT[g] =max(0,min(255,F(g));

举例：

图像灰度Gama校正
void RmwGryImageGamaCorrect_Slow(BYTE *pGryImg,double gama,int nSize)
{
for(int i=0;i<nSize;i++)
{
*(pGryImg+i) = min(255,pow(*(pGryImg+i),gama));
}
return;
}
void RmwGryImageGamaCorrect_Fast(BYTE *pGryImg,double gama,int nSize)
{
int LUT[256],i;
for(i=1;i<256;i++) LUT[i]=min(255,pow(i,gama));
for(i=0;i<nSize;i++)
{
*(pGryImg+i) = LUT[*(pGryImg+i)];
}
return;
}
当nSize=768*576时,在IBMX31笔记本上， RmwGryImageGamaCorrect_Slow用时
114ms， RmwGryImageGamaCorrect_Fast用时1.6ms。 (VC6.0用release版编译)(速度提高了71倍)。

彩色图像到灰度图像转换的典型程序
void RmwRGB2GryImg(BYTE *pRGBImg,int width,int height,BYTE *pGryImg)
{
BYTE *pRGB,*pGry,*pEnd=pRGBImg+3*width*height;
double gry;
for(pRGB=pRGBImg,pGry=pGryImg;pRGB<pEnd;)
{
gry = *(pRGB++)*0.114; //B*0.114
gry += *(pRGB++)*0.587; //G*0.587
gry += *(pRGB++)*0.299; //R*0.229
*(pGry++) = (int)(gry);
}
return;
}

彩色图像到灰度图像转换使用查找表
void RmwRGB2GryImgByLUT(BYTE *pRGBImg,int width,int height,BYTE *pGryImg)
{
BYTE *pRGB,*pGry,*pEnd=pRGBImg+3*width*height;
int LUTR[256],LUTG[256],LUTB[256];
int i,gry;
for(i=0;i<256;i++)
{
LUTR[i]=(int)(0.299*i*1024);
LUTG[i]=(int)(0.587*i*1024);
LUTB[i]=(int)(0.114*i*1024);
}
for(pRGB=pRGBImg,pGry=pGryImg;pRGB<pEnd;)
{
gry = LUTB[*(pRGB++)];
gry += LUTG[*(pRGB++)];
gry += LUTR[*(pRGB++)];
*(pGry++) = gry>>10;
}
return;
}

Sobel算子的典型程序
void RmwSobelGryImg(BYTE *pGryImg,int width,int height,BYTE *pSbImg)
{
int x,y,g1,g2;
BYTE *pGry,*pSb;
memset(pSbImg,0,width); //首行不做
for(y=1,pGry=pGryImg+width,pSb=pSbImg+width;y<height-1;y++)
{
pGry++;
*(pSb++)=0; //首列不做
for(x=1;x<width-1;x++)
{
g1 = *(pGry-width-1)+ (*(pGry-width)*2) + *(pGry-width+1);
g1 -= *(pGry+width-1)+ (*(pGry+width)*2) + *(pGry+width+1);
g2 = *(pGry-1-width)+ (*(pGry-1)*2) + *(pGry-1+width);
g2 -= *(pGry+1-width)+ (*(pGry+1)*2) + *(pGry+1+width);
*(pSb++)=min(255,abs(g1)+abs(g2)); //此步需要优化
}
pGry++;
*(pSb++)=0; //尾列不做
}
memset(pSb,0,width); //尾行不做
}

Sobel算子的快速实现
void RmwGryImgSobelLUT(BYTE *pGryImg,int width,int height,BYTE *pSbImg)
{ //注意：需要10K的查表空间
//绝对值表
int absLUT[2048],*pAbs;
//求和表
int sumLUT[512];
//运算
BYTE *pCur,*pPre,*pNxt,*pSb;
int g1,g2;
int x,y;
int i;
// step.1--------初始化绝对值查找表,超过255取255-----------------//
pAbs=absLUT+1024; //梯度范围为[-255*4,255*4]
for(i=0;i<256;i++) pAbs[i]=pAbs[-i]=i;
for(i=256;i<1024;i++) pAbs[i]=pAbs[-i]=255;
// step.2--------初始化求和表,超过255取255-----------------------//
for(i=0;i<256;i++) sumLUT[i]=i; //范围为[0,255*2]
for(i=256;i<512;i++) sumLUT[i]=255;
/ / step.3--------Sobel,注意查表和饱和运算------------------------//
memset(pSbImg,0,width); //首行不做
for(y=1,pCur=pGryImg+width,pPre=pCur-width, pNxt=pCur+width,
pSb=pSbImg+y*width;
y<height-1;y++
)
{
//首列不做
pCur++; pPre++; pNxt++;
*(pSb++)=0;
//中间
for(x=1;x<width-1;x++,pCur++,pPre++,pNxt++,pSb++)
{
g1 = *(pPre-1) + (*pPre)*2 + *(pPre+1) - *(pNxt-1) - (*pNxt)*2 - *(pNxt+1);
g2 = *(pPre-1) + *(pCur-1)*2+ *(pNxt-1) - *(pPre+1) - *(pCur+1)*2- *(pNxt+1);
*pSb=sumLUT[pAbs[g1]+pAbs[g2]];//*pSb=min(255,abs(g1)+abs(g2));
}
//尾列不做
pCur++; pPre++; pNxt++;
*(pSb++)=0;
}
memset(pSb,0,width); //尾行不做
}

查找表实际上是一种以空间换时间的策略。 通常图像处理系统都有成组的查找表供
编程使用。 在通用的计算机上， 没有查找表可以直接利用， 则通过一个线性数组来
实现， 这时查找表仅表现为一种数据结构。
同时， 查找又是一个间接寻址的过程， 每个像素均增加了一次计算机对内存的访问
， 有时当计算机内存的吞吐能力很小时， 反而降低了速度， 即查表不如直接计算更
快。 在很多DSP中， 可以把表放在一级Cache中(L1中) ， 更能提高速度。 PC的高速缓存是操作系统管理的， 用户很难确定地知道表是否在高速缓存中。 查表量和计算量， 要达到最佳的平衡， 才能提高程序的执行速度。
查找表适合于表达自变量范围很小时的函数。 图像处理中的灰度变换、 图像增强、
图像空间的变换、 梯度运算、 差图像等等， 都可借鉴查找表。
查找表太大时， 会非常浪费空间， 这时可以使用部分查找表， 来平衡时间和空间的
效率。 比如梯度算子中的开方运算， y=sqrt(x),x的值域为[0， 255*255*2],此时若
构造平方根表， 则会非常浪费空间；考虑到x一般都是小于1000的， 可以判断x小于
1000就查表， 否则就求平方根。
查找表一般都是一维的， 二维表由于占内存空间大和寻址不连续导致缓冲区效率下
降和增加额外的计算量而不用。

Histogram（直方图）

直方图（histogram）是灰度级的函数，反映图像中每种灰度出现的频率。由于频
率是浮点数，其与频数之间只是相差一个归一化常数，所以图像编程中直方图就是
灰度级出现的频数，是个32位的长整数(防止计数上溢)。
直方图在图像处理中应用广泛： 图像增强中的直方图的均衡化、规定化，图像分割
中基于直方图的阈值选取，图像检索中的图像相似度测量。
从编程的观点看，直方图是一种很有效的数据结构： 所占内存空间很少（1024个字
节） ， 能反映出图像中的灰度分布和目标特性等等，基于直方图的有序性可以容易
高效地得到图像的亮度、对比度、最大亮度、最小亮度及亮度中值。
直方图的常用编程：

直方图的构造
void RmwGetHistogram(BYTE *pImg,int imgSize,int *histogram)
{
BYTE *pCur, *pEnd;
//for(int i=0;i<256;i++) histogram[i]=0;
memset(histogram,0,sizeof(int)*256);
for(pCur=pImg, pEnd=pImg+ImgSize; pCur<pEnd;)
histogram[*(pCur++)]++;
return;
}

图像的亮度和对比度
void RmwBrightContrast(int *histogram,double *bright, double *contrast)
{
int g,sumGry,imgSize;
for(sumGry=imgSize=0,g=0;g<256;g++)
{
sumGry+=histogram[g]*g;
imgSize += histogram[g];
}
*brightness=1.0*sumGry/imgSize; //亮度（平均值）
for(sumGry =0,g=0;g<256;g++)
{
if (histogram[g])
{
sumGry+=histogram[g] * (g-*brightness) * (g-*brightness);
}
}
*contrast = sqrt(sumGry/(imgSize-1.0)); //对比度（方差）
return;
}

图像的最大亮度
void RmwImageMaxGry(int *histogram, int *maxGry)
{ int i;
for(i=255;i>=0;i--) //从最大值开始找
if (histogram[i]) break; // 灰度的最大值
*maxGry=i;
return;
}

图像的最小亮度
void RmwImageMinGry(int *histogram, int *minGry)
{
int i;
for(i=0;i<256;i++)
if (histogram[i]) break; //灰度的最小值
*minGry=i;
return;
}

图像的亮度中值
void RmwImageMedianGry(int *histogram, int imgSize, int *medGry)
{
int i;
for(sum=0,i=0; i<256; i++)
{
sum += histogram[i];//直方图本来就是按照灰度排好顺序的
//if (sum>imgSize/2) break;
if (sum*2>imgSize) break; //序列的一半就是中值
}
*medGry=i;
return;
}

直方图均衡化
void RmwHistogramEqualize(BYTE *pImg,int width,int height)
{
BYTE *pCur,*pEnd=pImg+width*height;
unsigned int hist[256];
int LUT[256],i,sum;
memset(hist,0,sizeof(int)*256);
for(pCur=pImg;pCur<pEnd;) hist[*(pCur++)]++;
for(sum=hist[0],LUT[0]=0,i=1;i<256;i++)
{
sum = sum+hist[i];
LUT[i]=255*sum/(width*height);
}
for(pCur=pImg;pCur<pEnd;) *(pCur++)=LUT[*pCur];
return;
}

应用举例: 中值滤波:
直方图本身的有序特性使得求图像的中值非常简单，大大减少比较运算的次数，
与图像的大小无关，最多比较256次。因此可以利用直方图，实现高效的图像中
值滤波。
另外，在滑动窗口滤波中，对于相邻像素的两个滑动窗口而言，变化的内容只是
丢掉了最左边的列，而增加了一个新的右边的列，窗口中间大面积的重叠区中的
像素并没有变化，并不需要重新访问和处理。
应用举例: 中值滤波
假定中值滤波的窗口为（nx*2+1） *（ny*2+1） ， 即当nx=1， ny=1时， 滤波窗口
就是3x3。 pGryImg原始灰度图像， pResImg是中值滤波后的图像， width是图像的
宽度， height是图像的高度。
不考虑边界和换行处理（是需要费心考虑的， 此处忽略） ， 一个简单的中值滤波
程序如下：

size=(nx*2+1)*(ny*2+1);
for(y=ny,pRes=pResImg+y*width;y<height-ny;y++)
{
//跳过左边界
pRes+=nx;
//step1.1初始化直方图
x=nx;
memset(histogram,0,sizeof(int)*256);
for(i=y-ny;i<=y+ny;i++)
for(j=x-nx;j<=x+nx;j++)
histogram[*(pGryImg+i*width+j)]++;
//step1.2初始化中值
for(num=0,med=0;i<256;med++)
{
num+=histogram[med];
if (num*2>=size) break;
}
//step1.3中值赋值
*(pRes++)=med;
//step.2对于该行中的每个像素(跳过左边界和右边界者)进行中值滤波
for(x=nx+1;x<width-nx;x++)
{
//step2.1更新直方图
for( i=y-ny,
pLeft=pGryImg+i*width+(x-nx-1),
pRight=pGryImg+i*width+(x+nx);
i<=y+ny;
i++,pLeft+=width,pRight+=width
)
{
//-Left
if (*pLeft<=med) num--;
histogram[*pLeft]--;
//+Right
if (*pRight<=med) num++;
histogram[*pRight]++;
}
//step2.2更新中值
if (num*2>size)
{
for(;med>=0;med--)
{ num -= histogram[med];
if (num*2<size) break;
}
num += histogram[med];
}
else if (num*2<size)
{
for(med=med+1;med<256;med++)
{ num += histogram[med];
if (num*2>=size) break;
}
}
//step2.3中值赋值
*(pRes++)=med;
} //end of x循环
//跳过右边界
pRes+=nx;
} //end of y循环
//跳过右边界
pRes+=nx;
} //end of y

Project（投影）

对图像处理稍微熟悉就可以知道，图像中的加法运算具有模糊的效果，即表现从整
体上看是怎么样，比如一个班级的成绩累加之和，突出的是这个班级的整体水平；
而图像中的减法运算突出的个体差别，是个体之间的差异比较，比如是张三同学比
李四同学成绩高了多少或者少了多少分。
图像中的均值滤波、高斯滤波都是典型的加法运算；边缘检测算子，无论一阶微分
算子还是二阶微分算子，则是典型的减法运算；当然，为了对噪声的抑制，边缘检
测算子中也往往会有加法运算，即“先平滑后求导”的思想，比如Sobel算子。
图像中还有一种典型的加法运算，即本节要讲的投影。 何谓投影，简单地说就是
在某个方向上的所有像素的灰度值之和，或者直线L上的所有像素的灰度之和。当
然，水平投影和垂直投影是两种典型投影。
自然界中的很多目标受到重力的影响、很多人造目标受到加工的限制，都具有典型
的水平或者垂直特征，比如车辆的尾部有一些水平线条，车头是流线型的，但也有
一些水平线条；树木都是直的；建筑都是方的等等；因为合理使用投影，能得到大
量的信息，且计算速度快，由于采用的是加法运算，性能又非常鲁棒。




梯度投影
按照基本的定义，投影是在某条直线上的图像像素灰度之和(成为灰度投影)。若图
像是个梯度图像呢？那么就不再是灰度之和，而是梯度之和(称为梯度投影)；若是
边缘图像呢？就是像素点数之和(称为边缘投影)。

IntegralImage（积分图）

积分图的概念最早是由Paul Viola等人(Detection Using a Boosted Cascade of
Simple Features. Computer Vision and Pattern Recognition, 2001, Volume
1, 8-14) 提出的，并被应用到实时的对象检测框架中。
对于一个灰度图像而言，其积分图也是一张图，只不过这个图跟普通的灰度图不
同,该图上任意一点(x,y)的值是指从灰度图像的左上角与当前点所围成的矩形区
域内所有像素的灰度值之和。

积分图的程序实现
void RmwDoSumGryImg(BYTE *pGryImg,int width,int height,int *pSumImg)
{
BYTE *pGry;
int *pSum;
int x,y;
// step.1-----------------先求第一行--------------------//
pGry=pGryImg;
pSum=pSumImg;
*(pSum++)=*(pGry++);
for(x=1;x<width;x++) *(pSum++)=*(pSum-1)+(*(pGry++));
// step.2-----------------再求其他行--------------------//
for(y=1;y<height;y++)
{
*(pSum++)=*(pSum-width)+(*(pGry++)); //先求第一列
for(x=1;x<width;x++) //再求其他列
{
*(pSum++)=*(pGry++)+(*(pSum-1))+(*(pSum-width))-(*(pSum-1-width));
}
}
}

应用举例: 均值滤波：

除法变成乘法和移位

图像间的除法运算
void RmwGryImageADivB_Slow(BYTE *pImgA,BYTE *pImgB,int width,int height)
{
int B;
BYTE *pA,*pB,*pEnd;
for(pA=pImgA,pEnd=pImgA+width*height,pB=pImgB;pA<pEnd;)
{
B=*(pB++);
if (B) *(pA++) = *pA/B;
else *(pA++)=255;
}
return;
}

void RmwGryImageADivB_Fast(BYTE *pImgA,BYTE *pImgB,int width,int height)
{
int B;
BYTE *pA,*pB,*pEnd;
int cLUT[256];
for(B=1;B<256;B++) cLUT[B]=(1<<22)/B;
for(pA=pImgA,pEnd=pImgA+width*height,pB=pImgB;pA<pEnd;)
{
B=*(pB++);
if (B) *(pA++) = (*pA*cLUT[B])>>22;
else *(pA++)=255;
}
return;
}

当width=768,height=576， 且pGryImgB图像中的灰度值均大于0 时,
RmwGryImageADivB_Slow用时6.2ms， RmwGryImageADivB_Fast用时1.6ms。
(VC6.0用release版编译)。

彩色图像到灰度图像转换使用乘法和移位
void RmwRGB2GryImgFast(BYTE *pRGBImg,int width,int height,BYTE *pGryImg)
{
BYTE *pRGB,*pEnd;
BYTE *pGry;
int sum;
pEnd=pRGBImg+3*width*height;
for(pRGB=pRGBImg,pGry=pGryImg;pRGB<pEnd;)
{
sum = *(pRGB++)*7471; //B*0.114
sum += *(pRGB++)*38469; //G*0.587
sum += *(pRGB++)*19595; //R*0.299
*(pGry++) = sum>>16;
}
return;
}

b*0.114 变化为 b*0.114*65536/65536, 即b*7471/65536， (b*7471） >>16
g*0.587 变化为 g*0.587*65536/65536, 即g*38469/65536， (g*38469） >>16
r*0.229 变化为 r*0.229*65536/65536, 即r*19595/65536， (r*19595） >>16

其他

在图像分析的编程中，还有更多的编程优化方法，比如使用SIMD指令集，在Intel
的CPU中采用MMX， SSE等指令，在ARM中使用Neon指令，都能大大提高程序的执行速
度。