2010-05-16 12:17:00 qytan36 阅读数 974
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4125 人正在学习 去看看 贾志刚

实验环境:
1,Linux操作系统
2,QT3编程开发环境
3,C++编程语言

傅立叶变换和傅立叶反 变换

1.1. 主要源代码

readImage() 从图像中读取数据
writeImage() 往图像中写入数据
fft() 快速傅立叶变换
ifft() 快速傅立叶反变换
adjustImageSize() 调整图像大小
fourier() 傅立叶变换
ifourier() 傅立叶反变换

1.1.1 从图像中读取数据

void ImageProcess::readImage(complex<double> data[], const QImage &srcImage)

{

byte *pImageBytes = srcImage.bits(); //数据首地址

int depth = srcImage.depth(); //每个像素的bit数

int lineBytes = srcImage.bytesPerLine(); //每行的字节数

int w = srcImage.width(); //宽

int h = srcImage.height(); //高

byte *pByte;

//遍历读取每个像素,并转换为灰度值

int i, j;

for(i = 0; i < h; i++)

{

for(j = 0; j < w; j++)

{

if(8 == depth) //采用了256色调色板,8位颜色索引

{

pByte = pImageBytes + i * lineBytes + j;

data[i * w + j] = complex<double>( *pByte, 0);

}

else if(32 == depth)//32位表示,数据格式为0xFFBBGGRR或0xAABBGGRR

{

pByte = pImageBytes + i * lineBytes + j * 4;

//根据RGB模式转化成YIQ色彩模式的方式,取Y作为灰 度值

byte pixelValue = (byte)(0.299 * (float)pByte[0] + 0.587 * (float)pByte[1]

+ 0.114 * (float)pByte[2]);

data[i * w + j] = complex<double>( pixelValue, 0);

}

else

{

cout << "invalid format. depth = " << depth << "/n";

return;

}

}

}

}

1.1.2 将数据写入图像

//coef为比例系数,主要用来调整灰度值以便于观察

void ImageProcess::writeImage(QImage &destImage, const complex<double> data[], double coef)

{

int lineBytes = destImage.bytesPerLine();

int depth = destImage.depth();

int w = destImage.width();

int h = destImage.height();

byte *pImageBytes = destImage.bits();

byte *pByte;

for(int i = 0; i < h; i++)

{

for(int j = 0; j < w; j++)

{

double spectral = abs(data[i * w + j]) * coef; //灰度值调整

spectral = spectral > 255 ? 255 : spectral;

//根据图像格式写数据

if(8 == depth)

{

pByte = pImageBytes + i * lineBytes + j;

*pByte = spectral;

}

else if(32 == depth)

{

pByte = pImageBytes + i * lineBytes + j * 4;

pByte[0] = pByte[1] = pByte[2] = spectral;

}

else

{

return;

}

}

}

}

1.1.3 递归形式的快速傅立叶变换

// 数组a为输入,数组y为输出,2的power次方为数组的长度

void ImageProcess::fft(const complex<double> a[], complex<double> y[], int power)

{

if(0 == power)

{

y[0] = a[0];

return;

}

int n = 1 << power;

double angle = 2 * PI / n;

complex<double> wn(cos(angle), sin(angle));

complex<double> w(1, 0);

complex<double> *a0 = new complex<double>[n / 2];

complex<double> *a1 = new complex<double>[n / 2];

complex<double> *y0 = new complex<double>[n / 2];

complex<double> *y1 = new complex<double>[n / 2];

for(int i = 0; i < n / 2; i ++)

{

a0[i] = a[2 * i];

a1[i] = a[2 * i + 1];

}

//分开成两个子fft过程

fft(a0, y0, power - 1);

fft(a1, y1, power - 1);

complex<double> u;

for(int k = 0; k < n / 2; k++) //蝶形算法

{

u = w * y1[k];

y[k] = y0[k] + u;

y[k + n / 2] = y0[k] - u;

w = w * wn;

}

delete[] a0;

delete[] a1;

delete[] y0;

delete[] y1;

}

1.1.4 快速傅立叶反变换

//y为输入,a为输出,2的power次方为数组的长度

void ImageProcess::ifft(const complex<double> y[], complex<double> a[], int power)

{

int count = 1 << power;

complex<double> *x = new complex<double>[count];

memcpy(x, y, sizeof(complex<double>) * count);

int i;

for(i = 0; i < count; i++)

{

x[i] = complex<double>(x[i].real(), -x[i].imag()); //共轭复数

}

fft(x, a, power); //调用快速傅立叶变换算法

for(i = 0; i < count; i++)

{

a[i] = complex<double>(a[i].real() / count, -a[i].imag() / count); //共轭复数

}

delete[] x;

}

1.1.5 调整图像的 大小

//宽和高都截取为2的指数倍

void ImageProcess::adjustImageSize(QImage &image)

{

int w = 1;

int h = 1;

int width = image.width();

int height = image.height();

wp = 0, hp = 0;

while(w * 2 <= width) { w *= 2; wp++; }

while(h * 2 <= height) {h *= 2; hp++;}

QImage adjustedImage(w, h, image.depth(), image.numColors(), image.bitOrder());

byte *destBytes = adjustedImage.bits();

byte *srcBytes = image.bits();

int lineBytes = image.bytesPerLine();

int bytesPerPixel = image.depth() / 8; //每个象素的字节数

for(int i = 0; i < h; i++) //拷贝数据

{

memcpy(destBytes + i * w * bytesPerPixel, srcBytes + i * lineBytes,

sizeof(byte) * w * bytesPerPixel);

}

image = adjustedImage; //更新图像

}

1.1.6 傅立 叶变换的主过程

void ImageProcess::fourier()

{

int w = currentImage.width();

int h = currentImage.height();

if(needAdjust) //调整图像的大小为2的幂次以便于快速傅立叶变换

{

adjustImageSize(currentImage); //调整大小

needAdjust = false;

if(currentImageData)

{

delete[] currentImageData;

}

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage); //读取数据

}

else if(NULL == currentImageData)

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage); //读取数据

}

w = currentImage.width(); //更新宽和高

h = currentImage.height();

complex<double> *TD = currentImageData; //当前读取的数据为时域

complex<double> *FD = new complex<double>[w * h]; //申请空间保存变换结果

int i, j;

for(i = 0; i < h; i++) //在x方向上对按行进行快速傅立叶变换

{

fft(&TD[w * i], &FD[w * i], wp);

}

memcpy(TD, FD, sizeof(complex<double>) * w * h);

complex<double> *columnt = new complex<double>[h];

complex<double> *columnf = new complex<double>[h];

for(i = 0; i < w; i++) //调整行列数据,在y方向上按列进行快速傅立叶变换

{

for(j = 0; j < h; j++)

{

columnt[j] = TD[j * w + i];

}

fft(columnt, columnf, hp);

for(j = 0; j < h; j++)

{

FD[j * w + i] = columnf[j];

}

}

delete[] columnt;

delete[] columnf;

writeImage(currentImage, FD, 0.02); //写入数据

delete[] currentImageData;

currentImageData = FD;

pDispLabel->setPixmap(QPixmap(currentImage));

}

1.1.7 傅立叶反变换

傅立叶反变换的思想与傅立叶变化相似,只是时 域和频域互换,然后调用快速傅立叶反变换ifft而不是快速傅立叶变换fft。

1.2. 运行截图

1.2.1 正方形

输入一个256*256的图形,背景为白色,中间有一黑色的正方形,如图1-1所示。经过傅立叶变换后的 结果如图1-2所示(注:没有采用平移到中心的方法)。

clip_image002

图1-1

clip_image004

图1-2

1.2.2 旋 转45

将图1-1旋转45度后的输入如图1-3所示。其傅立叶变换结果如图1-4所示。

clip_image006

图1-3

clip_image008

图1-4

1.2.3 输 入长方形图像

输入图像如图1-5所示。傅立叶变换结果如图1-6所示。

clip_image010

图1-5

clip_image012

图1-6

1.2.4 傅 立叶反变换

对傅立叶变换结果图1-2进行傅立叶反变换,其结果与原图1-1相同,如图1-7所示:

clip_image014

图1-7


图像增强

图像增强是一种很重要的图像处理技术,为了方便人们观察以及机器处理而去处理给定的一幅图像。有很多图像增 强的方法,以下这部分实现了其中的平滑和锐化这两种方法。

2.1. 主要源码

2.1.1 平滑

平滑采用的模板是clip_image016 , 实现如下:

void ImageProcess::smooth()

{

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要重新读取数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//拷贝一份数据便于计算

complex<double> *buffer = new complex<double>[w * h];

memcpy(buffer, currentImageData, sizeof(complex<double>) * w * h);

//根据模板进 行计算

//为了简化编码忽略了图像边界(i =0 or h, j =0 or w),对于整体效果没有影响

int i, j;

for(i = 1; i < h - 1; i++)

{

for(j = 1; j < w - 1; j++)

{

complex<double> k;

k = buffer[(i - 1) * w + j - 1];

k += buffer[(i - 1) * w + j];

k += buffer[(i - 1) * w + j + 1];

k += buffer[i * w + j - 1];

k += buffer[i * w + j];

k += buffer[i * w + j + 1];

k += buffer[(i + 1) * w + j - 1];

k += buffer[(i + 1) * w + j];

k += buffer[(i + 1) * w + j + 1];

k = complex<double>(k.real() / 9, 0);

currentImageData[i * w + j] = k;

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

2.1.2 锐化

采用拉普拉斯锐化,其模板为clip_image018 ,其实现如下:

void ImageProcess::sharp()

{

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要读取数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//拷贝一份数据便于计算

complex<double> *buffer = new complex<double>[w * h];

memcpy(buffer, currentImageData, sizeof(complex<double>) * w * h);

//根据模板进 行计算

//为了简化编码忽略了图像边界(i =0 or h, j =0 or w),对于整体效果没有影响

int i, j;

complex<double> k;

for(i = 1; i < h - 1; i++)

{

for(j = 1; j < w - 1; j++)

{

k = buffer[i * w + j];

k = complex<double>(k.real() * 5, 0);

k -= buffer[(i - 1) * w + j];

k -= buffer[i * w + j - 1];

k -= buffer[i * w + j + 1];

k -= buffer[(i + 1) * w + j];

currentImageData[i * w + j] = k;

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

2.2. 运行截图

输入图像2-1,其平滑结果为图2-2,锐化结果为 图2-3。

clip_image020

图2-1原来的图像

clip_image022

图2-2 平滑后的图像

clip_image024

图2-3 锐化后的图像


图像分析

这部分主要实现了图像的模板匹配。模板匹配是一种非常原始的模式识别方法。有很多模板匹配的算法。这里采用 的算法是计算二者之间的相似度,在目标图像中选取一个坐标,将以该坐标为左上角选定一块区域,计算该区域与模板的相似度,相似度最大的点即为匹配之处。通 过二者之间的差异度来判断其相似程度,差异度的计算:m = clip_image026 。即将累加其像素之间的差值,为了提高计算速度,可以设置阀值,当m大 于阀值时,认定该块区域不匹配,继续寻找下一区域。

3.1. 主要源码

void ImageProcess::match()

{

//让用户选取模板

QString fileName = QFileDialog::getOpenFileName("/home/tanqiyu", "Images (*.png *.xpm

.jpg)", this, "open file dialog", "Choose a model image");

if(QString::null == fileName)

{

return;

}

//读取模板数据

QImage modelImage(fileName);

int mw = modelImage.width();

int mh = modelImage.height();

complex<double> *modelImageData = new complex<double>[mw * mh];

readImage(modelImageData, modelImage);

unsigned long t = mw * mh * 8; //根据匹配模板的大小设置一定的阀值

unsigned long m = t; //初始差异度

int ri = -1; //z左上角坐标(ri, rj)

int rj = -1;

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要读取目标图像数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//遍历目标图像,选取左上角坐标,考虑到模板图像的大小注意不要越界

int i, j;

for(i = 0; i < h - mh + 1; i++ )

{

for(j = 0; j < w - mw + 1; j++)

{

//下面开始对点(i, j)为左上角的mw * mh区域进行匹配

bool overFlag = false;

unsigned long k = 0; //差异值的累加和

int u, v;

for(u = 0; u < mh && !overFlag; u++)

{

for(v = 0; v < mw && !overFlag; v++)

{

k += abs(currentImageData[(i + u) * w + j + v].real()

- modelImageData[u * mw + v].real()); //计算差值并累加

if(k >= t) //判断是否大于阀值

{

overFlag = true;

}

}

}

if(k < m) //判断是否找到更加匹配的区域

{

ri = i;

rj = j;

m = k;

}

}

}

//找到匹配区域,则将目标图像匹配区域之 外的点置成白色以便于观察结果

if(ri != -1)

{

for(i = 0; i < h; i++)

{

for(j = 0; j < w; j++)

{

if(i < ri || j < rj || i > ri + mh - 1 || j > rj + mw - 1)

{

currentImageData[i * w + j] = complex<double>(255, 0);

}

}

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

3.2. 运行截图

目标图像为图3-1,图3-2位匹配模板,匹配结果为图3-3。

clip_image027

图3-1 目标图像

clip_image028

图3-2 匹配模板

clip_image030

匹配结果

2016-01-09 21:35:13 qq_27789527 阅读数 3964
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4125 人正在学习 去看看 贾志刚

摘要

本文主要总结了进行目标跟踪、检测中经常使用到的图像相似度测量和模板匹配方法,并给出了具体的基于OpenCV的代码实现。

引言

模板匹配是一种在源图像中寻找与图像patch最相似的技术,常常用来进行目标的识别、跟踪与检测。其中最相似肯定是基于某种相似度准则来讲的,也就是需要进行相似度的测量。另外,寻找就需要在图像上进行逐行、逐列的patch窗口扫描,当然也不一定需要逐行逐列的扫描,当几个像素的误差比计算速度来的不重要时就可以设置扫描的行列步进值,以加快扫描和计算的时间消耗。下面就对相似度测量和模板匹配进行介绍(所有的图像都假定是灰度图)。

正文

图像相似度测量

直方图方法

方法描述:有两幅图像patch(当然也可是整幅图像),分别计算两幅图像的直方图,并将直方图进行归一化,然后按照某种距离度量的标准进行相似度的测量。

方法的思想:基于简单的向量相似度来对图像相似度进行度量。

优点:直方图能够很好的归一化,比如256个bin条,那么即使是不同分辨率的图像都可以直接通过其直方图来计算相似度,计算量适中。比较适合描述难以自动分割的图像。

缺点:直方图反应的是图像灰度值得概率分布,并没有图像的空间位置信息在里面,因此,常常出现误判;从信息论来讲,通过直方图转换,信息丢失量较大,因此单一的通过直方图进行匹配显得有点力不从心。

基于opencv的实现,我把它封装为函数(输入输出都是灰度图像,且设定的灰度级为8,及0-255),有需要的可以直接拿去使用

<span style="font-size:18px;"><span style="font-size:18px;">double getHistSimilarity(const Mat& I1, const Mat& I2)
{
	int histSize = 256;
	float range[] = {0,256};
	const float* histRange = {range};
	bool uniform = true;
	bool accumulate = false;

	Mat hist1,hist2;

	calcHist(&I1,1,0,Mat(),hist1,1,&histSize,&histRange,uniform,accumulate);
	normalize(hist1,hist1,0,1,NORM_MINMAX,-1,Mat());

	calcHist(&I2,1,0,Mat(),hist2,1,&histSize,&histRange,uniform,accumulate);
	normalize(hist2,hist2,0,1,NORM_MINMAX,-1,Mat());

	return compareHist(hist1, hist2, CV_COMP_CORREL);

}</span></span>
其中直方图比较的方法在opencv中的实现的有:

/* Histogram comparison methods */
enum
{
    CV_COMP_CORREL        =0,
    CV_COMP_CHISQR        =1,
    CV_COMP_INTERSECT     =2,
    CV_COMP_BHATTACHARYYA =3,
    CV_COMP_HELLINGER     =CV_COMP_BHATTACHARYYA
};
分别是:相关度,卡方,相交系数和巴氏距离,其计算公式如下图所示:

技术分享

直方图方法的一些改进的出发点大致是为了改善直方图无法利用空间位置信息的弱点,比如FragTrack算法,将图像再进行分割成横纵的patch,然后对于每一个patch搜索与之匹配的直方图,来计算两幅图像的相似度,从而融入了直方图对应的位置信息,但计算效率肯定不高。


矩阵分解的方法

方法描述:将图像patch做矩阵分解,比如SVD奇异值分解和NMF非负矩阵分解等,然后再做相似度的计算。

方法思想:因为图像本身来讲就是一个矩阵,可以依靠矩阵分解获取一些更加鲁棒的特征来对图像进行相似度的计算。

基于SVD分解的方法优点:奇异值的稳定性,比例不变性,旋转不变性和压缩性。即奇异值分解是基于整体的表示,不但具有正交变换、旋转、位移、镜像映射等代数和几何上的不变性,而且具有良好的稳定性和抗噪性,广泛应用于模式识别与图像分析中。对图像进行奇异值分解的目的是得到唯一、稳定的特征描述,降低特征空间的维度,提高抗干扰能力。

基于SVD分解的方法缺点是:奇异值分解得到的奇异矢量中有负数存在,不能很好的解释其物理意义。


基于NMF分解的方法:将非负矩阵分解为可以体现图像主要信息的基矩阵与系数矩阵,并且可以对基矩阵赋予很好的解释,比如对人脸的分割,得到的基向量就是人的“眼睛”、“鼻子”等主要概念特征,源图像表示为基矩阵的加权组合,所以,NMF在人脸识别场合发挥着巨大的作用。
基于矩阵特征值计算的方法还有很多,比如Trace变换,不变矩计算等。

基于特征点方法

方法描述:统计两个图像patch中匹配的特征点数,如果相似的特征点数比例最大,则认为最相似,最匹配
方法思想:图像可以中特征点来描述,比如sift特征点,LK光流法中的角点等等。这样相似度的测量就转变为特征点的匹配了。
以前做过一些实验,关于特征点匹配的,对一幅图像进行仿射变换,然后匹配两者之间的特征点,选取的特征点有sift和快速的sift变形版本surf等。
方法优点:能被选作特征点的大致要满足不变性,尺度不变性,旋转不变等。这样图像的相似度计算也就具备了这些不变性。
方法缺点:特征点的匹配计算速度比较慢,同时特征点也有可能出现错误匹配的现象。

基于峰值信噪比(PSNR)的方法

当我们想检查压缩视频带来的细微差异的时候,就需要构建一个能够逐帧比较差视频差异的系统。最
常用的比较算法是PSNR( Peak signal-to-noise ratio)。这是个使用“局部均值误差”来判断差异的最简单的方法,假设有这两幅图像:I1和I2,它们的行列数分别是i,j,有c个通道。每个像素的每个通道的值占用一个字节,值域[0,255]。注意当两幅图像的相同的话,MSE的值会变成0。这样会导致PSNR的公式会除以0而变得没有意义。所以我们需要单独的处理这样的特殊情况。此外由于像素的动态范围很广,在处理时会使用对数变换来缩小范围。
技术分享

基于OpenCV的PSNR方法实现,同样将其封装为函数,可以直接拿来调用:

<span style="font-size:18px;">double getPSNR(const Mat& I1, const Mat& I2)
{
	Mat s1;
	absdiff(I1, I2, s1);       // |I1 - I2|
	s1.convertTo(s1, CV_32F);  // cannot make a square on 8 bits
	s1 = s1.mul(s1);           // |I1 - I2|^2

	Scalar s = sum(s1);         // sum elements per channel

	double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels

	if( sse <= 1e-10) // for small values return zero
		return 0;
	else
	{
		double  mse =sse /(double)(I1.channels() * I1.total());
		double psnr = 10.0*log10((255*255)/mse);
		return psnr;
	}
}</span>

考察压缩后的视频时,这个值大约在30到50之间,数字越大则表明压缩质量越好。如果图像差异很明显,就可能会得到15甚至更低的值。PSNR算法简单,检查的速度也很快。但是其呈现的差异值有时候和人的主观感受不成比例。所以有另外一种称作 结构相似性 的算法做出了这方面的改进。

基于结构相似性(SSIM,structural similarity (SSIM) index measurement)的方法


结构相似性理论认为,自然图像信号是高度结构化的,即像素间有很强的相关性,特别是空域中最接近的像素,这种相关性蕴含着视觉场景中物体结构的重要信息;HVS的主要功能是从视野中提取结构信息,可以用对结构信息的度量作为图像感知质量的近似。结构相似性理论是一种不同于以往模拟HVS低阶的组成结构的全新思想,与基于HVS特性的方法相比,最大的区别是自顶向下与自底向上的区别。这一新思想的关键是从对感知误差度量到对感知结构失真度量的转变。它没有试图通过累加与心理物理学简单认知模式有关的误差来估计图像质量,而是直接估计两个复杂结构信号的结构改变,从而在某种程度上绕开了自然图像内容复杂性及多通道去相关的问题.作为结构相似性理论的实现,结构相似度指数从图像组成的角度将结构信息定义为独立于亮度、对比度的,反映场景中物体结构的属性,并将失真建模为亮度、对比度和结构三个不同因素的组合。用均值作为亮度的估计,标准差作为对比度的估计,协方差作为结构相似程度的度量。

基于OpenCV的SSIM方法实现,同样将其封装为函数,可以直接拿来调用:

<span style="font-size:18px;">Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
	const double C1 = 6.5025, C2 = 58.5225;
	/***************************** INITS **********************************/
	int d     = CV_32F;

	Mat I1, I2;
	i1.convertTo(I1, d);           // cannot calculate on one byte large values
	i2.convertTo(I2, d);

	Mat I2_2   = I2.mul(I2);        // I2^2
	Mat I1_2   = I1.mul(I1);        // I1^2
	Mat I1_I2  = I1.mul(I2);        // I1 * I2

	/*************************** END INITS **********************************/

	Mat mu1, mu2;   // PRELIMINARY COMPUTING
	GaussianBlur(I1, mu1, Size(11, 11), 1.5);
	GaussianBlur(I2, mu2, Size(11, 11), 1.5);

	Mat mu1_2   =   mu1.mul(mu1);
	Mat mu2_2   =   mu2.mul(mu2);
	Mat mu1_mu2 =   mu1.mul(mu2);

	Mat sigma1_2, sigma2_2, sigma12;

	GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
	sigma1_2 -= mu1_2;

	GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
	sigma2_2 -= mu2_2;

	GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
	sigma12 -= mu1_mu2;

	///////////////////////////////// FORMULA ////////////////////////////////
	Mat t1, t2, t3;

	t1 = 2 * mu1_mu2 + C1;
	t2 = 2 * sigma12 + C2;
	t3 = t1.mul(t2);              // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))

	t1 = mu1_2 + mu2_2 + C1;
	t2 = sigma1_2 + sigma2_2 + C2;
	t1 = t1.mul(t2);               // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))

	Mat ssim_map;
	divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;

	Scalar mssim = mean( ssim_map ); // mssim = average of ssim map
	return mssim;
}
</span>

图像模板匹配

一般而言,源图像与模板图像patch尺寸一样的话,可以直接使用上面介绍的图像相似度测量的方法;如果源图像与模板图像尺寸不一样,通常需要进行滑动匹配窗口,扫面个整幅图像获得最好的匹配patch。

在OpenCV中对应的函数为:matchTemplate():函数功能是在输入图像中滑动窗口寻找各个位置与模板图像patch的相似度。相似度的评价标准(匹配方法)有:CV_TM_SQDIFF平方差匹配法(相似度越高,值越小),CV_TM_CCORR相关匹配法(采用乘法操作,相似度越高值越大),CV_TM_CCOEFF相关系数匹配法(1表示最好的匹配,-1表示最差的匹配)。

CV_TM_SQDIFF计算公式:

技术分享

CV_TM_CCORR计算公式:

技术分享


有一种新的用来计算相似度或者进行距离度量的方法:EMD,Earth Mover‘s Distances

EMD is defined as the minimal cost that must be paid to transform one histograminto the other, where there is a “ground distance” between the basic featuresthat are aggregated into the histogram。

光线变化能引起图像颜色值的漂移,尽管漂移没有改变颜色直方图的形状,但漂移引起了颜色值位置的变化,从而可能导致匹配策略失效。而EMD是一种度量准则,度量怎样将一个直方图转变为另一个直方图的形状,包括移动直方图的部分(或全部)到一个新的位置,可以在任意维度的直方图上进行这种度量。

在OpenCV中有相应的计算方法:cvCalcEMD2()。结合着opencv支持库,计算直方图均衡后与原图的HSV颜色空间直方图之间的EMD距离。
*******************************************************************************************************************************************************************************************************

2015-7-24

版权声明:本文为博主原创文章,未经博主允许不得转载。

图像相似度测量与模板匹配总结

标签:图像相似度   图像匹配   直方图匹配   特征点匹配   图像的矩阵分解   

2011-02-19 07:56:00 liminlu0314 阅读数 4524
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4125 人正在学习 去看看 贾志刚

     上文说到使用OpenCV进行模板匹配的函数matchTemplate,下面就matchTemplate函数的内部处理过程做一个简单的说明。matchTemplate函数的源代码在OpenCV的源代码目录下的 modules/imgproc/src/templmatch.cpp 文件中。其核心函数代码如下(其中的注释是我添加的):

void matchTemplate( const Mat& _img, const Mat& _templ, Mat& result, int method )
{
    CV_Assert( CV_TM_SQDIFF <= method && method <= CV_TM_CCOEFF_NORMED );
    //numType用来表示模板匹配的方式,0表示相关匹配法,1表示相关系数匹配法,2表示平方差匹配法
    //isNormed表示是否进行归一化处理,true表示进行归一化,false表示不进行归一化处理
    int numType = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 :
                  method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2;
    bool isNormed = method == CV_TM_CCORR_NORMED ||
                    method == CV_TM_SQDIFF_NORMED ||
                    method == CV_TM_CCOEFF_NORMED;
    //判断两幅图像的大小关系,如果输入的原始图像比匹配图像要小,则将原始图像作为模板,原来的模板图像作为搜索图
    Mat img = _img, templ = _templ;
    if( img.rows < templ.rows || img.cols < templ.cols )
        std::swap(img, templ);
    
    CV_Assert( (img.depth() == CV_8U || img.depth() == CV_32F) &&
               img.type() == templ.type() );

   //crossCorr函数是将输入图像做了一次DFT变换(离散傅里叶变换),将空间域的图像转换到频率域中来进行处理,并将处理的结果存放在result中
    int cn = img.channels();
    crossCorr( img, templ, result,
               Size(img.cols - templ.cols + 1, img.rows - templ.rows + 1),
               CV_32F, Point(0,0), 0, 0);

    //如果是相关匹配方法,此处已经计算完毕,返回
if( method == CV_TM_CCORR ) return;

    //将模板看作单位1,计算每一个像元所占的百分比(也可以理解为整个模板面积为1,计算每个像元的面积)
double invArea = 1./((double)templ.rows * templ.cols); Mat sum, sqsum; Scalar templMean, templSdv; double *q0 = 0, *q1 = 0, *q2 = 0, *q3 = 0; double templNorm = 0, templSum2 = 0;

    //相关系数匹配算法
if( method == CV_TM_CCOEFF ) { integral(img, sum, CV_64F);//对原始图像进行求和 templMean = mean(templ);//计算模板图像的均值向量 } else//其他匹配算法 { integral(img, sum, sqsum, CV_64F);//计算原始图像的和以及平方和 meanStdDev( templ, templMean, templSdv );//计算模板图像的均值向量和方差向量 templNorm = CV_SQR(templSdv[0]) + CV_SQR(templSdv[1]) + CV_SQR(templSdv[2]) + CV_SQR(templSdv[3]);//计算所有通道的方差和 if( templNorm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED ) {//如果所有通道的方差的和等于0,并且使用的方法是归一化相关系数匹配方法,则返回 result = Scalar::all(1); return; } templSum2 = templNorm + CV_SQR(templMean[0]) + CV_SQR(templMean[1]) + CV_SQR(templMean[2]) + CV_SQR(templMean[3]);//计算所有通道的均值的平方和 if( numType != 1 )//匹配方式不是相关系数,对模板均值向量和templNorm重新赋值 { templMean = Scalar::all(0); templNorm = templSum2; } templSum2 /= invArea; templNorm = sqrt(templNorm); templNorm /= sqrt(invArea); // care of accuracy here q0 = (double*)sqsum.data; q1 = q0 + templ.cols*cn; q2 = (double*)(sqsum.data + templ.rows*sqsum.step); q3 = q2 + templ.cols*cn; }

    //下面就是在结果图像中进行查找匹配的结果位置,代码略去,具体可参考OpenCV源代码

 

    在测试过程中,发现直接用OpenCV打开超过1G(可能比这个大小还要小)的图像就会导致错误,查看OpenCV代码,发现其一次将所有的图像数据全部载入内存,导致内存不够出错,而且OpenCV只支持普通的常用的图像数据格式,并不能直接打开Erdas的img格式或者PCI的pix格式,针对这两个问题,我提出的解决思路是实用GDAL来作为图像数据的读取库,然后分块读取一部分数据(指定一个大小,我指定的是每次读64M的数据),如果图像数据超过该大小,则分块处理,每次处理一块,OpenCV可以建立一个内存数据,就是将图像的RGB值按照特定的顺序读取到内存中,然后将该指针交给OpenCV,OpenCV可以通过该内存数据来构建一个图像,然后再调用模板匹配算法,计算出的结果在还原到原来大的图像中的位置。

    在实际的执行过程中,上述方式是可行的,执行效率大概在一个1G的img图像查找一个128*128的模板需要1分钟左右,相比按照前一篇博客中的方式速度提升了很多倍(1G的数据大概计算了4个小时还没有算完)。对于1min中查找出结果应该是可以接受的,但是感觉还是有点慢,于是想到首先将模板和原始图像同时按照一定的采样比例建立金字塔,然后从金字塔的顶层开始进行匹配,根据顶层找到的位置,再到下一层对应的区域中查找,依次知道查找到原始图像级别。这样就可以减少很多的运算。经过测试,使用这种方式,还是上面的1G的数据,查找只要10S钟即可完成,而且查找的结果也是很准确的。

    在编写程序的过程中,可以使用GDAL的RasteIO函数,直接使用最后面的三个参数,可以直接将图像中的数据读取成OpenCV中的图像存储方式 BGRBGRBGR…

   对于OpenCV在进行第一步匹配的时候,是在频率域中进行处理的,相关的函数也在OpenCV的源代码目录下的 modules/imgproc/src/templmatch.cpp 文件中,关于DFT变换(离散傅里叶变换)比较复杂,理论知识尤其是需要很好的数学功底才能够理解,感兴趣的童鞋可以google或者参考相关的数字图像处理的书籍,一般的数字图像处理的书籍中都会有的,因为DFT和FFT在数字图像处理以及信号处理方面使用的非常广泛。

Technorati 标签: ,

2018-01-05 19:15:09 m0_37402140 阅读数 17301
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4125 人正在学习 去看看 贾志刚


源代码及毕业论文参考模板下载

指纹识别的一般步骤为指纹采集、预处理、特征点提取、特征点匹配。指纹分为螺旋形、弓形、环形。指纹的处理效果影响着后面特征点的提取和识别效果,所以图像的预处理占有重要的地位。指纹采集一般有专业的设备,所以这一步骤一搬不关注。

①预处理

因为采集指纹时力度和各种因素所以采集的指纹灰度图会有很大不同,首先对图像进行归一化处理,归一化主要针对两个步骤:大小和灰度值。把采集到的指纹图统一调整到特定大小。灰度值我会根据整幅图的均值方差调整到某一范围内。

归一化处理完毕后会对图像进行分割处理,目的是区分出前景色和背景色。我采用的分割为根据多区域阈值分割。多区域分割的效果取决于区域的大小,而指纹的区域分为一脊一谷最好,所以我选择3x3的区域大小。我会根据对区域多次进行求均值和方差进行分割。采集到的指纹图背景的灰度值大于前景色,背景主要为低频,所以背景的方差小于前景的方差。我分别求得背景和前景的均值和方差然后会得到背景为白色 脊线为黑色。然后保存在矩阵e(二值图)中,我会根据e中位置等于1的点的八邻域点的和小于四得到背景色,达到背景和前景分离(e矩阵)。然后黑白反转让感兴趣的前景色变为白色(保存在Icc中),灰度图(gray)的背景值替换为小区域块的和的均值(G1.但是得到的脊线方向并不能达到准确识别指纹。所以下一步会沿脊线方向增强指纹纹路,采用的方法为基于脊线方向场的增强方法。为了估计脊线的方向场,把脊线的方向场划分为八个方向,然后根据八个方向的灰度值的总和来得到脊线的方向。并对图像进行二值化。此时脊线还是为黑色。因为各种采集原因(油脂水分等)会使指纹粘连断裂,会影响后续的特征提取和识别,接下来会去除指纹中的空洞和毛刺,如果当前位置点值为0(背景)该点的四邻域点(上下左右)的和大于3则为毛刺,空洞的判断方法为该点为白色(背景)的四周为黑色(前景)八领域点两的和为0,则为空洞。我们得到的图像的纹线仍具有一定的宽度,而指纹的识别只与纹线的走向有关。所以我们只需要纹线的宽度为一个像素宽度即可。下面我执行了黑白反转使感兴趣的区域(纹线)变为白色。在执行开操作和闭操作使边界平滑,消除细小的尖刺,断开窄小的连接。执行细化(bwmorph)得到细画图(thin)。

②特征点提取

特征点提取的点为端点和交叉点,遍历细化图的每一个像素点,端点的判别方法为八领域点两两相减取绝对值求和如果值为2则为端点(周围只有一个为1的白色点)和为6时为交叉点(周围有三值为1的白色点)。把找到的端点位置存储在txyNx3,第一列存储x的位置,第二列存储y的位置,第三列标识点的类别,如果为交叉点则值为6,如果端点则值为2.然后对纹线再次进行光滑处理。原理:找到每个端点,使其沿着纹线的方向移动五个像素,如果在五个像素之内遇到交叉点则此点为毛刺,去除此点。判断离端点num个距离内是否有另一个端点的函数为walk判断的主要方法先将该点置为0然后根据八领域点和;和为0或大于2则证明该点距离num内有端点,然后循环向前遍历。惊醒光滑处理后的特征点比以前少了。但是我们采集指纹时由于采集器的关系,图像的边缘会有很多端点,这已然会影响后续的识别工作,所以我们需要熊特征点中去除这些端点,cut函数主要做这些工作,边缘的主要特征就是黑色多白色少,也就是均值小,所以我主要里用灰度图的分区域(31*31)求均值如果灰度值小于70则在该区域的特征点去除,然后会得到更少的特征点,但这些点依然不够少说明不够特殊,下面定义了combinethinrtxynum)函数可以找出周围半径为r个像素的圆内没有任何端点和交叉点,沿纹线走num个距离没有交叉点和端点,cominewalksingle_point函数的综合,walk函数上面已经介绍,single_point函数主要是找出独特的点作为特征点,原理是根据两个端点之间的距离。求每个端点距离其他端点的距离,找取距离大于r的端点。执行完combine后会得到更少的端点(实验结果为3个端点)。

③特征点匹配

特征点匹配主要采用三个方法

1. 根据距离判断

找到某一个特征点,从该特征点沿着纹线走num个距离,并计算出每走一步距离该特征点的距离,最后会得到哟个装有长度信息的数组,如果两幅指纹相同则他们含有相同的特征点而且得到的数组对应的位置的数据基本相等

2三角形边长匹配,找到一个特征点以后,可以找出距离最近的两个端点与原特征点构成三角形,若两幅图的三角形的边长比例相等则说明两幅图匹配

3. 点类型匹配

4. 找到一个特征点以后,找出距离最近的num个端点,统计num个端点中端点和交叉点的个数,若两幅图匹配,则端点占的比例大致相同

读入指纹图片,通过预处理,特征点提取来识别指纹

二、现实意义:

指纹识别技术在现实生活中的意义主要体验在隐私安全保护上,可通过指纹特征点的多个特征来识别不同人的指纹来达到保护隐私安全,是具有及其重要的现实作用的。

三、涉及知识内容:

1、中值滤波

2、开运算

3、闭运算

4、二值化

5、图像细化

四、实例分析及截图效果:

(1)代码显示:

 

 

1、程序中定义图像重要变量说明

 

1image--------------------------------------------------------------原图变量;

2gray-------------------------------------------------------归一化后的图像;

3Icc-------------------------背景前景分离后的二值图;

4thin---细化图

5txy--------------------------特征点存储的矩阵;

6pxy----------------------特征点筛选后的特征点;

2重要实现代码:

%找出独特的端点作为特征点  一个端点的周围半径为r个像素的圆内没有任何端点或交叉点

%随着r的逐渐变大,这样的点会变少,也就会越来越独特

function [pxy2,error]=single_point(txy,r)
error=0;
x=txy(:,1);
y=txy(:,2);
n=length(x);
d(1:n,1:n)=0;
for j=1:n
    for i=1:n
        if(i~=j)
            %计算两点间距离
            d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
        else
            d(i,j)=2*r;
        end
    end
end
[a,b]=min(d);
c=find(a>r);
pxy2=txy(c,:);
pxy2=pxy2(find(pxy2(:,3)==2),:);
t=size(pxy2,1);
if t==0
    error=1;
else
    plot(x,y,'b.');
    hold on
    plot(pxy2(:,1),pxy2(:,2),'r.');
End

 

 

function [pxy2,error]=single_point(txy,r)
error=0;
x=txy(:,1);
y=txy(:,2);
n=length(x);
d(1:n,1:n)=0;
for j=1:n
    for i=1:n
        if(i~=j)
            %计算两点间距离
            d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
        else
            d(i,j)=2*r;
        end
    end
end
[a,b]=min(d);
c=find(a>r);
pxy2=txy(c,:);
pxy2=pxy2(find(pxy2(:,3)==2),:);
t=size(pxy2,1);
if t==0
    error=1;
else
    plot(x,y,'b.');
    hold on
    plot(pxy2(:,1),pxy2(:,2),'r.');
End

 

 

 

 

3、运行效果截图:

 

第一步:读取原图,并显示

f=imread(image);
f=imresize(f,[363,312]);
figure;imshow(f);

 

第二步:归一化,将指纹灰度值调整到特定范围内

%归一化,灰度值限制在某一范围
M=0;var=0;
%均值
for x=1:m
    for y=1:n
        M=M+gray(x,y);
    end
end
M1=M/(m*n);
%方差
for x=1:m
    for y=1:n
        var=var+(gray(x,y)-M1).^2;
    end;
end;
var1=var/(m*n);
for x=1:m
    for y=1:n
        if gray(x,y)>M1
            gray(x,y)=150+sqrt(2000*(gray(x,y)-M1)/var1);
        else 
            gray(x,y)=150-sqrt(2000*(M1-gray(x,y))/var1);
        end
    end
end
figure;imshow(uint8(gray));



 

 

第三步:细化图

Theshold = graythresh(Image);%取得图象的全局域值
Image_BW = im2bw(Image,Theshold);%二值化图象
figure,imshow(Image_BW);
title(' 【初次二值化图像】');



 

四步特征点提取

%找出独特的端点作为特征点  一个端点的周围半径为r个像素的圆内没有任何端点或交叉点
%随着r的逐渐变大,这样的点会变少,也就会越来越独特
function [pxy2,error]=single_point(txy,r)
error=0;
x=txy(:,1);
y=txy(:,2);
n=length(x);
d(1:n,1:n)=0;
for j=1:n
    for i=1:n
        if(i~=j)
            %计算两点间距离
            d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
        else
            d(i,j)=2*r;
        end
    end
end
[a,b]=min(d);
c=find(a>r);
pxy2=txy(c,:);
pxy2=pxy2(find(pxy2(:,3)==2),:);
t=size(pxy2,1);
if t==0
    error=1;
else
    plot(x,y,'b.');
    hold on
    plot(pxy2(:,1),pxy2(:,2),'r.');
End

 

 

步:特征点的再次提取

% 综合walk和single_point函数,通过执行[pxy3,error2]=combine(thin,r,txy,num)可以找出周围半径为r个像素的院内没有任何交叉点或端点,并且沿纹线走num个
% 距离内没有任何哟个交叉点或端点
function [pxy3,error2]=combine(thin,r,txy,num)
error=0;
[pxy2,error]=single_point(txy,r);
n=size(pxy2,1);
k=1;
erroe2=0;
for i=1:n
    [error,a,b]=walk(thin,pxy2(i,1),pxy2(i,2),num);
    if error~=1
        pxy3(k,1)=pxy2(i,1);
        pxy3(k,2)=pxy2(i,2);
        pxy3(k,3)=pxy2(i,3);
        k=k+1;
        error2=0;
        plot(pxy2(i,1),pxy2(i,2),'r+');
    end
end

 

特征识别

% 特征点匹配
% 三角形边长匹配
% 找到一个特征点后,可以找出距离其最近的两个端点,与原特征点构成三角形,瑞两幅图像的三角形边长比例相等,则说明匹配
% find_point()找到最近的端点
function pxy=find_point(x0,y0,txy,num)
x=txy(:,1);
y=txy(:,2);
n=length(x);
k(1,n)=0;
lnn=0;
pxy(num,:)=[0,0,0];
for i=1:n
    k(i)=sqrt((x(i)-x0)^2+(y(i)-y0)^2);
end
kk=sort(k);
for i=1:num
    xiao=kk(i+lnn);
    nn=find(k==xiao);
    lnn=length(nn);
    pxy(i,:)=[x(nn(1)),y(nn(1)),txy(nn(1),3)];
end
plot(x0,y0,'bo');
x0;
y0;
hold on;
plot(pxy(:,1),pxy(:,2),'ro');
 
% ff=(sum(abs((dd1./dd2)-1))) ff越接近0 匹配度越高

 

 

 

% 特征点匹配
% 纹线长度匹配 对于找到的特征点和纹线 沿着纹线走5个像素到原始端点的距离
function d=distance(x0,y0,num,thin)
num2=fix(num/5);
for i=1:num2
    [error,a,b]=walk(thin,x0,y0,5*i);
    if error~=1
        d(i)=sqrt((a-x0)^2+(b-y0)^2);
    else
        break;
    end
end
 
% 最后会得到一个装有长度信息的数组,如果两幅指纹途中的指纹是一样的,则他们含有相同的特征点和从这个特征点出发画出的纹线
% 则这两个数组对应位置的数据基本相等(图像大小相同,则他们的比例接近1)
% f=(sum(abs((d1./d2)-1)))f越接近0 匹配度越高

 

 ③

% 找到一个特征点以后,在其周围找到四十个端点或交叉点,统计四十个点的交叉点和端点的个数
% 若两幅图断点占的比例近似相同则匹配
% fff=abs(f11-f21)/(f11+f12) 越接近于0 匹配度越高
close all;
tic;
clear;
thin1=tuxiangyuchuli('zhiwen.png');
thin2=tuxiangyuchuli('zhiwen.png');
figure;
txy1=point(thin1);
txy2=point(thin2);
[w1,txy1]=guanghua(thin2,txy2);
[w2,txy2]=guanghua(thin2,txy2);
thin1=w1;
thin2=w2;
txy1=cut(thin1,txy1);
txy2=cut(thin2,txy2);
[pxy31,error2]=combine(thin1,8,txy1,60);
[pxy32,error2]=combine(thin2,8,txy2,60);
error=1;
num=20;
cxy1=pxy31;
cxy2=pxy32;
d1=distance(cxy1(1,1),cxy1(1,2),num,thin1);
d2=distance(cxy2(1,1),cxy2(1,2),num,thin2);
f=(sum(abs((d1./d2)-1)));
if(f<0.5)
    error=0;
else
    error=1;
end

 

 

步:最终检测结果:

 

f =

 

     0

 

 

ff =

 

     0

 

 

error =

 

     0

 

 

fff =

 

     0

 

时间已过 13.628304 秒。

四、算法分析

1中值滤波

    利用中值滤波可以对图像进行平滑处理。其算法简单,时间复杂度低,但其对点、线和尖顶多的图像不宜采用中值滤波。很容易自适应化

2开运算
  先腐蚀后膨胀的过程称为开运算。用来消除小物体、在纤细点处分离物体、平滑较大物体的边界的同时并不明显改变其面积。

3闭运算
    先膨胀后腐蚀的过程称为闭运算。用来填充物体内细小空洞、连接邻近物体、平滑其边界的同时并不明显改变其面积。



 



                                
2013-06-28 19:11:56 JoeBlackzqq 阅读数 1975
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4125 人正在学习 去看看 贾志刚

From: http://www.cnblogs.com/qytan36/archive/2010/04/04/1704226.html

实验环境:
1,Linux操作系统
2,QT3编程开发环境
3,C++编程语言

傅立叶变换和傅立叶反变换

1.1. 主要源代码

readImage() 从图像中读取数据
writeImage() 往图像中写入数据
fft() 快速傅立叶变换
ifft() 快速傅立叶反变换
adjustImageSize() 调整图像大小
fourier() 傅立叶变换
ifourier() 傅立叶反变换

1.1.1 从图像中读取数据

void ImageProcess::readImage(complex<double> data[], const QImage &srcImage)

{

byte *pImageBytes = srcImage.bits(); //数据首地址

int depth = srcImage.depth(); //每个像素的bit数

int lineBytes = srcImage.bytesPerLine(); //每行的字节数

int w = srcImage.width(); //宽

int h = srcImage.height(); //高

byte *pByte;

//遍历读取每个像素,并转换为灰度值

int i, j;

for(i = 0; i < h; i++)

{

for(j = 0; j < w; j++)

{

if(8 == depth) //采用了256色调色板,8位颜色索引

{

pByte = pImageBytes + i * lineBytes + j;

data[i * w + j] = complex<double>( *pByte, 0);

}

else if(32 == depth)//32位表示,数据格式为0xFFBBGGRR或0xAABBGGRR

{

pByte = pImageBytes + i * lineBytes + j * 4;

//根据RGB模式转化成YIQ色彩模式的方式,取Y作为灰度值

byte pixelValue = (byte)(0.299 * (float)pByte[0] + 0.587 * (float)pByte[1]

+ 0.114 * (float)pByte[2]);

data[i * w + j] = complex<double>( pixelValue, 0);

}

else

{

cout << "invalid format. depth = " << depth << "\n";

return;

}

}

}

}

1.1.2 将数据写入图像

//coef为比例系数,主要用来调整灰度值以便于观察

void ImageProcess::writeImage(QImage &destImage, const complex<double> data[], double coef)

{

int lineBytes = destImage.bytesPerLine();

int depth = destImage.depth();

int w = destImage.width();

int h = destImage.height();

byte *pImageBytes = destImage.bits();

byte *pByte;

for(int i = 0; i < h; i++)

{

for(int j = 0; j < w; j++)

{

double spectral = abs(data[i * w + j]) * coef; //灰度值调整

spectral = spectral > 255 ? 255 : spectral;

//根据图像格式写数据

if(8 == depth)

{

pByte = pImageBytes + i * lineBytes + j;

*pByte = spectral;

}

else if(32 == depth)

{

pByte = pImageBytes + i * lineBytes + j * 4;

pByte[0] = pByte[1] = pByte[2] = spectral;

}

else

{

return;

}

}

}

}

1.1.3 递归形式的快速傅立叶变换

//数组a为输入,数组y为输出,2的power次方为数组的长度

void ImageProcess::fft(const complex<double> a[], complex<double> y[], int power)

{

if(0 == power)

{

y[0] = a[0];

return;

}

int n = 1 << power;

double angle = 2 * PI / n;

complex<double> wn(cos(angle), sin(angle));

complex<double> w(1, 0);

complex<double> *a0 = new complex<double>[n / 2];

complex<double> *a1 = new complex<double>[n / 2];

complex<double> *y0 = new complex<double>[n / 2];

complex<double> *y1 = new complex<double>[n / 2];

for(int i = 0; i < n / 2; i ++)

{

a0[i] = a[2 * i];

a1[i] = a[2 * i + 1];

}

//分开成两个子fft过程

fft(a0, y0, power - 1);

fft(a1, y1, power - 1);

complex<double> u;

for(int k = 0; k < n / 2; k++) //蝶形算法

{

u = w * y1[k];

y[k] = y0[k] + u;

y[k + n / 2] = y0[k] - u;

w = w * wn;

}

delete[] a0;

delete[] a1;

delete[] y0;

delete[] y1;

}

1.1.4 快速傅立叶反变换

//y为输入,a为输出,2的power次方为数组的长度

void ImageProcess::ifft(const complex<double> y[], complex<double> a[], int power)

{

int count = 1 << power;

complex<double> *x = new complex<double>[count];

memcpy(x, y, sizeof(complex<double>) * count);

int i;

for(i = 0; i < count; i++)

{

x[i] = complex<double>(x[i].real(), -x[i].imag()); //共轭复数

}

fft(x, a, power); //调用快速傅立叶变换算法

for(i = 0; i < count; i++)

{

a[i] = complex<double>(a[i].real() / count, -a[i].imag() / count); //共轭复数

}

delete[] x;

}

1.1.5 调整图像的大小

//宽和高都截取为2的指数倍

void ImageProcess::adjustImageSize(QImage &image)

{

int w = 1;

int h = 1;

int width = image.width();

int height = image.height();

wp = 0, hp = 0;

while(w * 2 <= width) { w *= 2; wp++; }

while(h * 2 <= height) {h *= 2; hp++;}

QImage adjustedImage(w, h, image.depth(), image.numColors(), image.bitOrder());

byte *destBytes = adjustedImage.bits();

byte *srcBytes = image.bits();

int lineBytes = image.bytesPerLine();

int bytesPerPixel = image.depth() / 8; //每个象素的字节数

for(int i = 0; i < h; i++) //拷贝数据

{

memcpy(destBytes + i * w * bytesPerPixel, srcBytes + i * lineBytes,

sizeof(byte) * w * bytesPerPixel);

}

image = adjustedImage; //更新图像

}

1.1.6 傅立叶变换的主过程

void ImageProcess::fourier()

{

int w = currentImage.width();

int h = currentImage.height();

if(needAdjust) //调整图像的大小为2的幂次以便于快速傅立叶变换

{

adjustImageSize(currentImage); //调整大小

needAdjust = false;

if(currentImageData)

{

delete[] currentImageData;

}

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage); //读取数据

}

else if(NULL == currentImageData)

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage); //读取数据

}

w = currentImage.width(); //更新宽和高

h = currentImage.height();

complex<double> *TD = currentImageData; //当前读取的数据为时域

complex<double> *FD = new complex<double>[w * h]; //申请空间保存变换结果

int i, j;

for(i = 0; i < h; i++) //在x方向上对按行进行快速傅立叶变换

{

fft(&TD[w * i], &FD[w * i], wp);

}

memcpy(TD, FD, sizeof(complex<double>) * w * h);

complex<double> *columnt = new complex<double>[h];

complex<double> *columnf = new complex<double>[h];

for(i = 0; i < w; i++) //调整行列数据,在y方向上按列进行快速傅立叶变换

{

for(j = 0; j < h; j++)

{

columnt[j] = TD[j * w + i];

}

fft(columnt, columnf, hp);

for(j = 0; j < h; j++)

{

FD[j * w + i] = columnf[j];

}

}

delete[] columnt;

delete[] columnf;

writeImage(currentImage, FD, 0.02); //写入数据

delete[] currentImageData;

currentImageData = FD;

pDispLabel->setPixmap(QPixmap(currentImage));

}

1.1.7 傅立叶反变换

傅立叶反变换的思想与傅立叶变化相似,只是时域和频域互换,然后调用快速傅立叶反变换ifft而不是快速傅立叶变换fft。

1.2. 运行截图

1.2.1 正方形

输入一个256*256的图形,背景为白色,中间有一黑色的正方形,如图1-1所示。经过傅立叶变换后的结果如图1-2所示(注:没有采用平移到中心的方法)。

clip_image002

图1-1

clip_image004

图1-2

1.2.2 旋转45度

将图1-1旋转45度后的输入如图1-3所示。其傅立叶变换结果如图1-4所示。

clip_image006

图1-3

clip_image008

图1-4

1.2.3 输入长方形图像

输入图像如图1-5所示。傅立叶变换结果如图1-6所示。

clip_image010

图1-5

clip_image012

图1-6

1.2.4 傅立叶反变换

对傅立叶变换结果图1-2进行傅立叶反变换,其结果与原图1-1相同,如图1-7所示:

clip_image014

图1-7


图像增强

图像增强是一种很重要的图像处理技术,为了方便人们观察以及机器处理而去处理给定的一幅图像。有很多图像增强的方法,以下这部分实现了其中的平滑和锐化这两种方法。

2.1. 主要源码

2.1.1 平滑

平滑采用的模板是clip_image016, 实现如下:

void ImageProcess::smooth()

{

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要重新读取数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//拷贝一份数据便于计算

complex<double> *buffer = new complex<double>[w * h];

memcpy(buffer, currentImageData, sizeof(complex<double>) * w * h);

//根据模板进行计算

//为了简化编码忽略了图像边界(i =0 or h, j =0 or w),对于整体效果没有影响

int i, j;

for(i = 1; i < h - 1; i++)

{

for(j = 1; j < w - 1; j++)

{

complex<double> k;

k = buffer[(i - 1) * w + j - 1];

k += buffer[(i - 1) * w + j];

k += buffer[(i - 1) * w + j + 1];

k += buffer[i * w + j - 1];

k += buffer[i * w + j];

k += buffer[i * w + j + 1];

k += buffer[(i + 1) * w + j - 1];

k += buffer[(i + 1) * w + j];

k += buffer[(i + 1) * w + j + 1];

k = complex<double>(k.real() / 9, 0);

currentImageData[i * w + j] = k;

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

2.1.2 锐化

采用拉普拉斯锐化,其模板为clip_image018,其实现如下:

void ImageProcess::sharp()

{

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要读取数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//拷贝一份数据便于计算

complex<double> *buffer = new complex<double>[w * h];

memcpy(buffer, currentImageData, sizeof(complex<double>) * w * h);

//根据模板进行计算

//为了简化编码忽略了图像边界(i =0 or h, j =0 or w),对于整体效果没有影响

int i, j;

complex<double> k;

for(i = 1; i < h - 1; i++)

{

for(j = 1; j < w - 1; j++)

{

k = buffer[i * w + j];

k = complex<double>(k.real() * 5, 0);

k -= buffer[(i - 1) * w + j];

k -= buffer[i * w + j - 1];

k -= buffer[i * w + j + 1];

k -= buffer[(i + 1) * w + j];

currentImageData[i * w + j] = k;

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

2.2. 运行截图

输入图像2-1,其平滑结果为图2-2,锐化结果为 图2-3。

clip_image020

图2-1原来的图像

clip_image022

图2-2 平滑后的图像

clip_image024

图2-3 锐化后的图像


图像分析

这部分主要实现了图像的模板匹配。模板匹配是一种非常原始的模式识别方法。有很多模板匹配的算法。这里采用的算法是计算二者之间的相似度,在目标图像中选取一个坐标,将以该坐标为左上角选定一块区域,计算该区域与模板的相似度,相似度最大的点即为匹配之处。通过二者之间的差异度来判断其相似程度,差异度的计算:m = clip_image026。即将累加其像素之间的差值,为了提高计算速度,可以设置阀值,当m大于阀值时,认定该块区域不匹配,继续寻找下一区域。

3.1. 主要源码

void ImageProcess::match()

{

//让用户选取模板

QString fileName = QFileDialog::getOpenFileName("/home/tanqiyu", "Images (*.png *.xpm

.jpg)", this, "open file dialog", "Choose a model image");

if(QString::null == fileName)

{

return;

}

//读取模板数据

QImage modelImage(fileName);

int mw = modelImage.width();

int mh = modelImage.height();

complex<double> *modelImageData = new complex<double>[mw * mh];

readImage(modelImageData, modelImage);

unsigned long t = mw * mh * 8; //根据匹配模板的大小设置一定的阀值

unsigned long m = t; //初始差异度

int ri = -1; //z左上角坐标(ri, rj)

int rj = -1;

int w = currentImage.width();

int h = currentImage.height();

if(NULL == currentImageData) //判断是否需要读取目标图像数据

{

currentImageData = new complex<double>[w * h];

readImage(currentImageData, currentImage);

}

//遍历目标图像,选取左上角坐标,考虑到模板图像的大小注意不要越界

int i, j;

for(i = 0; i < h - mh + 1; i++ )

{

for(j = 0; j < w - mw + 1; j++)

{

//下面开始对点(i, j)为左上角的mw * mh区域进行匹配

bool overFlag = false;

unsigned long k = 0; //差异值的累加和

int u, v;

for(u = 0; u < mh && !overFlag; u++)

{

for(v = 0; v < mw && !overFlag; v++)

{

k += abs(currentImageData[(i + u) * w + j + v].real()

- modelImageData[u * mw + v].real()); //计算差值并累加

if(k >= t) //判断是否大于阀值

{

overFlag = true;

}

}

}

if(k < m) //判断是否找到更加匹配的区域

{

ri = i;

rj = j;

m = k;

}

}

}

//找到匹配区域,则将目标图像匹配区域之外的点置成白色以便于观察结果

if(ri != -1)

{

for(i = 0; i < h; i++)

{

for(j = 0; j < w; j++)

{

if(i < ri || j < rj || i > ri + mh - 1 || j > rj + mw - 1)

{

currentImageData[i * w + j] = complex<double>(255, 0);

}

}

}

}

writeImage(currentImage, currentImageData);

pDispLabel->setPixmap(QPixmap(currentImage));

}

3.2. 运行截图

目标图像为图3-1,图3-2位匹配模板,匹配结果为图3-3。

clip_image027

图3-1 目标图像

clip_image028

图3-2 匹配模板

clip_image030

匹配结果


机器学习

阅读数 948

OCR字符识别

阅读数 43958

没有更多推荐了,返回首页