数字图像处理的基础

将彩色图转换为灰度图

代码:

1
2
3
4
RGB = imread('D:\reports\test.jpg')
subplot(121),imshow(RGB)
gray = rgb2gray(RGB)
subplot(122),imshow(gray)

转换的结果是:

image-20220731101025493

将灰度图转为为索引图

代码:

1
2
3
4
5
6
RGB = imread('D:\reports\test.jpg')
subplot(221),imshow(RGB)
gray = rgb2gray(RGB)
subplot(222),imshow(gray)
[X,map] = gray2ind(gray,8)
subplot(223),imshow(X,map)

结果:

image-20220731101632758

结论:变为索引值之后图像变得非常的不自然了,图中的显色有分块的现象,特别是在水中的那部分,这是为什么呢?将函数中索引值8增大,这个值越大发现图像之中的颜色表示就更加的自然(最后那幅图的索引值是64的)

索引值:图像颜色映射的个数,也表示图像的灰度级数

image-20220731102131472

将索引图转化为真彩色的图

代码:

1
2
3
4
5
6
7
8
RGB = imread('D:\reports\test.jpg')
subplot(221),imshow(RGB)
gray = rgb2gray(RGB)
subplot(222),imshow(gray)
[X,map] = gray2ind(gray,64)
subplot(223),imshow(X,map),
color = ind2rgb(X,map)
subplot(224),imshow(color)

这样子得到的rgb和原来的二值图是一样的,也不知道为什么,并且按照网上的方式读取tif文件,一直不能够读取成功,使用bmp格式的图像得到的结果仍然是一样的,所以二值图转换的原理是什么需要考究一下

将灰度图像转化为二值化图

代码:

1
2
3
4
5
6
7
8
color = imread('D:\reports\test.jpg')
subplot(221),imshow(color)
gray = rgb2gray(color)
subplot(222),imshow(gray)
BW1 = im2bw(gray,0.3)
subplot(223),imshow(BW1)
BW2 = im2bw(gray,0.7)
subplot(224),imshow(BW2)

运行的结果:

image-20220731105415590

思考:灰度转二值图的函数是im2bw,后面的阈值的意思是,判断哪些该转换成1,哪些该转换成0,阈值设置得越小,得到的二值图像的白色面积越多。也就是进行阈值处理。

图像的几何变换

图像的平移

原理说明:

image-20220731145937610

代码:

mymove.m(脚本文件,在这里面定义了一个用于移动的函数,在同一文件夹目录下可以直接调用该函数)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function OUT = mymove(Image,move_X,move_Y)
J = double(Image); % 将图像中的数据转换成双精度的形式
HW=size(J);% 获取原图像的大小
OUT = zeros(HW); % 新建新图像的矩阵
OUT(1:HW(1),1:HW(2)) = inf; % inf在matlab中表现的是无限大的形式,即这里将图像的内容初始化为全白的背景
if((move_X>=0)&&(move_Y>=0))
OUT(move_X+1:HW(1),move_Y+1:HW(2)) = J(1:HW(1)-move_X,1:HW(2)-move_Y);
elseif((move_X>0)&&(move_Y<0))
OUT(move_X+1:HW(1),1:HW(2)+move_Y) = J(1:HW(1)-move_X,1-move_Y:HW(2));
elseif((move_X<0)&&(move_Y>0))
OUT(1:HW(1)+move_X,1+move_Y:HW(2)) = J(1-move_X:HW(1),1:HW(2)-move_Y);
elseif((move_X<0)&&(move_Y<0))
OUT(1:HW(1)+move_X,1:HW(2)+move_Y) = J(1-move_X:HW(1),1-move_Y:HW(2));
end % 判断性语句结束的标志
OUT = uint8(OUT);
end

主程序代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = mymove(image,60,100);
subplot(232),imshow(image_1);
image_2 = mymove(image,-60,100);
subplot(234),imshow(image_2);
image_3 = mymove(image,60,-100);
subplot(235),imshow(image_3);
image_4 = mymove(image,-60,-100);
subplot(236),imshow(image_4);

代码运行的结果:

image-20220731145831362

图像的缩放

原理:

使用的双线性插值算法,遍历变换之后的每个像素点,并且定位和该像素点相关的四个像素点,利用这四个像素点的像素根据公式推导出该像素点的像素。计算的公式如下图所示

image-20220731162031267

代码:

脚本函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function OUT = changeSize(Image,zoom_X,zoom_Y)
J = double(Image); % 将图像中的数据转换成双精度的形式
HW = size(J);% 获取原图像的大小 返回的第一个参数是行数 第二个参数是列数
OUT = zeros(HW); % 新建新图像的矩阵
OUT(1:HW(1),1:HW(2)) = inf; % inf在matlab中表现的是无限大的形式,即这里将图像的内容初始化为全白的背景
changeHW = [floor(HW(1)*zoom_Y),floor(HW(2)*zoom_X)];%放大缩小后的图像的大小
%便利放大缩小之后的图像的每一个像素点
for i = 1:changeHW(1)
for j = 1:changeHW(2)
% 变化之后的图像中的像素点定位到原图像中的对应的位置
X = i/zoom_Y;
Y = j/zoom_X;
% 获得数值的小数部分
u = X-floor(X);
v = Y-floor(Y);
if X< 1
X =1;
end

if Y< 1
Y = 1;
end
OUT(i,j) = J(floor(X),floor(Y))*u*v + J(ceil(X),floor(Y))*u*(1-v) + J(floor(X),ceil(Y))*(1-u)*v + J(ceil(X),ceil(Y))*(1-u)*(1-v);
end
end
OUT = uint8(OUT(1:HW(1),1:HW(2)));
end

主函数

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = changeSize(image,0.5,0.5);
subplot(232),imshow(image_1),axis on;
image_2 = changeSize(image,1.5,1.5);
subplot(233),imshow(image_2),axis on;
image_3 = changeSize(image,0.5,0.8);
subplot(234),imshow(image_3),axis on;
image_4 = changeSize(image,1.0,1.5);
subplot(235),imshow(image_4),axis on;

代码运行的结果:

image-20220731161451618

图像旋转

新图像的高和宽的公式:

**nHW(1)=ceil(HW(1)*abs(cos(abs(alpha)))+HW(2)*abs(sin(abs(alpha)))); nHW(2)=ceil(HW(1)abs(sin(abs(alpha)))+HW(2)abs(cos(abs(alpha))));

换算到逆时针旋转之后的坐标公式:
x1=x0cos(β)-y0sin(β)
y1=x0sin(β)+y0cos(β)

图像:

image-20220731163223386

则由新图像中的坐标换算到原图像中的坐标公式:
x=icos(β)+jsin(β)
y=-isin(β)+jcos(β)

变换之后的图像中的每个坐标不只是有角度上的转变,还有一定的平移,因为要确保图像在中央

旋转函数代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
function OUT = myrotate(image,alpha) % alpha是旋转的角度
J = double(image);% 将图像中的值转换成双精度形式
HW = size(J);% 获得图像的大小 size
alpha = mod(alpha,360);% 将角度转变到0~360
alpha = (alpha*pi)/180;% 将角度转变成弧度制
nHW(1) = ceil(HW(1)*abs(cos(abs(alpha)))+HW(2)*abs(sin(abs(alpha))));% 旋转之后的图像的大小(长和宽都要重新定义)
nHW(2) = ceil(HW(1)*abs(sin(abs(alpha)))+HW(2)*abs(cos(abs(alpha))));
OUT = zeros(nHW);% 利用新的长和宽定义图像的矩阵
OUT(1:nHW(1),1:nHW(2)) = inf;% 对这个矩阵进行初始化(inf 初始化成白色的)
u0 = HW(2)*abs(sin(alpha));u2 = HW(1)*abs(cos(alpha));% 竖直方向的平移量
u1 = HW(1)*abs(sin(alpha));u3 = HW(2)*abs(cos(alpha)); % 水平方向的平移量
T = [cos(alpha),sin(alpha);-sin(alpha),cos(alpha)];% 对坐标进行换算的矩阵
% 遍历这个新的图像
for i = 1:nHW(1)
for j = 1:nHW(2)
if(alpha>=0 && alpha<pi/2)
XY = T * ([i;j]-[u0;0]);% 此时整个图像需要向上移动
elseif(alpha>=pi/2 && alpha<pi)
XY = T *([i;j]-[u0+u2;u3]);
elseif(alpha>=pi && alpha<3*pi/2)
XY = T *([i;j]-[u2;u1+u3]);
elseif(alpha>=3*pi/2 && alpha<2*pi)
XY = T * ([i;j]-[0;u1]);
end
X = XY(1);
Y = XY(2);
if X>=1 && X<=HW(1) && Y>=1 && Y<=HW(2)
u = X - floor(X);
v = Y - floor(Y);
OUT(i,j) = J(floor(X),floor(Y))*u*v + J(ceil(X),floor(Y))*u*(1-v) + J(floor(X),ceil(Y))*(1-u)*v + J(ceil(X),ceil(Y))*(1-u)*(1-v);
end
end
end
OUT = uint8(OUT(1:nHW(1),1:nHW(2)));
% 根据旋转的角度平移坐标值 并且 乘以换算的矩阵 得到在原图像中的坐标点
% 得到小数的部分
% 使用双线插入法获得对应的像素值
% 将OUT的值转变成uint8的形式
end

主函数代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image1 = myrotate(image,30);
subplot(232),imshow(image1),axis on;
image1 = myrotate(image,120);
subplot(233),imshow(image1),axis on;
image1 = myrotate(image,210);
subplot(234),imshow(image1),axis on;
image1 = myrotate(image,310);
subplot(235),imshow(image1),axis on;

运行的效果图:

image-20220731175016882

图像镜像

水平镜像的公式:(水平镜像操作是以原图像的垂直中轴线为中心,将图像分为左右两部分进行对称变换)
x1=x0
y1=width+1-y0

垂直镜像的公式:(垂直镜像操作是以原图像的水平中轴线为中心,将图像分为上下两部分进行对称变换)

x1=height+1-x0
y1=y0

width =》 第二个得到的数据

height =》第一个得到的数据

镜像函数的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function OUT = mymirror(Image,choice)
J= double(Image);% 将图像的数据进行双精度
HW = size(J);% 获得图像的大小 size
OUT = zeros(HW(1),HW(2));% 利用图像的大小创建对应的数组
OUT(1:HW(1),1:HW(2)) = inf;% 初始化矩阵 inf
% 遍历矩阵
for i = 1:HW(1)
for j = 1:HW(2)
if strcmp(choice,'level') % 水平镜像
OUT(i,j) = J(i,HW(2)+1-j);
elseif strcmp(choice,'vertical') % 垂直镜像
OUT(i,j) = J(HW(1)+1-i,j);
end
end
end
% 针对不同的选择进行设置
% 对相应位置图像的像素进行设置
OUT = uint8(OUT);
end

主函数的代码:

1
2
3
4
5
6
7
8
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = mymirror(image,'level')
subplot(232),imshow(image_1),axis on,title('水平镜像')
image_2 = mymirror(image,'vertical')
subplot(233),imshow(image_2),axis on,title('垂直镜像')

效果图:

image-20220731235124610

图像的转置

转置就是像素点的x坐标和y坐标进行交换,注意这样使得它的长和宽的大小也进行了交换

转置函数代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
function OUT = myTranspose(Image)
J= double(Image);% 将图像的数据进行双精度
HW = size(J);% 获得图像的大小 size
OUT = zeros(HW(2),HW(1));% 因为是转置所以它的行宽的大小进行了交换
OUT(1:HW(2),1:HW(1)) = inf;% 初始化矩阵 inf
% 遍历矩阵
for i = 1:HW(2)
for j = 1:HW(1)
OUT(i,j) = J(j,i); % 转置像素赋值的公式
end
end
OUT = uint8(OUT);
end

主函数代码:

1
2
3
4
5
6
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = myTranspose(image);
subplot(232),imshow(image_1),axis on,title('转置的图像')

效果图:
image-20220731235832423

图像的剪切

剪切函数代码:

1
2
3
4
5
6
function OUT = myCut(Image,strat_X,strat_Y)
J= double(Image);% 将图像的数据进行双精度
HW = size(J);
OUT = J(strat_X:HW(1),strat_Y:HW(2));% 因为是转置所以它的行宽的大小进行了交换
OUT = uint8(OUT);
end

主函数代码:

1
2
3
4
5
6
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = myCut(image,50,70);
subplot(232),imshow(image_1),axis on,title('剪切的图像')

效果图

image-20220801000929925

灰度变换

函数映射

使用imadjust这个灰度值映射函数

语法:g=imadjust(f,[low-in,hign-in],[low-out,high-out],gamma);

将输入图像指定灰度范围映射到想要得到的输出图像的灰度范围

gamma为映射函数曲线的特征参数,gamma小于1,则映射被加权至高的输出值,大于1,则被加权至较低的输出值,等于1为线性映射。gamma默认为1。

运行代码

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
image_1 = imadjust(image,[0,1],[1,0]);
subplot(232),imshow(image_1)
image_2 = imadjust(image,[0.3,0.8],[0,1]);
subplot(233),imshow(image_2)
image_3 = imadjust(image,[0.3,0.8],[0,1],0.5);
subplot(234),imshow(image_3)
image_4 = imadjust(image,[0.3,0.8],[0,1],1.5);
subplot(235),imshow(image_4)

效果图:

image-20220801002104962

分析:第二张图是将[0,1]的像素值投影到[1,0]的像素值区域,这样得到的图片的效果就是颜色取反,越黑的像素点转变之后就越白。第三张图使得图像黑色区域更黑,白色区域更白,容易使得图像整体的感觉变得更暗,所以有第四张图,将grammer设置为0.5(数值小于1了),能够提高图像的亮度,第五张图,将grammer设置为1.5,降低了图像的亮度。

直方图处理

绘制图像直方图

代码:

1
2
3
4
5
6
7
8
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
subplot(232),imshow(image)
imhist(image,16)
subplot(233),imshow(image)
plot(imhist(image))

代码的运行结果:

image-20220801004113238

均衡化调用的函数hiseq()

运行的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
subplot(231),imshow(image)
subplot(232),imshow(image)
imhist(image,16)
subplot(233),imshow(image)
plot(imhist(image))
g = histeq(image,256)
subplot(234),imshow(g)
subplot(235),imshow(g)
imhist(g,16)
subplot(236),imshow(g)
plot(imhist(g))

效果图

image-20220801004646948

分析:通过图我们可以看到均衡化之后的直方图中每部分像素值(每个分类部分)中的点的数量近似相同,图像显示也有了明显的差异

使用直方图图匹配

参考的文章:https://www.programminghunter.com/article/37032194016/

自己写的代码:(代码里面找最近的点的那个地方,使用了index()函数自己没有看懂)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image)
% 双峰高斯函数
r = 127;
x = -r:r+1;
sigma = 20;
y1=exp(-((x-80).^2)/(2*sigma^2));
y2=exp(-((x+80).^2)/(2*sigma^2));
y=y1+y2;
y = y /sum(y);% 归一化 概率
subplot(232),imshow(y)
plot(y)
% 该函数的累计直方图
G = [];
for i=1:256
G = [G sum(y(1:i))];
end

subplot(233),imshow(image)
imhist(image)
% 获得原图像的累计直方图
hist = imhist(image);% 获得原图像的直方图
imageP = hist/(HW(1)*HW(2));% 每一中像素的概率
% 对概率再累计求和
s=[];
for i=1:256
s=[s sum(imageP(1:i))];
end

% 找到两个概率累计图中最近的点
for i = 1:256
tem{i} = G-s(i);
tem{i} = abs(tem{i}); % 这里看的是距离,所以取绝对值进行比较就可以了
[a index(i)] = min(tem{i}); % 获得距离最小的点
end

newimage=zeros(HW(1),HW(2));
for i = 1:HW(1)
for j = 1:HW(2)
newimage(i,j) = index(image(i,j)+1)-1;
end
end

newimage = uint8(newimage);
subplot(234),imshow(newimage)
subplot(234),imshow(newimage)
imhist(newimage)

运行的结果是:
image-20220801104042622

分析:看这个曲线就能知道匹配了上面的图像了,这里面要使用累积概率进行计算才可以

图像的平滑处理

采用均值和高斯滤波器进行平滑处理

用来生成滤波器的函数

h = fspecial(type)
h = fspecial(type,para)

type参数来指定滤波器的种类,para来对具体的滤波器种类添加额外的参数信息

高斯滤波器

H = fspecial(‘gaussian’,hsize,sigma) :type = ‘gaussian’时就是高斯滤波器了,size指定滤波器的大小,默认值是3×3,sigma指定滤波器的标准差,默认值是0.5(决定了高斯模糊核的模糊程度),使用sigma的值很大的时候,然后运用此模糊核对图像处理,会使图像更加模糊。

例: h1 = fspecial(‘gaussian’,5,1)

均值滤波器

H =fspecial(‘average’,hsize) :hsize指定滤波器的尺寸

例:h1 = fspecial(‘average’,5) %5x5的核

圆形均值滤波器

H = fspecialL(‘disk’,RADIUS),radius代表区域半径。

例:h1 = fspecial(‘disk’,3)

拉普拉斯滤波器

H = fspecial(’laplacian’,ALPHA):alpha用于控制算子形状,取值范围为【0,1】,默认值为0.2,尺寸一定是3×3

例:h1 =fspecial(‘laplacian’)

拉普拉斯高斯滤波器

H = fspecial(‘log’,hsize,sigma) :hsize表示模板尺寸,默认值为3x3,sigma为滤波器的标准差

例:h1 =fspecial(‘log’,3,0.2)

prewitt

H = fspecial(‘prewitt’) :用于边缘增强,这个函数是没有参数的,默认的大小就是3x3,

sobel

H = fspecial(‘sobel’):这个函数是无参数的,用于边缘的提取。

例:h1 = fspecial(‘sobel’)

高斯和均值滤波器的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
average3image = filter2(fspecial('average',3),image)/255;% 采用均值滤波器进行平滑处理 3x3
subplot(232),imshow(average3image),title('3x3的均值滤波器')
average5image = filter2(fspecial('average',5),image)/255;% 5x5
subplot(233),imshow(average5image),title('5x5的均值滤波器')
average7image = filter2(fspecial('average',7),image)/255;% 7X7
subplot(234),imshow(average7image),title('7x7的均值滤波器')
gaussianImage1 = filter2(fspecial('gaussian'),image)/255;
subplot(235),imshow(gaussianImage1),title('0.5的高斯滤波器')
gaussianImage2 = filter2(fspecial('gaussian',5,10),image)/255;
subplot(236),imshow(gaussianImage2),title('10的高斯滤波器');

效果图:
image-20220801110943637

分析:可以看到越大的模板,处理之后的图像就会越模糊,sigma越大就会使得页面更加的模糊

边缘增强的两个滤波器prewitt和sobel:

代码

1
2
3
4
5
6
7
8
9
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
prewittImage = filter2(fspecial('prewitt'),image)/255
subplot(232),imshow(prewittImage),title('prewitt边缘增强')
sobelImage = filter2(fspecial('sobel'),image)/255
subplot(233),imshow(sobelImage),title('sobel边缘增强')

锐化处理的公式(把边缘部分进行突出)

  • 原图*2-平滑图像
  • 原图+边缘处理图像

锐化处理代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
average5image = filter2(fspecial('average',5),image)/255;% 5x5
subplot(234),imshow(average5image),title('5x5的均值滤波器')
prewittImage = filter2(fspecial('prewitt'),image)/255;
subplot(232),imshow(prewittImage),title('prewitt边缘增强')
sobelImage = filter2(fspecial('sobel'),image)/255;
subplot(233),imshow(sobelImage),title('sobel边缘增强')
subImage1 = image*2-uint8(average5image);
subplot(235),imshow(subImage1),title('原图*2-平滑图像锐化')
subImage2 = image+uint8(prewittImage);
subplot(236),imshow(subImage2),title('原图+边缘处理图像锐化')

效果图:
image-20220801112456585

分析:因为我用的是模糊程度偏大的图像,所以使得使用第一种方法锐化之后的结果,图像的整体的亮度会偏大,第二种方式锐化之后的,仔细看能够感觉细节上面增强了,但是直观上来看,我认为是不明显的。

边缘检测

边缘的两个性质:

  • 强度的一阶导数的模大于某个阈值的位置(在边缘处像素几句变化)
  • 强度的二阶导数有过零点的位置

使用edge()函数检测 edge()函数的各种用法

sobel和canny边缘检测的对比

代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
BW1 = edge(image,'sobel');
subplot(232),imshow(BW1),title('sobel边缘检测')
BW2 = edge(image,'canny');
subplot(233),imshow(BW2),title('canny边缘检测')
figure
imshowpair(BW1,BW2,'montage')
title('sobel边缘检测 canny边缘检测 ')

效果图:

image-20220801114928541

分析:canny边缘的提取能力高于其它的

边缘检测的方法:
robert:边缘定位精度较高,对于陡峭边缘且噪声低的图像效果较好,但没有进行平滑处理,没有抑制噪声的能力

sobel prewitt:进行了平滑处理,对噪声具有一定抑制能力,但容易出现多像素宽度

laplacian:对噪声较为敏感,使噪声能力成分得到加强,容易丢失部分边缘方向信息,造成一些不连续的检测边缘,同时抗噪声能力较差。

log:抗噪声能力较强,但会造成一些尖锐的边缘无法检测到。

canny:最优化思想的边缘检测算子,同时采用高斯函数对图像进行平滑处理,但会造成将高频边缘平滑掉,造成边缘丢失,采用双阈值算法检测和连接边缘。

上面几种边缘检测方法的对比:

image-20220801120104487

分析:对比上面的几个图,我们能够看到canny的边缘检测的效果是最好的。

图像的分割

这里采用的方法都是阈值分割,上面用于边缘检测也是可以用于图像分割的

全局分割

代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
subplot(232),imshow(image),title('原图直方图')
imhist(image,200)
% 采用全局阈值分割
grobeCutImage1 = image>100;
subplot(234),imshow(grobeCutImage1),title('分割图像一')
grobeCutImage2 = image>180;
subplot(235),imshow(grobeCutImage2),title('分割图像二')

效果图:
image-20220801151058135

Otsu阈值分割

代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
subplot(232),imshow(image),title('原图直方图')
imhist(image,200)
% 采用Otsu阈值分割
doubleImage = im2double(image);
T = graythresh(doubleImage) ;
J = im2bw(doubleImage,T);
subplot(234),imshow(J),title('Otsu阈值处理')

效果图:

image-20220801151842037

迭代式阈值分割法

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
subplot(232),imshow(image),title('原图直方图')
imhist(image,200)
% 采用迭代式阈值分割法
doubleImage = im2double(image);
thred = 0.01 ;% 精度
T1 = (min(doubleImage)+max(doubleImage));% 得到第一个分界点 最大值和最小值取二分之一的位置(T1)
sumRight = find(doubleImage<T1);% 找出小于这个分界点的值
sumLeft = find(doubleImage>=T1);% 找出大于这个分界点的值
T2 = (mean(doubleImage(sumRight))+mean(doubleImage(sumLeft)))/2;% 分别求这两边数字的平均值
% 对这两边的平均值求和再*1/2
while abs(T1-T2)<=thred
T1 = T2;
sumRight = find(doubleImage<T1);
sumLeft = find(doubleImage>=T1);
T2 = (mean(doubleImage(sumRight))+mean(doubleImage(sumLeft)))/2;
end % whiel循环,T1和T2之间的差在0.01之上时不断进行迭代
% T1 =T2 再继续用相同的方法求得T2(这样保证左右两边的像素总值近似相等时结束)
% 使用im2bw()进行阈值的划分
result = im2bw(doubleImage,T2);
subplot(232),imshow(result),title('迭代式阈值分割')

效果图:

image-20220801154509823

分析:这里的这个阈值点的确定的要求,要使得这个阈值以下的所有点像素相加的结果等于这个阈值以上的所有点像素相加的结果,但我觉得分割效果并不是很优秀

区域分割的分水岭分割

代码:

1
2
3
4
5
6
7
8
9
10
close all;
RGB = imread('D:\reports\test.jpg');
image = rgb2gray(RGB);
HW = size(image);
subplot(231),imshow(image),title('原图')
subplot(232),imshow(image),title('原图直方图')
imhist(image,200)
% 采用分水岭分割
newimage = watershed(image,1); % 这个地方能输入的参数值是 1, 4, 6, 8, 18, or 26 默认是8 指的是连通区的个数
subplot(234),imshow(newimage),title('分水岭分割')

效果图:

image-20220801155205365

感觉不到任何的效果,不知道这种分割方法的效果和意义所在

图像的形态学处理

图像腐蚀

运行的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
clear all;close all;
bw = zeros(9,9);
bw(2:6,3:7) = 1;% 给2到6排和3到7列的数据赋值为1
se = strel('square',5);% 创立5*5的结构体
se2 = strel('square',1);% 创立1*1的结构体
se3 = strel('square',3);% 创立3*3的结构体
bw2 = imerode(bw,se);% 用se模板对bw进行腐蚀的操作
bw3 = imerode(bw,se2);% 用se模板对bw进行腐蚀的操作
bw4 = imerode(bw,se3);% 用se模板对bw进行腐蚀的操作
subplot(221); imshow(bw);title('输入的图像')
subplot(222); imshow(bw3);title('1*1 腐蚀后的图像')
subplot(223); imshow(bw4);title('3*3 腐蚀后的图像')
subplot(224); imshow(bw2);title('5*5 腐蚀后的图像')

效果图:

image-20220804110512155

3*3的腐蚀模板:

image-20220804110649716

原图像的模板:

image-20220804110723961

3*3腐蚀之后的结果:

image-20220804110756214

图像膨胀

运行的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
clear all;close all;
bw = zeros(9,9);
bw(4:5,4:6) = 1;% 给4到5排和4到6列的数据赋值为1
se = strel('square',5);% 创立5*5的结构体
se2 = strel('square',1);% 创立1*1的结构体
se3 = strel('square',3);% 创立3*3的结构体
bw2 = imdilate(bw,se);% 用se模板对bw进行膨胀的操作
bw3 = imdilate(bw,se2);% 用se模板对bw进行膨胀的操作
bw4 = imdilate(bw,se3);% 用se模板对bw进行膨胀的操作
subplot(221); imshow(bw);title('输入的图像')
subplot(222); imshow(bw3);title('1*1 膨胀后的图像')
subplot(223); imshow(bw4);title('3*3 膨胀后的图像')
subplot(224); imshow(bw2);title('5*5 膨胀后的图像')

图像运行结果的截图:

image-20220804111645824

线性膨胀图像

代码:

1
2
3
4
5
6
7
8
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('line',11,90);% 制定一个线性的结构,长度是11,角度是90度
bw1 = imdilate(image,se);
subplot(232);imshow(bw1);title('垂直线性膨胀图像')

效果图:
image-20220804112706637

用于膨胀的模板

image-20220804112819360

球形膨胀

代码:

1
2
3
4
5
6
7
8
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('ball',5,5);% 制定一个球形膨胀的结构
bw1 = imdilate(image,se);
subplot(232);imshow(bw1);title('球形膨胀图像')

效果图:

image-20220804113013278

球形模板图:

image-20220804113314104

开运算

腐蚀运算,然后再做膨胀运算

运算代码:

1
2
3
4
5
6
7
8
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('disk',3);% 制定一个棋盘膨胀的结构
bw1 = imopen(image,se); % 进行开运算
subplot(232);imshow(bw1);title('开运算图像')

效果图:

image-20220804113806790

分析:

我们能够看到的是右边的图相比于左边的图已经不连贯而变得模糊了,是因为开运算的第一次腐蚀,就会将很多细节的地方去除掉。

闭运算

先做膨胀运算,然后再做腐蚀运算

代码:

1
2
3
4
5
6
7
8
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('disk',3);% 制定一个棋盘膨胀的结构
bw1 = imclose(image,se); % 进行闭运算
subplot(232);imshow(bw1);title('闭运算图像')

效果图:

image-20220804114257142

形态学的高帽滤波

高帽滤波定义

image-20220804114443334

A输入的图像,B为采用的结构元素,即从图像中减去形态学开操作后的图像,通过高帽滤波可以增强图像的对比度

代码:

1
2
3
4
5
6
7
8
9
10
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('disk',3);% 制定一个棋盘膨胀的结构
bw1 = imtophat(image,se); % 高帽滤波
adjustBw1 = imadjust(bw1);% 进行图像增强的操作
subplot(232);imshow(bw1);title('高帽滤波后的图像')
subplot(233);imshow(adjustBw1);title('图像增强后图像')

效果图:

image-20220804115244654

分析:因为高帽滤波的减法运算,所以图像的背景色就是0了

低帽滤波

定义:

image-20220804115641941

A输入的图像,B为采用的结构元素,即从图像中减去形态学闭操作后的图像。通过低帽滤波可以获取图像的边缘

代码:

1
2
3
4
5
6
7
8
9
10
11
12
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
HW = size(image);
subplot(231);imshow(image);title('原图像')
se = strel('disk',3);% 制定一个棋盘膨胀的结构
bw1 = imtophat(image,se); % 高帽滤波
bw2 = imbothat(image,se); % 低帽滤波
dealBw2 = imsubtract(imadd(image,bw1),bw2);% 进行图像处理的操作
subplot(232);imshow(bw1);title('高帽滤波后的图像')
subplot(233);imshow(bw2);title('低帽滤波后的图像')
subplot(234);imshow(dealBw2);title('处理后的图像')

效果图:

image-20220804120236113

imsubtract(imadd(image,bw1),bw2):其中bw1进行的是高帽滤波,bw2进行的是低帽滤波,这样处理的结果使得图像的对比度增强

填充操作

对灰度图像进行填充

代码:

1
2
3
4
5
6
7
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
% bw = im2bw(image);
fillBw = imfill(image,'holes');
subplot(121),imshow(image),title('二值化图像')
subplot(122),imshow(fillBw),title('填充操作之后的图像')

效果图:

image-20220804142511067

分析:右边的图相比于左边哦我们能看到a和e这种中间有镂空的空洞部分都被进行了填充

灰度图像中设定阈值的局部极大值

[输入像素值-周围值(10)]>X(设定阈值)时,这个像素点就会被特别标记

通过imhmax()对极大值进行抑制

代码:

1
2
3
4
5
6
7
8
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
R = imextendedmax(image,100); % 带阈值的局部区域最大
R2 = imextendedmax(image,150);
subplot(131),imshow(image),title('原图像')
subplot(132),imshow(R),title('设定阈值(100)');
subplot(133),imshow(R2),title('设定阈值(150)');

效果图:

image-20220804143530355

灰度图像边缘测定

使用膨胀腐蚀

代码:

1
2
3
4
5
6
7
8
9
10
11
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
se = strel('disk',3);% 创模板
erodeImage = imerode(image,se);% 腐蚀操作
dilateImage = imdilate(image,se);% 膨胀操作
dealImage1 = dilateImage - erodeImage;
dealImage2 = image - erodeImage;
subplot(221),imshow(image),title('原图像')
subplot(222),imshow(dealImage1),title('膨胀-腐蚀')
subplot(223),imshow(dealImage2),title('原图-腐蚀')

效果图:

image-20220804144738253

傅里叶变换

图像的频谱图展示

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
clear all;close all;
rgb = imread('D:\reports\label.jpg');
image = rgb2gray(rgb);
F = fft2(im2double(image)); % FFT傅里叶变换
F = fftshift(F); % FFT频谱平移
F = real(F);
T = log(F+1); %频谱对数变换
subplot(231),imshow(image),title('原图像')
subplot(234),imshow(T,[]),title('原图像的频率图')
noisyImage = imnoise(image ,'salt & pepper',0.04); % 向图像之中添加椒盐噪声
noisyF = fft2(im2double(noisyImage));
noisyF = fftshift(noisyF);
noisyF = real(noisyF);
noisyT = log(noisyF+1); % 频谱对数变换
subplot(232),imshow(noisyImage),title('椒盐噪声')
subplot(235),imshow(noisyT,[]),title('椒盐噪声的频率图')

效果图:
image-20220804153219688

高斯低通(高通)频域滤波

低通滤波:对频谱图做高斯低通滤波,使低频通过而使高频衰减,滤波之后的图像会变得模糊,比原始图像减少尖锐的细节部分而突出平滑过渡部分,

高通滤波:使高频通过而使低频衰减,结果发现滤波后的图像变锐化,比原始图像减少平滑过渡而突出边缘等细节部分。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
image = imread('D:\reports\blackAndWhite.jpg')
I = rgb2gray(image)
P = double(I);
M = im2double(I); % 转化为归一化二维矩阵
Q = fft2(P); % 转化为二维矩阵
N = fft2(M);
move1 = fftshift(N); % 转换到频谱的中心
move2 = fftshift(Q);
subplot(231),imshow(I),title('原图像')
HW = size(move2);
d0 = 50 ;% 截至频率,数值越小越平滑
row1 = fix(HW(1)/2)
col1 = fix(HW(2)/2)
for i = 1:HW(1)
for j = 1:HW(2)
d = sqrt((i-row1)^2+(j-col1)^2);
hl(i,j) = exp(-d^2 /(2*d0^2));% 高斯低通滤波器
hh(i,j) = 1-hl(i,j); % 高斯高通滤波器
gl(i,j) = hl(i,j) * move2(i,j); % 高斯低通滤波
gh(i,j) = hh(i,j) * move2(i,j); % 高斯高通滤波
end
end
% 对图像进行傅里叶的反变换
gl = ifftshift(gl);
g1 = ifft2(gl);
k1 = uint8(real(gl));
gh = ifftshift(ifft(gh));
k2 = uint8(real(gh));
subplot(232),imshow(k1),title('原图像低通滤波')
subplot(233),imshow(k2),title('原图像高通滤波')

效果图:
image-20220804165900350

总结:不知道为什么达不到平滑图像的效果,对比网上的代码,感觉也没有错

参考文章

频率域波的运用处理

图像的合成处理

对图像数据进行加减异或运算从而对图像进行处理,这里就简单的展示了一下加运算

加运算

数据信息:uint8 类型的数据容易溢出,饱和截取操作;或者预先转化为uint16数据类型。

代码:

1
2
3
4
5
6
7
8
9
10
image = imread('D:\reports\blackAndWhite.jpg');
image = rgb2gray(image);
image2 = imread('D:\reports\test.jpg');
image2 = rgb2gray(image2);
newImage = imadd(image,image2,'uint8');
newImage2 = imadd(image,image2,'uint16');
subplot(231),imshow(image),title('图像1')
subplot(232),imshow(image2),title('图像2')
subplot(234),imshow(newImage,[]),title('相加(饱和截取)')
subplot(235),imshow(newImage2,[]),title('uint16数据类型')

效果:

image-20220804172555954

彩色图像的处理

cat运算符可以堆叠出RGB图像

rgb2ind()函数生成8色的抖动和非抖动图像

代码:

1
2
3
4
5
6
image = imread('D:\reports\test.jpg');
image2 = imread('D:\reports\test.jpg');
[f1,map] = rgb2ind(image,8,'nodither');% 将rgb图像转化为索引图像,且索引图像是8种颜色(不抖动)
subplot(121),imshow(f1,map),title('rgb2ind无抖动图像');
[f2,map2] = rgb2ind(image,8,'dither');
subplot(122),imshow(f2,map2),title('rgb2ind抖动图像');

效果图:

image-20220804234124105

‘replicate’:图像大小通过复制外边界的值来扩展

平滑rgb图像

代码:

1
2
3
4
5
6
image = imread('D:\reports\test.jpg');
image2 = imread('D:\reports\test.jpg');
mode = fspecial('motion',50,40); % 创建一个滤波
filteredImage = imfilter(image,mode);
subplot(121),imshow(image),title('原图像')
subplot(122),imshow(filteredImage),title('平滑图像')

效果图:

image-20220804235002633

锐化rgb图像

这里是用原图像减去拉普拉斯算子处理之后的图像得到锐化图像

代码:

1
2
3
4
5
image = imread('D:\reports\test.jpg');
mode = [1 1 1 1 1;1 1 1 1 1;1 1 -24 1 1;1 1 1 1 1;1 1 1 1 1];% 拉普拉斯算子
filteredImage = image - imfilter(image,mode,'replicate');
subplot(121),imshow(image),title('原图像')
subplot(122),imshow(filteredImage),title('锐化图像')

效果图:

image-20220804235923711

将图像转化为hsi(色调,饱和度,亮度),再对其中一个分量进行处理得到的图像

rgb2hsi()函数代码(同目录之下创建rgb2hsi文件)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function hsi = rgb2hsi(rgb)
%hsi = rgb2hsi( rgb )将rgb转化为hsi
%输出hsi是double型的

% rgb =imread(rgb);%没有这句会报错
%提取图像RGB分量
rgb = im2double(rgb);
r = rgb(:,:,1);
g = rgb(:,:,2);
b = rgb(:,:,3);

%执行转化
num = 0.5*((r-g)+(r-b));
den = sqrt((r-g).^2 +(r-b).*(g-b));
theta = acos(num./(den + eps)); %eps极小值标示,防止除以0出错

H = theta;
H( b>g ) = 2*pi - H( b>g );
H = H/(2*pi);

num = min(min(r,g),b);
den = r + g + b;
den(den == 0) = eps;%eps极小值标示,防止除以0出错
S = 1 - 3.* num./den;
H( S==0 ) = 0;
I = ( r + g + b )/3;

%组合hsi图像,H、S、I是矩阵
hsi = cat(3,H,S,I);


end

hsi2rgb()函数脚本:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
%hsi2rgb源程序
function rgb=hsi2rgb(hsi)
%rgb=hsi2rgb(hsi)把一幅HSI图像转换为RGB图像
%其中HSI是double型
%hsi(:,:,1)是色度分量,它的范围是除以2*pi后在[0,1]内
%hsi(:,:,2)是饱和度分量,它的范围是在[0,1]内
%hsi(:,:,3)是亮度分量,它的范围是在[0,1]内

%输出图像的分量是
%rgb(:,:,1)为红
%rgb(:,:,2)为绿
%rgb(:,:,3)为蓝

%提取HSI的各个分量
hsi=im2double(hsi);%把hsi转化为双精度浮点类型
H=hsi(:,:,1)*2*pi;
S=hsi(:,:,2);
I=hsi(:,:,3);

%执行变换方程
R=zeros(size(hsi,1),size(hsi,2));
G=zeros(size(hsi,1),size(hsi,2));
B=zeros(size(hsi,1),size(hsi,2));

%RG区(0<=H<2*pi/3)
idx=find((0<=H) & (H<2*pi/3));%寻找0<=H<2*pi/3
B(idx)=I(idx).*(1-S(idx));
R(idx)=I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx)));
G(idx)=3*I(idx)-(R(idx)+B(idx));

%BG区(2*pi/3<=H<4*pi/3)
idx=find((2*pi/3<=H) & (H<4*pi/3));%寻找2*pi/3<=H<4*pi/3
R(idx)=I(idx).*(1-S(idx));
G(idx)=I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)./cos(pi-H(idx)));
B(idx)=3*I(idx)-(R(idx)+G(idx));

%BR区(4*pi/3<=H<=2*pi)
idx=find((4*pi/3<=H) & (H<=2*pi));%寻找4*pi/3<=H<=2*pi
G(idx)=I(idx).*(1-S(idx));
B(idx)=I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)./cos(5*pi/3-H(idx)));
R(idx)=3*I(idx)-(G(idx)+B(idx));

%将3个分量联合成为一个RGB图像
rgb=cat(3,R,G,B);
rgb=max(min(rgb,1),0);
end

主函数(这里是对亮度分量进行平滑处理的结果)

1
2
3
4
5
6
7
8
9
10
11
12
image = imread('D:\reports\test.jpg');
hsiImage = rgb2hsi(image)% 将原图像转化为hsi图像
% 分别得到hsi的三个分量
h = hsiImage(:,:,1);
s = hsiImage(:,:,2);
i = hsiImage(:,:,3);
mode = ones(25)./(25*25);
filterI = imfilter(i,mode,'replicate');% 对其中的i分量,亮度分量使用filter进行平滑处理
hsiImage = cat(3,h,s,filterI);% 将这三个方向的值组合在一起
image2 = hsi2rgb(hsiImage)% 将hsi图像转化成rgb图像
subplot(121),imshow(image),title('原图像')
subplot(122),imshow(image2),title('亮度分量平滑图像')

效果图:

image-20220805002223940

参考文章

小波变换

小波变换有着窗口自适应的特点,即高频信号分辨率高(但是频率分辨率差),低频信号频率分辨率高(但是时间分辨率差)

wavemngr(‘read’,1):进行小波的查看

waveinfo():查看小波的相关信息

小波在时间上的移动,逐一比较不同位置的窗口信号,小波系数,小波系数越大,则证明小波与该段信号的拟合程度越好,

小波的信息查询:waveinfo函数,例如waveinfo(‘wname’)

wname代表的小波有:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
'haar' : Haar wavelet.
'db' : Daubechies wavelets.
'sym' : Symlets.
'coif' : Coiflets.
'bior' : Biorthogonal wavelets.
'rbio' : Reverse biorthogonal wavelets.
'meyr' : Meyer wavelet.
'dmey' : Discrete Meyer wavelet.
'gaus' : Gaussian wavelets.
'mexh' : Mexican hat wavelet.
'morl' : Morlet wavelet.
'cgau' : Complex Gaussian wavelets.
'cmor' : Complex Morlet wavelets.
'shan' : Complex Shannon wavelets.
'fbsp' : Complex Frequency B-spline wavelets.
'fk' : Fejer-Korovkin orthogonal wavelets

image-20220805103058654

wfilters函数

[LOD,HID,LOR,HIR] = wfilters(‘wname’):计算正交小波或双正交小波(wname)有关联的四个滤波器

LO_D,分解低通滤波器;HI_D,分解高通滤波器;LO_R,重构低通滤波器;HI_R,重构高通滤波器

代码:

1
2
3
4
5
6
7
8
9
10
11
close all;
[LOD,HID,LOR,HIR] = wfilters('db45');% 得到该小波关联的四个滤波器
% 将这四个滤波器显示出来
subplot(221),stem(LOD),xlim([0,100]),title('分解低通滤波器')
xlabel('x'),ylabel('y')
subplot(222),stem(HID),xlim([0,100]),title('分解高通滤波器')
xlabel('x'),ylabel('y')
subplot(223),stem(LOR),xlim([0,100]),title('重构低通滤波器')
xlabel('x'),ylabel('y')
subplot(224),stem(HIR),xlim([0,100]),title('重构低通滤波器')
xlabel('x'),ylabel('y')

效果图:

image-20220805105124689

二维离散小波变换图像单层小波分解

使用dwt2()函数

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
close all;
clear all;
image = imread('D:\reports\test.jpg');% 获取照片
T = rgb2gray(image);% 将照片转化成灰度模式
[sm1,deHo,deVe,deDi] = dwt2(T,'bior3.7');% 将灰度照片单层分解成四个层面
% 将四个层面的图像分别输出
figure
subplot(141),imshow(uint8(sm1),title('近似分量')% 近似分量
subplot(142),imshow(deHo),title('细节水平分量')% 细节水平分量
subplot(143),imshow(deVe),title('细节垂直分量')% 细节垂直分量
subplot(144),imshow(deDi),title('细节对角分量')% 细节对角分量'
figure
subplot(121),imshow(image),title('原图')% 原图
subplot(122),imshow([sm1,deHo;deVe,deDi]),title('小波变换分量组合')% 小波变换分量组合图像

效果图:

image-20220805131621875

小波实现图像的重构

idwt2()函数

这里使用db1小波函数对图像进行分解,然后再用db4的小波函数对图像进行重构,做差之后就能直观的看到前后图像的差别

代码:

1
2
3
4
5
6
7
8
9
10
11
12
close all;
clear all;
image = imread('D:\reports\test.jpg');% 获取照片
T = rgb2gray(image);% 将照片转化成灰度模式
[sm1,deHo,deVe,deDi] = dwt2(T,'db1');% 将灰度照片单层分解成四个层面
imageS = size(T);%获得原图像的尺度
T0 = idwt2(sm1,deHo,deVe,deDi,'db4',imageS);% 利用小波分解的第一层系数去重构
T1 = double(T)-T0;
figure
subplot(131),imshow(uint8(T)),title('原图像')
subplot(132),imshow(uint8(T0)),title('重构图像')
subplot(133),imshow(uint8(T1)),title('差异图像')

效果图:

image-20220805132501969

多层小波重构

使用函数upcoef2()完成重构

appcoef2():获取小波分解的近似系数

[h,v,d] = detcoef2(‘all’,c,s,n):返回级别n处的水平h、垂直v和对角线d细节系数。’all‘这个参数也可以相应的修改而获得三个方向对应的细节系数

二层的分解需要将用小波函数分解出两层,然后再对每一层再进行小波变换的分解,重构的时候分方向重构。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
close all;
clear all;
image = imread('D:\reports\test.jpg');% 获取照片
T = rgb2gray(image);% 将照片转化成灰度模式
[c,s] = wavedec2(T,2,'db4'); % 对图像进行小波的二层分解 wavedec2 小波函数使用db4
siz = s(size(s,1),:);% 获得第二层小波分解的系数矩阵
ca1 = appcoef2(c,s,'db4',2); % 小波分解的近似系数 appcoef()
a1 = upcoef2('a',ca1,'db4',2,siz); % upcoef2()函数进行重构

chd1 = detcoef2('h',c,s,2); % 小波分解的细节系数水平分量detcoef2()'h
hd1 = upcoef2('h',chd1,'db4',2,siz);

cvd1 = detcoef2('v',c,s,2); % 细节系数垂直分量'v'
vd1 = upcoef2('v',cvd1,'db4',2,siz);

cdd1 = detcoef2('d',c,s,2); % 细节系数对角分量'd'
dd1 = upcoef2('d',cdd1,'db4',2,siz);

result1 = a1+hd1+vd1+dd1; % 将四个重构的分量相加得到图像的重构结果
figure
subplot(221),imshow(uint8(a1)),title('小波分解的近似系数')
subplot(222),imshow(hd1),title('小波分解的近似系数水平分量')
subplot(223),imshow(vd1),title('小波分解的细节系数垂直分量')
subplot(224),imshow(dd1),title('小波分解的细节系数对角分量')
figure
subplot(121),imshow(uint8(T)),title('原图像')
subplot(122),imshow(uint8(result1)),title('二层分解重构图像')

效果图:

image-20220805141517643

小波去噪去噪

wrcoef2():使用wname指定的小波返回基于图像的小波分解结构[c,s]的type重构系数矩阵

例:x = wrcoef2(type,c,s,wname) ;a2 = wrcoef2(“a”,c,s,”sym5”,2);

NC = wthcoef2(‘type’,C,S,N,T,SORH) :向量N和T中定义阈值,再从小波分解结构[C,S]中获得的水平、垂直或对角系数

参考文章

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
close all;
clear all;
image = imread('D:\reports\test.jpg');% 获取照片
T = rgb2gray(image);% 将照片转化成灰度模式
T = imnoise(T,'speckle',0.1);
subplot(131),imshow(uint8(T)),title('加噪图像')
% init=2055615866;
% randn('seed',init);
% T = double(T) + 2*randn(size(T));
[c,s] = wavedec2(T,2,'sym4'); % 对图像进行小波的二层分解 wavedec2 小波函数使用db4
result1 = wrcoef2('a',c,s,'sym4',2);% 重构第二层的近似系数,小波分解去噪的操作
subplot(132),imshow(result1),title('小波分解去噪')
n = [1,2];% 设置阈值N 尺度向量
r = [10.28, 24.08];% 设置阈值T 阈值向量
dealData0 = wthcoef2('t',c,s,n,r,'s');% 对高频小波进行阈值处理 wthcoef2
dealData1 = wthcoef2('t',dealData0,s,n,r,'s');% 对高频小波进行第二次阈值处理
result2 = waverec2(dealData1,s,'sym4');% 对阈值处理之后的小波进行重构 waverec2()
subplot(133),imshow(result1),title('小波阈值去噪')

实验结果:
image-20220805154151408

总结:能够看到噪声明显减弱了,但是图像也更模糊了,去噪的过程就是过滤高频小波,高频小波之中表现的是细节