图像本身数据特点
图像常用矩阵形式进行储存;但图像本身的数据量是极其大的。以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指令,都能大大提高程序的执行速
度。