YDYの博客

一只有理想的菜鸟

数字图像处理复习

发表于 分类于

数字图像处理复习

im2uint8

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%% im2uint8
clc
clear
f1 = [-0.5 0.5
     0.75 1.5]
g1 = im2uint8(f1)

f2 = uint8(f1*100)
g2 = im2uint8(f2)

f3 = f1*100
g3 = im2uint8(f3)

f4 = uint16(f1*50000)
g4 = im2uint8(f4)

灰度图象变为二值图像

1
2
3
4
5
6
7
8
9
clc
clear

I = imread('liftingbody.png');
imshow(I)

BW = im2bw(I,0.46);

figure,imshow(BW)

计算两幅图像间的绝对差

1
2
3
4
5
6
7
8
9
clc
clear

I = imread('cameraman.tif')
figure,imshow(I,[])
J = uint8(filter2(fspecial('gaussian'), I));
figure,imshow(J,[])
K = imabsdiff(I,J);
figure,imshow(K,[]) % [] = scale data automatically

对图像求补

1
2
3
4
5
6
7
8
%% imcomplement 对图像求补
clc
clear

bw = imread('text.png');
bw2 = imcomplement(bw)
subplot(1,2,1),imshow(bw)
subplot(1,2,2),imshow(bw2)

im2double

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
clc
clear
f1 = [-0.5 0.5
     0.75 1.5]
g1 = im2double(f1)

f22 = uint8(f1)
f2 = uint8(f1*100)
g2 = im2double(f2)

f3 = f1*100
g3 = im2double(f3)

f4 = uint16(f1*50000)
g4 = im2double(f4)

imadjust

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
g = imread('Fig0222(a)(face).tif')
imshow(g);

g1 = imadjust(g,[0 1],[1 0],1);% 反转
figure,imshow(g1);

g2 = imcomplement(g); 
figure,imshow(g2);

g2 = imadjust(g,[0.5 0.75],[0 1]);
g3 = imadjust(g,[ ],[ ],2);
figure,imshow(g3);

g4 = imadjust(g,stretchlim(g),[]);
figure,imshow(g4);

1
2
3
4
5
6
7
8
9
10
11
12
13
%% 使用函数 imadjust 目的：突出我们感兴趣的亮度带·压缩灰度级的低端并扩展灰度级的高端
clc
clear
f = imread('./CH2/Fig0222(a)(face).tif');
imshow(f)

g1 = imadjust(f,[0,1],[1,0]); % === imcomplement(f) 灰度反转 @ 灰度负片
imshow(g1)

g2 = imadjust(f,[0.5,0.75],[0,1]); % 突出我们感兴趣的亮度带
imshow(g2)

g3 = imadjust(f,[],[],2); % 压缩灰度级的低端并扩展灰度级的高端

分辨率

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
imagePath = '1.tif';
img = imread(imagePath);
%imshow(img)
% 读取图像信息
info = imfinfo(imagePath);
% 获取原始图片的分辨率（DPI）
xDPI = info.XResolution;
yDPI = info.YResolution;
% 打印原始图片的DPI
fprintf('原始图片的分辨率：\n');
fprintf('水平分辨率：%.2f DPI\n', xDPI);
fprintf('垂直分辨率：%.2f DPI\n', yDPI);

figure(1)
image_A=imresize(img,0.05);
imshow(image_A)
title('75dpi')

figure(2)
image_B=imresize(img,0.1);
imshow(image_B)
title('125dpi')

figure(3)
image_C=imresize(img,0.2);
imshow(image_C)
title('250dpi')

插值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
imagePath = '1.tif';
img = imread(imagePath);
%imshow(img)

image_A=imresize(img,0.1);
figure(1)
imshow(image_A)
title('125dpi')

figure(2)
image_B=imresize(img,0.2);
imshow(image_B)
title('理想：250dpi')

figure(3)
image_C=imresize(image_A,2,'nearest');
imshow(image_C)
title('最近邻插值：250dpi')

figure(4)
image_D=imresize(image_A,2,'bilinear');
imshow(image_D)
title('双线性插值：250dpi')

直方图均衡

1
2
3
4
5
6
7
8
9
10
I=imread('Fig0316(1).tif');
figure(1);
imshow(I);
figure(2);
imhist(I) % 计算和显示图象直方图
leq=histeq(I,256); % 直方图均衡化处理
figure(3);
imshow(leq);
figure(4);
imhist(leq)
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
I=imread('Fig0316(4).tif');
imshow(I);
[M,N]=size(I);
Id=double(I);
% 遍历所有像素，统计图像灰度直方图
IHist=zeros(1,256);
for i=1:M
    for j=1:N
         IHist(Id(i,j)+1)=IHist(Id(i,j)+1)+1;
    end
end
figure,plot(IHist);
%直方图归一化处理
IHist=IHist./(M*N);

% 采用输入图像概率累积函数进行灰度级映射计算
Sk=zeros(1,256);
for k=0:255
  Sk(k+1)=sum(IHist(1:k+1));
end
figure,plot(Sk)
% 灰度级的重新量化
Smin=min(Sk)
Sk=uint8(255*(Sk-Smin)./(1-Smin)+0.5);
figure,plot(Sk)

%输出经直方图均衡化的图像
Ieq=zeros(M,N);
for i=1:M
    for j=1:N
        L=double(I(i,j))+1;
        Ieq(i,j)=Sk(L);       
    end
end
Ieq=mat2gray(Ieq);
figure,imshow( Ieq);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%% 例3.5 直方图均衡化 默认为64
clc
clear
f = imread('./CH2/Fig0222(a)(face).tif');
subplot(121),imshow(f),subplot(122),imhist(f)
ylim('auto')

g = histeq(f,256);
figure,subplot(121),imshow(g),subplot(122),imhist(g)
ylim('auto')


g = histeq(f,128);
figure,subplot(121),imshow(g),subplot(122),imhist(g)
ylim('auto')


g = histeq(f);  % 默认为64
figure,subplot(121),imshow(g),subplot(122),imhist(g)
ylim('auto')

g = histeq(f,8);
figure,subplot(121),imshow(g),subplot(122),imhist(g)
ylim('auto')

对数变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
clc
clear

f = imread('./CH2/Fig0222(a)(face).tif');
subplot(121),imshow(f),subplot(122),imhist(f),axis tight


g = im2uint8(mat2gray(log(1 + double(f))));
figure,subplot(121),imshow(g),subplot(122),imhist(g),axis tight
title('使用对数变换减小动态范围')
% axis([0 255 0 4000])

% 对比度拉伸变换
m = 5;
E = 10;
h = im2uint8(mat2gray(1./(1 + (m./(double(f) + eps)).^E)));
figure,subplot(121),imshow(h),subplot(122),imhist(h),axis tight

线性滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%% imfilter 线性空间滤波（空间卷积）
clc
clear
f = imread('./CH2/Fig0222(a)(face).tif')
f = im2double(f);
figure,imshow(f)

w = ones(31);

gd = imfilter(f, w);
figure,imshow(gd,[])


gr = imfilter(f, w, 'replicate'); % 边界处理方法为复制边界值
figure,imshow(gr,[])

gc = imfilter(f, w, 'circular');
figure,imshow(gc,[])

f8 = im2uint8(f);
gr8 = imfilter(f8, w, 'replicate');
figure,imshow(gr8,[])

均值滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%使用自带函数
Picture=imread('fig3.tif');
subplot(1,4,1);imshow(Picture,[]);title('原图');
Tem_Smooth=fspecial('average', [3, 3]);
P_Smoothl=imfilter(Picture,Tem_Smooth);
subplot(1,4,2);imshow(P_Smoothl,[]);title('3*3均值滤波一次后图像');
%P_Smooth2=imfilter(P_Smoothl, Tem_Smooth);
%subplot(1,4,3);imshow(P_Smooth2,[]);title('均值滤波二次后图像');
%P_Smooth3=imfilter(P_Smooth2,Tem_Smooth);
%subplot(1,4,4);imshow(P_Smooth3,[]);title('均值滤波三次后图像');
Tem_Smooth=fspecial('average',[5,5]);
P_Smooth2=imfilter(Picture,Tem_Smooth);
subplot(1,4,3);imshow(P_Smooth2,[]);title('5*5均值滤波后图像');
Tem_Smooth=fspecial('average',[15,15]);
P_Smooth3=imfilter(Picture,Tem_Smooth);
subplot(1,4,4);imshow(P_Smooth3,[]);title('15*15均值滤波后图像');

中值滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f=imread( 'img\Fig0318.tif' );

% 非线性空间滤波器
fn=imnoise(f, 'salt & pepper', 0.2); % 椒盐噪声
gm=medfilt2(fn); % 中值滤波
gms=medfilt2(fn, 'symmetric');

subplot(2,2,1)
imshow(f, [ ])

subplot(2,2,2)
imshow(fn, [ ])

subplot(2,2,3)
imshow(gm, [])

subplot(2,2,4)
imshow(gms, [ ])

锐化

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64

%   骨骼图像增强  

im=imread('fig4.tif');  
im = im2double(im);  
%im=rgb2gray(im);  
  
%---原始图像  
subplot(2,4,1);   
imshow(im);  
title('1：原始图像');  
  
%---为图2，使用模板为[-1,-1,-1;-1,8,-1;-1,-1,-1]的滤波器对原图像进行拉普拉斯操作，为了便于显示，对图像进行了标定，这一步先对图像进行初步的锐化滤波。  
subplot(2,4,2);   
h =[-1,-1,-1;-1,8,-1;-1,-1,-1];                                    
im1 =imfilter(im,h);   
imshow(im1);   
title('2：拉普拉斯操作后图像');  
  
%---图3，由于使用的模板如上，让常数c=1，简单的将原图和图2相加就可以得到一幅经过锐化过的图像。  
 subplot(2,4,3);  
 im2=im+im1;  
 imshow(im2)  
 title('3:1图和2图相加后图像');  
  
%---图4，对原图像试用Sobel梯度操作，分量gx为[-1,-2,-1;0,0,0;1,2,1],而分量gy为[-1,0,1;-2,0,2;-1,0,1]的模板。  
  subplot(2,4,4);  
  hx=[-1,-2,-1;0,0,0;1,2,1]; %生产sobel垂直梯度模板  
  hy=[-1,0,1;-2,0,2;-1,0,1]; %生产sobel水平梯度模板  
  gradx=filter2(hx,im,'same');  
  gradx=abs(gradx); %计算图像的sobel垂直梯度  
  grady=filter2(hy,im,'same');  
  grady=abs(grady); %计算图像的sobel水平梯度  
  im3=gradx+grady;  %得到图像的sobel梯度  
  imshow(im3,[]);  
  title('4:1图sobel梯度处理后图像');  
  
% ---图5，使用大小为5*5的一个均值滤波器得到平滑后的Sobel梯度图像。  
subplot(2,4,5);  
h1 = fspecial('average',5) ;                                                       
im4 = imfilter(im3,h1);   
imshow(im4);  
title('5:使用5*5均值滤波器平滑后的sobel图像');  
  
% --图6，将拉普拉斯图像（即图3）与平滑后的梯度图像（即图5）进行点乘。  
subplot(2,4,6);  
% im5=immultiply(im2,im4);  
im5=im2.*im4;  
imshow(im5);  
title('6:3图和5图相乘相乘的掩蔽图像');  
  
% --图7，将乘积图像（即图6）与原图像相加就产生一幅需要的锐化图像。  
subplot(2,4,7);  
im6=im+im5;  
imshow(im6);  
title('7:1图和6图求和得到的锐化图像');  
  
% --图8，扩展灰度范围，对图7进行幂率变换处理，r=0.5，c=1，然后即可对图像进行幂率变换  
subplot(2,4,8);  
gamma=0.5;  
c=1;  
im7=c.*im6.^gamma;  
imshow(im7);  
title('8:图7进行幂率变换后的最终图像');

低通滤波

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
51
52
53
54
55
56
57
f=imread( 'img\Fig0405.tif' );
imshow(f, [])


% 获取图像的大小（行数 M 和列数 N）
[M, N] = size(f);

% 对图像进行二维快速傅里叶变换 (FFT)
F = fft2(f);

% 设置高斯低通滤波器的标准差
sig = 10;

% 生成高斯低通滤波器
H = lpfilter('gaussian', M, N, sig);

% 在频域中应用滤波器 H
G = H.*F;

% 对滤波后的图像进行逆傅里叶变换，并取其实部
g = real(ifft2(G));

% 显示滤波后的图像 g
figure;
imshow(g, [])

% 计算图像适用于频域滤波的填充尺寸
PQ = paddedsize(size(f));

% 对填充后的图像进行二维快速傅里叶变换
Fp = fft2(f, PQ(1), PQ(2));

% 生成填充尺寸的高斯低通滤波器，标准差为 sig 的两倍
Hp = lpfilter('gaussian', PQ(1), PQ(2), 2*sig);

% 在频域中应用填充后的滤波器 Hp
Gp = Hp.*Fp;

% 对滤波后的填充图像进行逆傅里叶变换，并取其实部
gp = real(ifft2(Gp));

% 截取图像恢复原尺寸
gpc = gp(1:size(f,1), 1:size(f,2));

% 显示填充滤波后的图像 gp
figure;
imshow(gp, [])

% 生成大小为 15x15，标准差为 7 的空间域高斯滤波器
h = fspecial('gaussian', 15, 7);

% 应用空间域高斯滤波器对图像进行卷积
gs = imfilter(f, h);

% 显示空间域高斯滤波后的图像 gs
figure;
imshow(gs, [])

高斯低通滤波器

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
f=imread( 'img\Fig0413.tif' );
imshow(f, [])
% 计算图像适用于频域滤波的填充尺寸 PQ
PQ = paddedsize(size(f));

% 生成用于频域滤波器的坐标矩阵 U 和 V
[U, V] = dftuv(PQ(1), PQ(2));

% 定义高斯低通滤波器的截止频率 D0，为图像宽度的 5%
D0 = 0.05 * PQ(2);

% 对填充后的图像进行二维快速傅里叶变换 (FFT)
F = fft2(f, PQ(1), PQ(2));

% 生成高斯低通滤波器 H
H = exp(-(U.^2 + V.^2) / (2 * (D0^2)));

% 应用频域滤波器 H，并对结果进行逆傅里叶变换以恢复到空间域
g = dftfilt(f, H);

% 显示滤波器 H 的中心移动到图像中心后的结果
figure;
imshow(fftshift(H), [])

% 显示滤波器 H 的对数幅度谱，以便更清楚地看到低频和高频的分布
figure;
imshow(log(1 + abs(fftshift(H))), [])

% 显示滤波后的图像 g
figure;
imshow(g, [])

Sobel 滤波器

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
f=imread( 'img\Fig0409.tif' );
% 显示读取的图像 f
imshow(f, [])

% 对图像进行二维快速傅里叶变换 (FFT)
F = fft2(f);

% 计算傅里叶变换的频谱图，并进行对数变换和中心化
S = fftshift(log(1 + abs(F)));

% 归一化频谱图以便显示
S = gscale(S);

% 显示频谱图 S
figure;
imshow(S, [])

% 创建 Sobel 滤波器
h = fspecial('sobel')';

% 显示 Sobel 滤波器的频率响应
freqz2(h);

% 计算图像适用于频域滤波的填充尺寸 PQ
PQ = paddedsize(size(f));

% 计算 Sobel 滤波器在频域的频率响应 H
H = freqz2(h, PQ(1), PQ(2));

% 将频域滤波器的零频点移动到左上角
H1 = ifftshift(H);

% 显示频域滤波器 H 和 H1
imshow(abs(H), []);
figure;
imshow(abs(H1), [])

% 在空间域应用 Sobel 滤波器
gs = imfilter(double(f), h);

% 在频域应用 Sobel 滤波器
gf = dftfilt(f, H1);

% 显示空间域滤波后的图像 gs
figure;
imshow(gs, [])

% 显示频域滤波后的图像 gf
figure;
imshow(gf, [])

% 显示空间域滤波后的图像 gs 的绝对值
figure;
imshow(abs(gs), [])

% 显示频域滤波后的图像 gf 的绝对值
figure;
imshow(abs(gf), [])

% 显示空间域滤波后的二值化图像，阈值为最大值的 20%
figure;
imshow(abs(gs) > 0.2 * abs(max(gs(:))), [])

% 显示频域滤波后的二值化图像，阈值为最大值的 20%
figure;
imshow(abs(gf) > 0.2 * max(abs(gf(:))), [])

% 计算空间域和频域滤波结果的差异图像
d = abs(gs - gf);

% 输出差异图像的最大值和最小值
max(d(:))
min(d(:))

fftshift

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
%% fftshift 对数变换 fftshift 主要用来演示 H 使用
clc
clear
f = imread('./CH4/Fig0403(a)(image).tif');
figure,imshow(f)
title('原始图像')
F = fft2(f);
S = abs(F);
% S(1:5,1:5)
figure,imshow(S,[])
title('傅立叶频谱图像')


Fc = fftshift(F);
S = abs(Fc);
% S(1:5,1:5)
figure,imshow(S,[])
title('居中的傅立叶频谱图像')


S = abs(F);
% S(1:5,1:5)
S2 = log(1+S);
figure,imshow(S2,[])
title('使用对数变换进行视觉增强后的傅立叶频谱图像')


Fc = fftshift(F);
S = abs(Fc);
% S(1:5,1:5)
S2 = log(1+S);
figure,imshow(S2,[])
title('使用对数变换进行视觉增强并居中后的傅立叶频谱图像')

高斯滤波

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
f=imread( 'img\Fig0405.tif' );
imshow(f, [])

[M, N] = size(f);
F = fft2(f);
sig = 10;
H = lpfilter('gaussian', M, N, sig);
G = H.*F;
g = real(ifft2(G));
figure;
imshow(g, [])

PQ = paddedsize(size(f));
Fp = fft2(f,PQ(1), PQ(2));
Hp = lpfilter('gaussian', PQ(1), PQ(2),2*sig);
Gp = Hp.*Fp;
gp = real(ifft2(Gp));
gpc = gp(1:size(f,1), 1:size(f,2));
figure;
imshow(gp, [])

h = fspecial('gaussian', 15, 7);
gs = imfilter(f, h);
figure;
imshow(gs, [])

使用填充和不使用填充的滤波效果

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

clc
clear

f = imread('./CH4/Fig0405(a)(square_original).tif');
f = im2double(f);
figure,imshow(f,[])
title('原始图像')

[M, N] = size(f);
F = fft2(f);  % max(F(:)) = 128000
figure,imshow(F,[])
title('傅立叶频谱（复数）图像')

sig = 10;
H = lpfilter('gaussian',M,N,sig); % max(H(:)) = 1
figure,imshow(1-H,[]) % 显示表明滤波器图像未置中
title('滤波器频谱（取反）图像')

G = H.*F;
g = real(ifft2(G));
figure,imshow(g,[])
title('不使用填充的频域低通滤波处理后的图像')

PQ = paddedsize(size(f)); % size(f)=[256 256]
Fp = fft2(f, PQ(1),PQ(2));
Hp = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512 512]
figure,imshow(fftshift(Hp),[])
Gp = Hp.*Fp;
gp = real(ifft2(Gp));
figure,imshow(gp,[])
title('使用填充的频域低通滤波处理后的图像')

gpc = gp(1:size(f,1),1:size(f,2));
figure,imshow(gpc,[])
title('使用填充的频域低通滤波处理后的（截取原始大小）图像')


h = fspecial('gaussian',15,7);
gs = imfilter(f,h);
figure,imshow(gs,[])
title('使用空间滤波器处理后的图像')
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
function H = lpfilter(type, M, N, D0, n)
%LPFILTER Computes frequency domain lowpass filters.
%   H = LPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a lowpass filter, H, of the specified TYPE and size (M-by-N). To
%   view the filter as an image or mesh plot, it should be centered
%   using H = fftshift(H). 
%
%   Valid values for TYPE, D0, and n are:
%
%   'ideal'    Ideal lowpass filter with cutoff frequency D0. n need
%              not be supplied.  D0 must be positive.
%
%   'btw'      Butterworth lowpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0.  D0 must be
%              positive.
%
%   'gaussian' Gaussian lowpass filter with cutoff (standard
%              deviation) D0.  n need not be supplied.  D0 must be
%              positive. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.8 $  $Date: 2004/11/04 22:33:16 $

% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances. 
[U, V] = dftuv(M, N);

% Compute the distances D(U, V).
D = sqrt(U.^2 + V.^2);

% Begin filter computations.
switch type
case 'ideal'
   H = double(D <= D0);
case 'btw'
   if nargin == 4
      n = 1;	
   end
   H = 1./(1 + (D./D0).^(2*n));
case 'gaussian'
   H = exp(-(D.^2)./(2*(D0^2)));
otherwise
   error('Unknown filter type.')
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function [U, V] = dftuv(M, N)
%DFTUV Computes meshgrid frequency matrices.
%   [U, V] = DFTUV(M, N) computes meshgrid frequency matrices U and
%   V.  U and V are useful for computing frequency-domain filter
%   functions that can be used with DFTFILT.  U and V are both
%   M-by-N.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.3 $  $Date: 2003/04/16 22:30:34 $

% Set up range of variables.
u = 0:(M - 1);
v = 0:(N - 1);

% Compute the indices for use in meshgrid.
idx = find(u > M/2);
u(idx) = u(idx) - M;
idy = find(v > N/2);
v(idy) = v(idy) - N;

% Compute the meshgrid arrays.
[V, U] = meshgrid(v, u);
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
function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering. 
%   PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
%   computes the two-element size vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
%   vectors, computes the two-element size vector PQ.  The elements
%   of PQ are the smallest even integers greater than or equal to 
%   AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]). 
    
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $

if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 & ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.
   
   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
elseif nargin == 3
   m = max([AB CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
else 
   error('Wrong number of inputs.')
end

fftshift ifftshift

1
2
3
4
5
6
7
8
9
10
11
12
13
clc
clear

A = [2 0 0 1
     0 0 0 0
     0 0 0 0
     3 0 0 4]

fftshift(A)

fftshift(fftshift(A))

ifftshift(fftshift(A))

高通滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clc
clear

f = imread('./CH4/Fig0413(a)(original_test_pattern).tif');
figure,imshow(f)
title('原始图像')

PQ = paddedsize(size(f));

D0 = 0.05*PQ(1); % 半径是 D0
H = hpfilter('gaussian',PQ(1),PQ(2),D0);
g = dftfilt(f,H);
figure,imshow(g,[])
title('高斯高通滤波后的图像')
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
function H = hpfilter(type, M, N, D0, n)
%HPFILTER Computes frequency domain highpass filters.
%   H = HPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a highpass filter, H, of the specified TYPE and size (M-by-N).
%   Valid values for TYPE, D0, and n are: 
%
%   'ideal'    Ideal highpass filter with cutoff frequency D0.  n
%              need not be supplied. D0 must be positive.
%
%   'btw'      Butterworth highpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0. D0 must be
%              positive.
%
%   'gaussian' Gaussian highpass filter with cutoff (standard
%              deviation) D0. n need not be supplied. D0 must be
%              positive.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/08/25 14:28:22 $

% The transfer function Hhp of a highpass filter is 1 - Hlp, 
% where Hlp is the transfer function of the corresponding lowpass 
% filter.  Thus, we can use function lpfilter to generate highpass 
% filters.

if nargin == 4
   n = 1; % Default value of n.
end

% Generate highpass filter.
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function g = dftfilt(f, H)
%DFTFILT Performs frequency domain filtering.
%   G = DFTFILT(F, H) filters F in the frequency domain using the
%   filter transfer function H. The output, G, is the filtered
%   image, which has the same size as F.  DFTFILT automatically pads
%   F to be the same size as H.  Function PADDEDSIZE can be used to
%   determine an appropriate size for H.
%
%   DFTFILT assumes that F is real and that H is a real, uncentered
%   circularly-symmetric filter function. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $

% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));

% Perform filtering. 
g = real(ifft2(H.*F));

% Crop to original size.
g = g(1:size(f, 1), 1:size(f, 2));

将高频滤波和直方图均衡化结合起来

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
clc
clear
f = imread('./CH4/Fig0419(a)(chestXray_original).tif');
figure,imshow(f)
title('原始图像')


PQ = paddedsize(size(f));
D0 = 0.05*PQ(1);
HBW = hpfilter('btw',PQ(1),PQ(2),D0,2);
H = 0.5+2*HBW;
gbw = dftfilt(f,HBW);
% 使用了 gscale(gbw) 之后，imshowMy(gbw) 等价于 imshowMy(gbw,[])
gbw = gscale(gbw); 
figure,imshow(gbw,[])
title('高通滤波后的图像')


gbe = histeq(gbw,256);
figure,imshow(gbe,[])
title('高通滤波并经过直方图均衡化后的图像')


ghf = dftfilt(f,H);
ghf = gscale(ghf);
figure,imshow(ghf,[])
title('高频强调滤波后的图像')


ghe = histeq(ghf,256);
figure,imshow(ghe,[])
title('高频强调滤波并经过直方图均衡化后的图像')
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
51
52
53
54
55
56
57
58
function g = gscale(f, varargin)
%GSCALE Scales the intensity of the input image.
%   G = GSCALE(F, 'full8') scales the intensities of F to the full
%   8-bit intensity range [0, 255].  This is the default if there is
%   only one input argument.
%
%   G = GSCALE(F, 'full16') scales the intensities of F to the full
%   16-bit intensity range [0, 65535]. 
%
%   G = GSCALE(F, 'minmax', LOW, HIGH) scales the intensities of F to
%   the range [LOW, HIGH]. These values must be provided, and they
%   must be in the range [0, 1], independently of the class of the
%   input. GSCALE performs any necessary scaling. If the input is of
%   class double, and its values are not in the range [0, 1], then
%   GSCALE scales it to this range before processing.
%
%   The class of the output is the same as the class of the input.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:36:09 $

if length(varargin) == 0 % If only one argument it must be f.
   method = 'full8';
else 
   method = varargin{1};
end

if strcmp(class(f), 'double') & (max(f(:)) > 1 | min(f(:)) < 0)
   f = mat2gray(f);
end

% Perform the specified scaling.
switch method
case 'full8'
   g = im2uint8(mat2gray(double(f)));
case 'full16'
   g = im2uint16(mat2gray(double(f)));
case 'minmax'
   low = varargin{2}; high = varargin{3};
   if low > 1 | low < 0 | high > 1 | high < 0
      error('Parameters low and high must be in the range [0, 1].')
   end
   if strcmp(class(f), 'double')
      low_in = min(f(:));
      high_in = max(f(:));
   elseif strcmp(class(f), 'uint8')
      low_in = double(min(f(:)))./255;
      high_in = double(max(f(:)))./255;
   elseif strcmp(class(f), 'uint16')
      low_in = double(min(f(:)))./65535;
      high_in = double(max(f(:)))./65535;    
   end
   % imadjust automatically matches the class of the input.
   g = imadjust(f, [low_in high_in], [low high]);   
otherwise
   error('Unknown method.')
end

各种噪声

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
clc
clear

r = imnoise2('gaussian',100000,1,0,1);
bins = 100;
hist(r,bins)
title('gaussian')

r = imnoise2('uniform',100000,1,0,1);
bins = 100;
figure,hist(r,bins)
title('uniform')

r = imnoise2('salt & pepper',1000,1,0.1,0.27);
bins = 100;
figure,hist(r,bins)
title('salt & pepper')

r = imnoise2('lognormal',100000,1);
bins = 100;
figure,hist(r,bins)
title('lognormal')

r = imnoise2('rayleigh',100000,1,0,1);
bins = 100;
figure,hist(r,bins)
title('rayleigh')

r = imnoise2('exponential',100000,1);
bins = 100;
figure,hist(r,bins)
title('exponential')

r = imnoise2('erlang',100000,1);
bins = 100;
figure,hist(r,bins)
title('erlang')
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
function R = imnoise2(type, M, N, a, b)
%IMNOISE2 Generates an array of random numbers with specified PDF.
%   R = IMNOISE2(TYPE, M, N, A, B) generates an array, R, of size
%   M-by-N, whose elements are random numbers of the specified TYPE
%   with parameters A and B. If only TYPE is included in the
%   input argument list, a  single random number of the specified
%   TYPE and default parameters shown below is generated. If only
%   TYPE, M, and N are provided, the default parameters shown below
%   are used.  If M = N = 1, IMNOISE2 generates a single random
%   number of the specified TYPE and parameters A and B.
%
%   Valid values for TYPE and parameters A and B are:
% 
%   'uniform'       Uniform random numbers in the interval (A, B).
%                   The default values are (0, 1).  
%   'gaussian'      Gaussian random numbers with mean A and standard
%                   deviation B. The default values are A = 0, B = 1.
%   'salt & pepper' Salt and pepper numbers of amplitude 0 with
%                   probability Pa = A, and amplitude 1 with
%                   probability Pb = B. The default values are Pa =
%                   Pb = A = B = 0.05.  Note that the noise has
%                   values 0 (with probability Pa = A) and 1 (with
%                   probability Pb = B), so scaling is necessary if
%                   values other than 0 and 1 are required. The noise
%                   matrix R is assigned three values. If R(x, y) =
%                   0, the noise at (x, y) is pepper (black).  If
%                   R(x, y) = 1, the noise at (x, y) is salt
%                   (white). If R(x, y) = 0.5, there is no noise
%                   assigned to coordinates (x, y). 
%   'lognormal'     Lognormal numbers with offset A and shape
%                   parameter B. The defaults are A = 1 and B =
%                   0.25.
%   'rayleigh'      Rayleigh noise with parameters A and B. The
%                   default values are A = 0 and B = 1. 
%   'exponential'   Exponential random numbers with parameter A.  The
%                   default is A = 1.
%   'erlang'        Erlang (gamma) random numbers with parameters A
%                   and B.  B must be a positive integer. The
%                   defaults are A = 2 and B = 5.  Erlang random
%                   numbers are approximated as the sum of B
%                   exponential random numbers.

%   Copyright 2002-2006 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2006/07/15 20:44:52 $

% Set default values.
if nargin == 1
   a = 0; b = 1;
   M = 1; N = 1;
elseif nargin == 3
   a = 0; b = 1;
end

% Begin processing. Use lower(type) to protect against input being
% capitalized. 
switch lower(type)
case 'uniform'
   R = a + (b - a)*rand(M, N);
case 'gaussian'
   R = a + b*randn(M, N);
case 'salt & pepper'
   if nargin <= 3
      a = 0.05; b = 0.05;
   end
   % Check to make sure that Pa + Pb is not > 1.
   if (a + b) > 1
      error('The sum Pa + Pb must not exceed 1.')
   end
   R(1:M, 1:N) = 0.5;
   % Generate an M-by-N array of uniformly-distributed random numbers
   % in the range (0, 1). Then, Pa*(M*N) of them will have values <=
   % a. The coordinates of these points we call 0 (pepper
   % noise). Similarly, Pb*(M*N) points will have values in the range
   % > a & <= (a + b).  These we call 1 (salt noise). 
   X = rand(M, N);
   c = find(X <= a);
   R(c) = 0;
   u = a + b;
   c = find(X > a & X <= u);
   R(c) = 1;
case 'lognormal'
   if nargin <= 3
      a = 1; b = 0.25;
   end
   R = exp(b*randn(M, N) + a);
case 'rayleigh'
   R = a + (-b*log(1 - rand(M, N))).^0.5;
case 'exponential'
   if nargin <= 3
      a = 1;
   end
   if a <= 0
      error('Parameter a must be positive for exponential type.')
   end
   k = -1/a;
   R = k*log(1 - rand(M, N));
case 'erlang'
   if nargin <= 3
      a = 2; b = 5;
   end
   if (b ~= round(b) | b <= 0)
      error('Param b must be a positive integer for Erlang.')
   end
   k = -1/a;
   R = zeros(M, N);
   for j = 1:b         
      R = R + k*log(1 - rand(M, N));
   end
otherwise
   error('Unknown distribution type.')
end

复原图像

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
51
52
53
54
% 生成一个8x8的棋盘图像并显示
f = checkerboard(8);
imshow(f, []);

% 创建一个长度为7、角度为45度的运动模糊点扩散函数 (PSF)
PSF = fspecial('motion', 7, 45);

% 使用圆周边界条件对棋盘图像应用运动模糊
gb = imfilter(f, PSF, 'circular');
figure;
imshow(gb, []);

% 生成均值为0、方差为0.001的高斯噪声，并且尺寸与图像相同
noise = imnoise(zeros(size(f)), 'gaussian', 0, 0.001);
figure;
imshow(noise, []);

% 将生成的噪声添加到模糊图像中
g = gb + noise;
figure;
imshow(g, []);

% 在没有噪声信号比的情况下执行维纳去卷积
fr1 = deconvwnr(g, PSF);

% 计算噪声的功率谱密度
Sn = abs(fft2(noise)) .^ 2;
nA = sum(Sn(:)) / prod(size(noise));

% 计算原始图像的功率谱密度
Sf = abs(fft2(f)) .^ 2;
fA = sum(Sf(:)) / prod(size(f));

% 计算噪声与信号功率比
R = nA / fA;

% 使用估计的噪声与信号比执行维纳去卷积
fr2 = deconvwnr(g, PSF, R);

% 计算噪声和原始图像的自相关函数
NCORR = fftshift(real(ifft2(Sn)));
ICORR = fftshift(real(ifft2(Sf)));

% 使用计算出的自相关函数执行维纳去卷积
fr3 = deconvwnr(g, PSF, NCORR, ICORR);

% 显示不同去卷积方法的结果
figure;
imshow(fr1, []);
figure;
imshow(fr2, []);
figure;
imshow(fr3, []);

空间噪声滤波器

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
clc
clear
f = imread('Fig0326(a)(embedded_square_noisy_512).tif');
figure;
imshow(f)
title('原始图像')

[M,N] = size(f);
R = imnoise2('salt & pepper',M,N,0.1,0);
c = find(R == 0);
gp = f;
gp(c) = 0;
figure;
imshow(gp)
title('被概率为0.1的胡椒噪声污染的图像')

R = imnoise2('salt & pepper',M,N,0,0.1);
c = find(R == 1);
gs = f;
gs(c) = 255;
figure;
imshow(gs)
title('被概率为0.1的盐粒噪声污染的图像')

fp = spfilt(gp,'chmean',3,3,1.5);
figure;
imshow(fp)
title('用阶为Q=1.5的3*3反调和滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果')

fs = spfilt(gs,'chmean',3,3,-1.5);
figure;
imshow(fs)
title('用阶为Q=-1.5的3*3反调和滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果')

fpmax = spfilt(gp,'max',3,3);
figure;
imshow(fpmax)
title('用3*3最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果')

fsmin = spfilt(gs,'min',3,3);
figure;
imshow(fsmin)
title('用3*3最小滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果')

自适应中值滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
clc
clear
f = imread('Fig0326(a)(embedded_square_noisy_512).tif');
figure;
imshow(f)
title('原始图像')

g = imnoise(f,'salt & pepper',0.25);% 噪声点有黑有白
figure;
imshow(g)
title('被概率为0.25椒盐噪声污染的图像')

f1 = medfilt2(g,[7 7],'symmetric');
figure;
imshow(f1)
title('用7*7中值滤波器对[被概率为0.25椒盐噪声污染的图像]滤波的结果')

f2 = adpmedian(g,7);
figure;
imshow(f2)
title('用Smax=7的自适应中值滤波器对[被概率为0.25椒盐噪声污染的图像]滤波的结果')

RGB | HSI

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
51
original_image=imread('1.tif');
% 显示原始RGB图像
figure;
subplot(2, 4, 1);
imshow(original_image);
title('Original RGB Image');

% 将RGB图像转换为HSI图像
hsi_image = rgb2hsv(original_image);

% 显示原始HSI分量图像
subplot(2, 4, 2);
imshow(hsi_image(:,:,1));
title('Original Hue');
subplot(2, 4, 3);
imshow(hsi_image(:,:,2));
title('Original Saturation');
subplot(2, 4, 4);
imshow(hsi_image(:,:,3));
title('Original Intensity');

% 修改HSI分量
modified_hsi_image = hsi_image;

% 色调增加30度
modified_hsi_image(:,:,1) = mod(hsi_image(:,:,1) + 30/360, 1); 

% 修改饱和度，将饱和度减半
modified_hsi_image(:,:,2) = hsi_image(:,:,2) * 0.5;

% 修改强度，将强度减半
modified_hsi_image(:,:,3) = hsi_image(:,:,3) * 0.5;

% 显示修改后的HSI分量图像
subplot(2, 4, 5);
imshow(modified_hsi_image(:,:,1));
title('Modified Hue');
subplot(2, 4, 6);
imshow(modified_hsi_image(:,:,2));
title('Modified Saturation');
subplot(2, 4, 7);
imshow(modified_hsi_image(:,:,3));
title('Modified Intensity');

% 将修改后的HSI图像转换回RGB图像
modified_rgb_image = hsv2rgb(modified_hsi_image);

% 显示修改后的RGB图像
subplot(2, 4, 8);
imshow(modified_rgb_image);
title('Modified RGB Image');

密度分层法

1
2
3
4
5
6
7
8
9
clc;
I=imread('1.tif');
I=rgb2gray(I);
GS8=grayslice(I,8);
GS64=grayslice(I,64);
subplot(1,3,1),imshow(I),title('原始灰度图像');
subplot(1,3,2),imshow (GS8, hot (8)) ,title('分成8层伪彩色');
subplot(1,3,3),imshow(GS64, hot(64)) ,title('分成64层伪彩色');

灰度变换法

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
clc;
I=imread('1.tif');
I=rgb2gray(I);
figure(1),imshow(I);
I=double(I);
[M,N]=size(I);
L=256;
for i=1:M
    for j=1:N
        if I(i,j)<=L/4
            R(i,j)=0;
            G(i,j)=4*I(i,j);
            B(i,j)=L;
        else if I(i,j)<=L/2
            R(i,j)=0;
            G(i,j)=L;
            B(i,j)=-4*I(i,j)+2*L;
            else if I(i,j)<=L/2
            R(i,j)=4*I(i,j)-2*L;
            G(i,j)=L;
            B(i,j)=0;
        else
                 R(i,j)=4*I(i,j)-2*L;
                 G(i,j)=L;
                 B(i,j)=0;
             end
            end
        end
    end
end
for i=1:M
    for j=1:N
        OUT(i,j,1)=R(i,j);
         OUT(i,j,2)=G(i,j);
          OUT(i,j,3)=B(i,j);
    end
end
OUT=OUT/256;
figure(2),imshow(OUT)

频域法

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
clc;
I=imread('1.tif');
I=rgb2gray(I);
figure(1),imshow(I);
[M,N]=size(I);
F=fft2(I);
fftshift(F);
REDcut=100;
GREENcut=200;
BLUEcenter=150;
BLUEwidth=100;
BLUEu0=10;
BLUEv0=10;
for u=1:M
    for v=1:N
        D(u,v)=sqrt(u^2+v^2);
        REDH(u,v)=1/(1+(sqrt(2)-1)*(D(u,v)/REDcut)^2);%红色滤波器为低通
        GREENH(u,v)=1/(1+(sqrt(2)-1)*(GREENcut/D(u,v))^2);%绿色滤波器为高通
        BLUED(u,v)=sqrt((u-BLUEu0)^2+(v-BLUEv0)^2);
        BLUEH(u,v)=1-1/(1+BLUED(u,v)*BLUEwidth/((BLUED(u,v))^2-(BLUEcenter)^2)^2);%蓝色滤波器为带通
    end
end
RED=REDH.*F;
REDcolor=ifft2(RED);
GREEN=GREENH.*F;
GREENcolor=ifft2(GREEN);
BLUE=BLUEH.*F;
BLUEcolor=ifft2(BLUE);
REDcolor=real(REDcolor)/256;
GREENcolor=real(GREENcolor)/256;
BLUEcolor=real(BLUEcolor)/256;
for i=1:M
    for j=1:N
        OUT(i,j,1)=REDcolor(i,j);
        OUT(i,j,2)=GREENcolor(i,j);
        OUT(i,j,3)=BLUEcolor(i,j);
    end
end
OUT=abs(OUT);
figure,imshow(OUT);

彩色图像平滑滤波

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
%彩色图像空间滤波
clc;
clear all;
close all;
disp('彩色图像空间滤波开始.......');
%%
%提取3个分量图像
f=imread('1.tif');   %加载彩色图像
%显示原图像
figure;imshow(f);title('彩色原图像');
fr=f(:,:,1);    %提取R通道分量图像
fg=f(:,:,2);    %提取G通道分量图像      
fb=f(:,:,3);    %提取B通道分量图像
%显示三通道图像
figure;
subplot(2,2,1);imshow(fr);title('R');
subplot(2,2,2);imshow(fg);title('G');
subplot(2,2,3);imshow(fb);title('B');
%%
%分别过滤每个分量图像
w = fspecial('average', 3);
fr_filter=imfilter(fr,w,'replicate');   %平滑红色分量图像
fg_filter=imfilter(fg,w,'replicate');   %平滑绿色分量图像
fb_filter=imfilter(fb,w,'replicate');   %平滑蓝色分量图像
%显示滤波后的三通道图像
figure;
subplot(2,2,1);imshow(fr_filter);title('R滤波后');
subplot(2,2,2);imshow(fg_filter);title('G滤波后');
subplot(2,2,3);imshow(fb_filter);title('B滤波后');

%%
%重建滤波后的RGB图像
ff=cat(3,fr_filter,fg_filter,fb_filter); %构造多维数组，即合并3分量图像为一副彩色图像
%显示重建后的图像
figure;imshow(ff);title('重建后');
%%
%使用与单色图像相同的语法来执行RGB 图像的线性滤波，可以把前三步合并为一步：
gg = imfilter(f, w, 'replicate');
figure;imshow(gg);title('步骤合并的滤波结果');

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
clc
clear
fc  = imread('1.tif');
figure;imshow(fc)
title('原始真彩色（256*256*256色）图像')

fr = fc(:,:,1);
fg = fc(:,:,2);
fb = fc(:,:,3);
% 
figure;imshow(fr)
title('红色分量图像')
figure;imshow(fg)
title('绿色分量图像')
figure;imshow(fb)
title('蓝色分量图像')

h = rgb2hsi(fc);
H = h(:,:,1);
S = h(:,:,2);
I = h(:,:,3);
figure;imshow(H)
title('色调分量图像')
figure;imshow(S)
title('饱和度分量图像')
figure;imshow(I)
title('亮度分量图像')

w = fspecial('average',15);
I_filtered = imfilter(I,w,'replicate');
h = cat(3,H,S,I_filtered);
f = hsi2rgb(h);
f = min(f,1);
figure;imshow(f)
title('仅平滑HSI图像的亮度分量所得到的RGB图像')

fc_filtered = imfilter(fc,w,'replicate');
figure;imshow(fc_filtered)
title('分别平滑R、G、B图像分量平面得到的RGB图像')

h_filtered = imfilter(h,w,'replicate');
f = hsi2rgb(h_filtered);
f = min(f,1);
figure;imshow(f)
title('分别平滑H、S、I图像分量平面得到的RGB图像')
% 
h_filtered = imfilter(h,w,'replicate');
figure;imshow(h_filtered)
title('分别平滑H、S、I图像分量平面得到的HSI图像')

彩色图像锐化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clc
clear
fc  = imread('1.tif');
figure;imshow(fc)
title('原始真彩色（256*256*256色）图像')

w = fspecial('average',15);
fc_filtered = imfilter(fc,w,'replicate');
figure;imshow(fc_filtered)
title('分别平滑R、G、B图像分量平面得到的RGB模糊图像')

lapmask = [1 1 1; 1 -8 1; 1 1 1];

fen = imsubtract(fc_filtered,imfilter(fc_filtered,lapmask,'replicate'));
figure;imshow(fen)
title('用拉普拉斯算子增强模糊图像（效果好象不是很明显！）')

LPA = imfilter(fc,lapmask,'replicate');
figure;imshow(LPA)
title('对原始真彩色图像用拉普拉斯算子提取出的图像')

fen = imsubtract(fc,imfilter(fc,lapmask,'replicate'));
figure;imshow(fen)
title('用拉普拉斯算子增强原始真彩色图像（采用提高边缘亮度手段）')

imdilate 膨胀

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
clc
clear

A = imread('Fig0929(a)(text_image).tif');
info = imfinfo('Fig0929(a)(text_image).tif')
B = [0 1 0
     1 1 1
     0 1 0];
 A2 = imdilate(A, B);
 A3 = imdilate(A2, B); % 二次膨胀
 A4 = imdilate(A3, B); % 三次膨胀

figure;
imshow(A)
title('原始图像')

figure;
imshow(A2)
title('使用结构元素[B]一次膨胀后的图像')

figure;
imshow(A3)
title('使用结构元素[B]二次膨胀后的图像')

figure;
imshow(A4)
title('使用结构元素[B]三次膨胀后的图像')

figure;
imshow(A2-A) % 显示增加的部分
title('使用结构元素[B]一次膨胀后和原图像相比较增加的部分')

腐蚀

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clc
clear

A = imread('FigP0918(left).tif');
se = strel('disk', 10);
figure;
imshow(A)
title('原始图像')

A2 = imerode(A, se);
figure;
imshow(A2)
title('使用结构元素[disk（10）]腐蚀后的图像')

se = strel('disk', 5);
A3 = imerode(A, se);
figure;
imshow(A3)
title('使用结构元素[disk（5）]腐蚀后的图像')

A4 = imerode(A, strel('disk', 20));
figure;
imshow(A4)
title('使用结构元素[disk（20）]腐蚀后的图像')

imopen imclose

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
clc
clear

f = imread('Fig0911(a)(noisy_fingerprint).tif');
% se = strel('square', 5);  % 结构元素 方型
se = strel('disk', 1);  % 结构元素 圆盘形
figure;
imshow(f)  % 原始图  图1
title('原始图像')

fo = imopen(f, se);  % 开 图2
figure;
imshow(fo) 
title('使用结构元素[disk]开操作后的图像')

fc = imclose(f, se); % 闭 图3
figure;
imshow(fc)
title('使用结构元素[disk]闭操作后的图像')

foc = imclose(fo, se);  % 先开再闭 图4
figure;
imshow(foc)
title('使用结构元素[disk]先开操作再闭操作后的图像')

fco = imopen(fc, se);   % 先闭再开 图5
figure;
imshow(fco)  
title('使用结构元素[disk]先闭操作再开操作后的图像')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% 先膨胀再腐蚀 图6
fse = imdilate(f, se);
figure,set(gcf,'outerposition',get(0,'screensize'))
subplot(211),imshow(fse)
title('使用结构元素[disk（5）]先膨胀后的图像')
fes = imerode(fse, se);
subplot(212),imshow(fes)  
title('使用结构元素[disk（5）]先膨胀再腐蚀后的图像')

% 先腐蚀再膨胀 图7
fse = imerode(f, se);
figure, set(gcf,'outerposition',get(0,'screensize'))
subplot(211),imshow(fse)
title('使用结构元素[disk（5）]先腐蚀后的图像')
fes = imdilate(fse, se);
subplot(212),imshow(fes)  
title('使用结构元素[disk（5）]先腐蚀再膨胀后的图像')

指纹

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
%% 例9.4.2 imopen imclose 指纹 
clc
clear

f = imread('Fig0911(a)(noisy_fingerprint).tif');
se = strel('square', 3);  % 结构元素
% se = strel('disk', 2);  % 结构元素 圆盘形

figure;
imshow(f)  % 原始图
title('原始图像')

A = imerode(f, se); % 腐蚀
figure;
imshow(A)
title('使用结构元素[square（3）]腐蚀后的图像')

fo = imopen(f, se);
figure;
imshow(fo)    % 开
title('使用结构元素[square（3）]开操作后的图像')

fc = imclose(f, se);   % 闭
figure;
imshow(fc)
title('使用结构元素[square（3）]闭操作后的图像')

foc = imclose(fo, se);  % 先开再闭
figure;
imshow(foc)
title('使用结构元素[square（3）]先开操作再闭操作后的图像')

fco = imopen(fc, se);   % 先闭再开
figure;
imshow(fco)
title('使用结构元素[square（3）]先闭操作再开操作后的图像')

endpoints

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
clc
clear

f = imread('Fig0914(a)(bone-skel).tif');
figure;
imshow(f)
title('原始形态骨骼的图像')

g = endpoints(f);
figure;
imshow(g)
title('使用函数[endpoints]后得到的端点图像')

f = imread('Fig0916(a)(bone).tif');
figure;
imshow(f)
title('原始骨头图像')

g = endpoints(f);
figure;
imshow(g)
title('使用函数[endpoints]后得到的端点图像（什么也没有）')
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
function g = endpoints(f)
%ENDPOINTS Computes end points of a binary image.
%   G = ENDPOINTS(F) computes the end points of the binary image F
%   and returns them in the binary image G. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.3 $  $Date: 2003/04/22 01:18:18 $

persistent lut

if isempty(lut)
   lut = makelut(@endpoint_fcn, 3);
end

g = applylut(f,lut);

%-------------------------------------------------------------------%
function is_end_point = endpoint_fcn(nhood)
%   Determines if a pixel is an end point.
%   IS_END_POINT = ENDPOINT_FCN(NHOOD) accepts a 3-by-3 binary
%   neighborhood, NHOOD, and returns a 1 if the center element is an
%   end point; otherwise it returns a 0. 

is_end_point = nhood(2,2) & (sum(nhood(:)) == 2);

bwmorph

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
%% bwmorph   'bridge' 'clean' 'hbreak'
clc
clear

f = imread('Fig0911(a)(noisy-fingerprint).tif');

figure;
imshow(f)  % 原始图
title('原始图像')

BW3 = bwmorph(f,'bridge',Inf);
figure;
imshow(BW3)

BW3 = bwmorph(BW3,'hbreak',Inf); % 极其细微的变化
figure;
imshow(BW3)

BW3 = bwmorph(f,'clean',Inf);
figure;
imshow(BW3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clc
clear

f = imread('Fig0911(a)(noisy-fingerprint).tif');
% f = imread('..\Pictures\Beautiful\hehe1.tif');
figure;
imshow(f)
title('原始指纹图像')

g1 = bwmorph(f, 'thin', 1);
figure;
imshow(g1)
title('使用函数[bwmorph]细化一次后的图像')

g2 = bwmorph(f, 'thin', 2);
figure;
imshow(g2)
title('使用函数[bwmorph]细化两次后的图像')

% 细化到稳定状态
ginf = bwmorph(f, 'thin', Inf); 
figure;
imshow(ginf)
title('使用函数[bwmorph]细化到稳定状态后的图像')
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
clc
clear

f = imread('Fig0911(a)(noisy-fingerprint).tif');
% f = imread('..\Pictures\Beautiful\hehe1.tif');
figure;
imshow(f)
title('原始指纹图像')

g1 = bwmorph(f, 'skel', 1);
figure;
imshow(g1)
title('使用函数[bwmorph]骨骼化一次后的图像')

g2 = bwmorph(f, 'skel', 2);
figure;
imshow(g2)
title('使用函数[bwmorph]骨骼化两次后的图像')

% 骨骼化到稳定状态
fs = bwmorph(f, 'skel', Inf); 
figure;
imshow(fs)
title('使用函数[bwmorph]骨骼化到稳定状态后的图像')

for k = 1:5
    fs = fs & ~endpoints(fs);
end
figure;
imshow(fs)
title('使用函数[endpoints]五次修剪骨骼端点后的图像')