数字图像处理-常考算法代码-详细注释

求一幅图像的灰度直方图

%% my 直方图
clc, clear, close all;

I = imread('circuit.tif');

P = zeros(1, 256);

for i = 1:size(I, 1)
    for j = 1:size(I, 2)
        k = I(i, j); 
        P(k+1) = P(k+1) + 1; % k 范围为 0~255，由于 matlab 数组下标从 1 开始，所以 k+1
    end
end

P = P / numel(I); % numel() 返回数组中元素的个数

plot(0:255, P);
xlabel('k');
ylabel('P(k)');
set(gca, 'xlim', [0 255]); % 设置x轴坐标域
set(gca, 'ylim', [0 max(P)*1.05]); % 设置y轴坐标域

对一幅图像的直方图均衡化

%% my 直方图均衡化
clc, clear, close all;

I = imread('circuit.tif');

subplot(2, 2, 1);
imshow(I);
title('原图');

subplot(2, 2, 2);
imhist(I); % imhist() 显示图像直方图
title('直方图');

I1 = histeq(I); % histeq() 直方图均衡化

subplot(2, 2, 3);
imshow(I1);
title('均衡化后原图');

subplot(2, 2, 4);
imhist(I1);
title('均衡化后直方图');

使用各种模板的中值滤波对图像进行处理

中值滤波：空域滤波的非线性滤波，平滑去噪（椒盐噪声）

%% my 中值滤波（空域滤波的非线性滤波，平滑去噪）
clc, clear, close all;
I = imread('circuit.tif');
J = imnoise(I, 'salt & pepper', 0.02); % imnoise() 添加噪声

subplot(2, 3, 1);
imshow(I);
title('原图像');

subplot(2, 3, 2);
imshow(J);
title('椒盐噪声');

k1 = medfilt2(J); % medfilt2(J) 中值滤波，默认 3 * 3 模板
subplot(2, 3, 3);
imshow(k1);

k2 = medfilt2(J, [5 5]);
subplot(2, 3, 4);
imshow(k2);

k3 = medfilt2(J, [7 7]);
subplot(2, 3, 5);
imshow(k3);

k4 = medfilt2(J, [9 9]);
subplot(2, 3, 6);
imshow(k4);

对图像进行傅里叶变换

（注：结果还是灰度图像）

%% my 傅里叶变换
clc, clear, close all;

d = zeros(32, 32);
d(13:20, 13:20) = 1;
subplot(1, 2, 1);
imshow(d);

D = fft2(d); % fft2() 二维快速傅里叶变换
subplot(1, 2, 2);
imshow(abs(D), [-1 5]);

使用巴特沃斯低通滤波器对图像进行处理

%% my 巴特沃斯低通滤波器 平滑，去噪
clc, clear, close all;

I = imread('saturn.tif');
J = imnoise(I, 'salt & pepper', 0.02);
subplot(1, 2, 1);
imshow(J);
title('椒盐噪声');
J = double(J);

f = fft2(J); % 傅里叶变换
g = fftshift(f); % 数据矩阵平衡
[M,N] = size(f); % 获取频谱矩阵大小
n = 3; % 巴特沃斯阶数
d0 = 20; % 截至频率 20
n1 = floor(M/2); % 求频谱中心
n2 = floor(N/2); % 求频谱中心

for i = 1:M
    for j = 1:N
        d = sqrt((i-n1)^2 + (j-n2)^2); % 计算距离频谱中心距离 d
        h = 1/(1+(d/d0)^(2*n));
        g(i, j) = h * g(i, j);
    end
end
g = ifftshift(g); % 逆数据矩阵平衡
g = uint8(real(ifft2(g)));  % 逆傅里叶变换,real返回数组 Z 中每个元素的实部
subplot(1, 2, 2);
imshow(rgb2gray(g), []); % 显示图像，并将显示范围缩放到图像中的像素值。图像以完整范围的灰度值显示。

使用迭代阈值分割图像

%% my 迭代阈值分割
clc;clear;close;

f = imread('barbara.bmp'); 
subplot(1,2,1); 
imshow(f); 
title('原始图像');

f = double(f);
T = (min(f(:)) + max(f(:)))/2; % 1.求初始 T0
done = false;

while ~done
    r1 = find(f <= T); % 2.划分 r1, r2
    r2 = find(f > T);
    Tnew = (mean(f(r1)) + mean(f(r2))) / 2; % 3. 更新阈值 Tnew
    done = abs(Tnew - T) < 1; % 4. 判断终止条件
    T = Tnew;
    i = i + 1;
end
f(r1) = 0;
f(r2) = 1;
subplot(1, 2, 2);
imshow(f);
title('迭代阈值处理后的图像');

使用区域生长法分割图像

%% my 区域生长法分割
clc, clear, close all;
f1 = imread('pepper.bmp');
f = double(rgb2gray(f1));
subplot(1,2,1);
imshow(f1);
seedx = [50,76,86]; %选择种子像素点
seedy = [150,81,110]; %选择种子像素点
hold on; 
plot(seedx, seedy, 'gs', 'linewidth', 1); 
title('原图像及种子点位置');
f = double(f);

markerim = f == f(seedy(1),seedx(1)); % 种子像素点

for i=2:length(seedx)
    markerim = markerim | (f==f(seedy(i),seedx(i))); % 种子像素点
end

thresh = [12, 6, 12];
maskim = zeros(size(f));
for i = 1:length(seedx)
    g = abs(f-f(seedy(i),seedx(i)))<=thresh(i); % 扩展区域，区域生长
    maskim = maskim | g;
end

bw1 = imreconstruct(markerim, maskim);

[g,nr] = bwlabel(bw1,8);
g = mat2gray(g);
subplot(1,2,2);
imshow(bw1); 
title('三个种子点区域生长结果');

使用背景差分法分割图像

该方法局限性大，不是很实用

%% my 背景差分法分割
clc, clear, close all;

a = imread('cat.bmp'); % imread() 读取图片，一般读取到的都是 uint8; 
subplot(2, 2, 1);
imshow(a); % imshow() 显示图片，一般为 uint8 类型，可以使用 im2uint8() 将图像转换为 8位无符号整数
title('原图');

b = imread('background.bmp');
subplot(2, 2, 2);
imshow(b);
title('背景');

aa = im2double(a); % 对图像进行加减处理时，要使用 im2double() 将图像从 8位无符号整数 转换为 双精度
bb = im2double(b);

c = aa - bb;
d = im2uint8(c);
subplot(2, 2, 3);
imshow(d);
title('差分后图像');

T = 50;
T = T/255;
i = find(abs(c)>=T); % find() 寻找满足关系式的数组下标
c(i) = 1;
i = find(abs(c)<T);
c(i) = 0;
e = im2uint8(c);
subplot(2, 2, 4);
imshow(e);
title('二值化后的图像');

祝大家考个好成绩👍👍👍