matlab去除横条纹-巴斯沃特低通滤波+中值滤波

clc; clear; close all; % 清空环境，开始一个干净的会话
% 步骤1：读取并显示原始图像
img_original = imread('40-NQ-15.jpg');
% 确保为灰度图像，如果原图是彩色的，需要转换
if size(img_original, 3) == 3
img_gray = rgb2gray(img_original);
else
img_gray = img_original;
end
figure;
imshow(img_gray);
title('原始带条纹噪声图像', 'FontSize', 11);
axis on;
% 步骤2：进行二维傅里叶变换并中心化（得到频域图像）
F = fft2(double(img_gray)); % fft2计算二维离散傅里叶变换
F_shifted = fftshift(F); % fftshift将零频率分量移到频谱中心
% 计算对数变换后的幅度谱，便于显示（因为动态范围太大）
magnitude_spectrum = log(abs(F_shifted) + 1);
subplot(1,3,2);
imshow(magnitude_spectrum, []);
title('中心化幅度谱（对数尺度）', 'FontSize', 11);
colormap jet; colorbar;
axis on;
% 步骤3：分析频谱，定位噪声频率
[M, N] = size(img_gray); % M行，N列
center_row = round(M/2);
center_col = round(N/2);
% 提取通过频谱中心水平线的幅度值
horizontal_profile = abs(F_shifted(center_row, :));
subplot(1,3,3);
plot(1:N, horizontal_profile, 'LineWidth', 1.5);
xlabel('水平频率坐标 (u)');
ylabel('幅度值');
title('频谱中心水平线剖面图');
grid on;
hold on;
% 标记出中心点（零频率）
plot(center_col, horizontal_profile(center_col), 'ro', 'MarkerSize', 8, 'LineWidth', 2);
% 步骤4：根据定位的噪声频率，设计巴特沃斯带阻滤波器
% 假设通过步骤3的剖面图，我们确定 u0 = 56 (举例)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u0 = 27; % 噪声频率点（距中心）
v0 = 9; % 假设噪声在竖直方向无变化，故在中心线上
D0 = 30; % 阻带宽度，需要根据噪声尖峰的宽度微调
n = 6; % 滤波器阶数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 创建网格坐标。注意：经过fftshift后，频谱中心在(M/2+1, N/2+1)
% 我们的u0是相对于中心的偏移，所以需要构建以频谱中心为原点的坐标
[U, V] = meshgrid(1:N, 1:M);
% 将坐标原点平移到频谱中心
U_centered = U - center_col;
V_centered = V - center_row;
% 计算到两个噪声点的距离
D1 = sqrt( (U_centered - u0).^2 + (V_centered - v0).^2 );
D2 = sqrt( (U_centered + u0).^2 + (V_centered - v0).^2 );
% 避免分母为零（虽然在实际中，D1或D2为零的点极少，但为计算稳定）
D1(D1 == 0) = eps;
D2(D2 == 0) = eps;
% 计算巴特沃斯带阻滤波器传递函数
H = 1 ./ (1 + ( (D0^2) ./ (D1 .* D2) ).^n );
% 可视化滤波器
figure('Position', [100, 100, 1000, 300]);
subplot(1,3,1);
imshow(H, []);
title('带阻滤波器（图像表示）', 'FontSize', 11);
colorbar;
subplot(1,3,2);
mesh(U_centered(1:10:end, 1:10:end), V_centered(1:10:end, 1:10:end), H(1:10:end, 1:10:end));
title('带阻滤波器三维透视图', 'FontSize', 11);
xlabel('u'); ylabel('v'); zlabel('H(u,v)');
axis tight;
% 绘制滤波器水平剖面，看其形状
subplot(1,3,3);
plot(U_centered(center_row, :), H(center_row, :), 'b-', 'LineWidth', 2);
xlabel('水平频率 u (中心为0)');
ylabel('滤波器响应 H');
title('滤波器中心水平剖面 (v=0)');
grid on;
hold on;
% 标记阻带中心
plot([-u0, u0], [0, 0], 'r^', 'MarkerSize', 10, 'LineWidth', 2);
legend('滤波器响应', '噪声频率位置');
% 步骤5：应用滤波器并反变换回空间域
G_shifted = F_shifted .* H; % 频域滤波：点乘
g_filtered = real(ifft2(ifftshift(G_shifted))); % 反中心化，反变换，取实部
% 步骤6：结果可视化与对比
figure('Position', [100, 100, 1400, 400]);
% 子图1：原始图像
subplot(2,4,1);
imshow(img_gray, []);
title('(a) 原始噪声图像');
% 子图2：原始频谱
subplot(2,4,2);
imshow(magnitude_spectrum, []);
title('(b) 原始幅度谱');
% 子图3：设计的滤波器
subplot(2,4,3);
imshow(H, []);
title('(c) 带阻滤波器 H(u,v)');
% 子图4：滤波后的频谱
filtered_magnitude = log(abs(G_shifted) + 1);
subplot(2,4,4);
imshow(filtered_magnitude, []);
title('(d) 滤波后幅度谱');
% 子图5：滤波后图像
subplot(2,4,5);
imshow(g_filtered, []);
title('(e) 滤波后图像');
% 子图6：局部放大对比 (例如，取图像左上角200x200区域)
rect_region = [1, 1, 199, 199]; % [x, y, width, height]
img_roi_original = imcrop(img_gray, rect_region);
img_roi_filtered = imcrop(g_filtered, rect_region);
subplot(2,4,6);
imshow(img_roi_original, []);
title('(f) 原始图像局部');
% 子图7：滤波后局部
subplot(2,4,7);
imshow(img_roi_filtered, []);
title('(g) 滤波后局部');
% 子图8：差异图（突出被去除的噪声）
difference = double(img_gray) - double(g_filtered);
% 将差异归一化到可显示范围
difference_normalized = mat2gray(difference);
subplot(2,4,8);
imshow(difference_normalized, []);
title('(h) 被滤除的成分（差异图）');
colormap(gca, 'parula'); colorbar;

%% step2:空间域处理-中值滤波(消除随机噪声+部分随机条纹)
medFiltSize = [10,1]; %垂直方向中值滤波，适用于水平条纹
medFiltered = medfilt2(g_filtered,medFiltSize);
imshow(medFiltered);

I22 = im2bw(medFiltered,0.5);
I22 = ~I22;
for i = 1:10
for j = 1:N
I22(i,j) = 0;
end
end
for i = M-10:M
for j = 1:N
I22(i,j) = 0;
end
end
figure;
imshow(I22);
figure;
imshow(img_original);
[L,num] = bwlabel(I22);
STATS = regionprops(L,'Area');
fprintf('%d');
for k = 1:num
fprintf('Region %d has %d pixels.\n', k, STATS(k).Area);
end
total_sum = 0;
for i = 1:length(STATS)
total_sum = total_sum + STATS(i).Area;
end
total_area_lv = total_sum *100/ M/N;