数字信号与图像处理实验四：图像增强技术

实验四：图像增强技术

通过本实验加深对数字图像增强操作的理解，熟悉MATLAB中的有关函数；了解直方图均衡化和卷积滤波的原理；熟悉低通和高通滤波模板的构造方法。

班级	姓名	实验名	完成时间
计算2014	lx	图像增强技术	12.8

实验要求：

1. 熟悉MATLAB的图像处理函数
2. 熟练掌握java/C#/C/C++的图像处理方法
3. 完成程序代码编写

1. 直方图均衡化

1.1 观察实验结果，写出和原图像有什么不同之处。

clear;
I = imread('Barbara.tif');
J = histeq(I);
figure; imshow(I)
figure; imshow(J)

左边是原始图像，右边是直方图均衡化后的图像，可以发现左边的图像看起来有点灰蒙蒙的（好像本来就是黑白照片），右边的图像看起来更清晰、细节更丰富，人物更加分明。

1.2 输入代码，比较I、J两幅图像的差别，写出其原因。

clear;
I=imread('Barbara.tif');
J=histeq(I);
subplot(1,2,1);imhist(I,64)
subplot(1,2,2);imhist(J,64)

左边的直方图对应了原始图像的灰度等级，右边的直方图对应了经过直方图均衡化后的图像的灰度等级。可以看出左边的原始图像灰度高频和低频信号比较少，灰度主要都集中在中间，所以图像看起来会比较平，缺少层次。

右边的图像灰度直方图范围变大，主要是利用了histeq函数对图像做了均衡化处理，图像的对比度相应也变大。

1.3 将程序中的图像换成其他图像，验证你的说法是否正确。

如图，符合上述的说法。

1.4 画出I、J的直方图，写出两者的差别。

原始图像的直方图主要集中在低频（0~20）和中高频（130~180）之间，右边的图像灰度直方图范围变大，灰度等级均匀分布在0~255之间。

由于本次实验所用图片为uint8格式，所以histeq默认使输出直方图J具有 255 个 bin 的直方图大致平坦。

2. 图像滤波

clear;
I=imread('Lena.tif');
h=ones(5,5)/25;
J=imfilter(I,h);
subplot(1,2,1); imshow(I);
subplot(1,2,2); imshow(J);

2.1 观察实验结果，写出和原图像有什么不同之处。

可以发现，右边的图像比左边的图像变得“模糊”。或者换一种说法，图像整体上变得更加平滑，边缘细节度变低。

2.2 将h分别修改下列的值，观察实验结果，写出产生此现象的原因。

h=ones(3,3)/9;

可以看出，图像同样变得“模糊”了，但是感觉没有2.1中图像模糊。

a=[1 1 1
   1 2 1
   1 1 1]
h=a/10;

同上，图像与原图几乎看不出差别。

---

要解释上面的现象，我们先来了解一下imfilter函数。

imfilter有两种用法：

• B = imfilter(A,h) 使用多维滤波器 h 对多维数组 A 进行滤波，并在 B 中返回结果。
• B = imfilter(A,h,options,...) 根据一个或多个指定的选项执行多维滤波。

我们在本题中采用的是第一种方法，即没有添加控制滤波运算的选项。可以发现，这里多次提到了**“滤波”**，下面我以2.1为例演示一下imfilter的具体流程：

首先，我们人为构造一个滤波器：

h=ones(5,5)/25;

这个滤波器长这样：

黄色的地方是这个滤波器的中心，然后我们将这个滤波器在原图像矩阵上进行滤波操作，以原图最左上角的值为例，打开matlab中的图像矩阵观察到的值如下：

对于指定偶数大小的核，核的中心为floor((size(h) + 1)/2)，即一个4*4的卷积核，其中心为[1，1]。

可以发现图像左上角的值均为162，滤波操作如下：

滤波器的中心与图像最左上角重合，然后各值与原图像上的值做积并求和，但是此时滤波器外围并没有图像填充，所以图像要向外填充，如下：

填充值默认为0，滤波器上每个值都为1/25，所以这里第一个值经过滤波后为162/25*9+0/25*16=58.32，由于图像类型为uint8，所以取整为‘58’。然后，滤波器向右移动，继续与图像上第二个值做滤波操作，如下：

用同样的方法计算得到162/25*12+0/25*13=77.76，向上取整为‘78’，然后，滤波器继续向右移动，当滤波器中心经过图片矩阵第一行所有元素后，跳转到第二行第一个元素继续上述操作。

对于第三个操作数options，官方文档给出的所有取值如下：

---

上述的是滤波器的基本操作，但是对于我们理解这道题有什么帮助呢？

首先，观察构造的5*5的滤波器，可以发现其和为1，且平均分布在滤波器中，作用就是把图像中每个像素的值转换为这个值附近5*5范围内的领域的平均值，也叫做均值滤波。所以，图像相比之前会变得更加平滑，但是细节度降低了很多。

第二个滤波器和第一个滤波器原理一致，但是范围变小了，作用就是把图像中每个像素的值转换为这个值附近3*3范围内的领域的平均值，也是均值滤波，但是由于滤波器减小了将近‘1/3’，所以图像的细节损失并不明显。

第三个滤波器和第二个滤波器基本一致，不同的地方在于中心的值变为了原来的两倍，这也就意味着图像中每个像素领域对该像素值的影响变小了，打个比方，如果图像中原来有这样一个区域：

[100,100,100,
100,200,100,
100,100,100]

其中心值在分别经过第二和第三个滤波器后的结果为100/9*8+200/9=111.1=111，100/10*8+200/5=120，可以发现后者在滤波操作后的变化更小，也就是相比之前的滤波操作保留了更多的细节。

虽然这里两者的差别并不大，但是图像上反映出来的差别却肉眼可见，可联系实验3中DCT变换的内容，图像的信息主要集中在靠近频谱中心的低频信号，所以对低频信号的细微改动都会很大程度上影响到图像。

2.3 构造一个低通滤波模板，并加以实现。

clc; clear;
img=imread('Lena.tif');
subplot(2,2,1),imshow(img); title('原图像');
% 将低频移动到图像的中心
s=fftshift(fft2(img)); 
subplot(2,2,3),imshow(log(abs(s)),[]); title('图像傅里叶变换取对数所得频谱');
[a,b] = size(img); low_filter = LowFil(a,b,s); % low_filter就是低通滤波器
subplot(2,2,4),imshow(log(abs(low_filter)),[]); title('低通滤波频谱');
new_img = uint8(real(ifft2(ifftshift(low_filter))));
subplot(2,2,2),imshow(new_img,[]); title('低通滤波后的图像');

% Parameters:
% H: 低通滤波器
% u_max: 图像第一维的长度
% v_max: 图像第二维的长度
% s: 

function [H] = LowFil(u_max, v_max, s)
% 计算图像中心点坐标以及截止频率D0
center_u = floor(u_max/2);
center_v = floor(v_max/2);
D0 = min(center_u,center_v)/12;
D0 = D0^2;

% 初始化
H=zeros(u_max, v_max);

for u=1:u_max
    for v=1:v_max
        % (u, v)和中心点间的距离
        dist = (u - center_u)^2 + (v - center_v)^2;
        % 判断是否要截止
        if(dist<=D0)
            h = 1;
        else
            h = 0;
        end
        H(u, v) = h * s(u, v);
    end
end

使用函数fft2将空间图像进行二维傅里叶变换后，转化到频域上，我们可以得到这个图像每个像素的相位和幅度值，对于相位而言，代表着图像中的位置形状信息，一般情况下我们不对其做处理，否则可能就无法恢复出原始图像的形状。
而幅度值则主要代表着能量的大小，也就是每一个频率上的能量大小。由于本题中我们要构造低通滤波器，即去除图像的高频部分，公式如下：
$$H(u,v)=\{\begin{array}{l}1D(u,v)\le D0,\\ 0D(u,v)>D0\end{array}$$

D0为截止频率，H(u,v)即为滤波器

$$G(u,v)=H(u,v)F(u,v)$$

F(u,v)为经过傅里叶变换后的频域图，G(u,v)为低通滤波频谱图。

需要注意的是，在进行傅里叶变换（fft2）后，低频部分集中分布在四个角落，通过移位（fftshift）之后，将其集中到图像的中心，这只是为了便于处理。当我们对幅值进行处理后，再进行反移位（ifftshift）恢复到原始的分布情况，

Y = fftshift(X) 通过将零频分量移动到数组中心，重新排列傅里叶变换 X。

如果 X 是矩阵，则 fftshift 会将 X 的第一象限与第三象限交换，将第二象限与第四象限交换。

之后结合相位计算出每个像素对应的值。对其进行傅里叶反变换（ifft）后，就可以恢复出时域图像，此时的图像就是经过低通滤波器处理后的图像了，如下：

需要注意，在经过傅里叶变换后，频域部分为复数，即使在傅里叶逆变换后，也没有变为实数，只是将复数部分变为0，所以最好要加一个real()函数，将频域图转为实数，不转的话会有如下警告。

2.4 构造2—3个高通滤波模板，并加以实现，和低通滤波的效果加以对比，写出高通滤波和低通滤波的差别，并写出低通滤波模板、高通滤波模板构造的规律。

首先是最简单的高通滤波，直接在上一道题的基础上改动一下：

for u=1:u_max
    for v=1:v_max
        % (u, v)和中心点间的距离
        dist = (u - center_u)^2 + (v - center_v)^2;
        % 判断是否要截止
        if(dist<=D0)
            h = 0;
        else
            h = 1;
        end
        H(u, v) = h * s(u, v);
    end
end

只要将原来的截止频率的判断改变一下，就成为了一个理想高通滤波器：
$$H(u,v)=\{\begin{array}{l}0D(u,v)\le D0,\\ 1D(u,v)>D0\end{array}$$

其余部分不用改变，生成的图片如下：

接着实现比较复杂的高通滤波器，这里以高斯高通滤波器为例，公式如下：
$$H(u,v)=1-e^{-D^2(u,v)/2D_0^2}$$

for u=1:u_max
    for v=1:v_max
        d = ((u-center_u)^2+(v-center_v)^2);
        H(u,v) = 1 - exp(-(d)/(2*(D0^2)));
    end
end
H = H .* s;

其余部分基本没有变动。

总的来说，模板构造显然是有规律的，主要步骤如下：

1. 将空间图像进行傅里叶变换转换为频域图，然后移位处理，将低频部分集中到图像中心。
2. 对频域图做滤波处理（自己编写滤波器）
3. 将滤波后的频域图进行反移位，傅里叶逆变换，恢复成原始图像。

3. 编程实现直方图均衡化函数（histeq）和卷积滤波函数（imfilter）

直方图均衡化：

#include<opencv2\opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;

Mat hist(Mat &image) {
	double gray[256] = { 0 };
	Mat temp = Mat(image.rows, image.cols, CV_8UC1);
	for (int i = 0; i < image.rows; i++) {
		for (int j = 0; j < image.cols; j++) {
			int k = image.at<uchar>(i, j);
			gray[k]++;
		}
	}
	for (int x = 0; x < 256; x++) {
		gray[x] = (double)gray[x] / (image.rows*image.cols);
	}

	for (int m = 0; m < 256; m++) {
		gray[m + 1] = gray[m + 1] + gray[m];
	}
	for (int a = 0; a < image.rows; a++) {
		for (int b = 0; b < image.cols; b++) {
			int k = image.at<uchar>(a, b);
			temp.at<uchar>(a, b) = (uchar)(255 * gray[k] + 0.5);
		}
	}
	imshow("hist", temp);
	return temp;
}
void main() {
	Mat image1 = imread("H:/大三上/数字信号/实验4/input/Cameraman.tif", 0);
	imshow("Origin",image1);
	hist(image1);
	waitKey(0);
	getchar();
}

效果如下：

卷积滤波函数：

#include <iostream>
#include <math.h>
#include<opencv2\opencv.hpp>
using namespace std;
using namespace cv;

void Sharpen(const Mat& myImage, Mat& Result)
{
	CV_Assert(myImage.depth() == CV_8U);        //判断函数CV_Assert
	const int nChannels = myImage.channels();

	for (int j = 1; j < myImage.rows - 1; ++j)
	{
		const uchar* precious = myImage.ptr<uchar>(j - 1);      //当前像素上一行指针
		const uchar* current = myImage.ptr<uchar>(j);           //当前像素行指针
		const uchar* next = myImage.ptr<uchar>(j + 1);          //当前像素下一行指针

		uchar* output = Result.ptr<uchar>(j);

		//利用公式和该中心点邻域的8个像素对当前像素值进行处理
		for (int i = nChannels; i < nChannels * (myImage.cols - 1); ++i)
		{
			// 3*3的均值滤波器
			*output++ = saturate_cast<uchar>
				(floor((current[i] + current[i - nChannels] + current[i + nChannels]
				+ precious[i] + precious[i - nChannels] +  precious[i + nChannels]
				+ next[i] + next[i - nChannels] + next[i + nChannels])/9.0));
		}
	}
	Result.row(0).setTo(Scalar(0));                 //设置第一行所有元素值为0
	Result.row(Result.rows - 1).setTo(Scalar(0));     //设置最后一行所有元素值为0
	Result.col(0).setTo(Scalar(0));                 //设置第一列所有元素值为0
	Result.col(Result.cols - 1).setTo(Scalar(0));     //设置最后一列所有元素值为0
}

void main()
{
	Mat srcImg = imread("H:/大三上/数字信号/实验4/input/Cameraman.tif", 0);
	imshow("Src", srcImg);
	Mat dstImg;
	dstImg.create(srcImg.size(), srcImg.type());
	Sharpen(srcImg, dstImg);
	imshow("Dst",dstImg);
	waitKey(0);
}

上述代码使用了3*3的均值滤波器对图像进行滤波处理，算法和2.3题讲解的一致，效果如下图：

可以看出右边的图像相比左边的图像变得更加模糊。

总结

本次实验对于我理解图像增强起到了很大的帮助，尤其是图像滤波哪个部分，在查阅了很多的文献后，发现滤波器的构造方法有规律可循的时候，内心还是比较激动的。当然，实验中也出现了很多的问题，在实现最后一问卷积滤波函数的时候，其实还是有不足的地方，没办法自定义滤波器，其实可以再花一些时间写个自定义滤波器的函数，但是由于期末时间比较紧张，留给我自由发挥的时间比较有限，所以便不了了之了。

总的来说，本次实验很好地激发了我对图像增强的兴趣，希望我以后有时间可以静下心来做一些这方面的研究。