《数字图像处理》实验3-频率域处理方法

% 读取彩色图像
img_rgb = imread('阿狸.jpg'); 

% 将彩色图像转换为灰度图
img_gray = rgb2gray(img_rgb); 

% 对灰度图进行二维傅里叶变换（需转换为double类型）
F = fft2(double(img_gray)); 

% 频谱中心化（将低频成分移到图像中心）
F_shift = fftshift(F); 

% 计算频谱的幅值，并取对数以便于显示（避免数值范围过大导致显示失真）
F_mag = log(1 + abs(F_shift)); 

% 创建figure窗口，展示结果
figure('Name','傅里叶变换实验','NumberTitle','off');
subplot(1,3,1); imshow(img_rgb); title('原图');
subplot(1,3,2); imshow(img_gray); title('灰度图');
subplot(1,3,3); imshow(F_mag, []); title('变换频谱'); % []自动映射灰度范围

% 读取彩色图像
img_rgb = imread('天火三玄变.jpg'); 

% 转换为灰度图
img_gray = rgb2gray(img_rgb);

% 转为double类型以进行傅里叶变换
img_gray_double = double(img_gray);

% 二维傅里叶变换
F = fft2(img_gray_double);

% 频谱中心化（用于可视化）
F_shift = fftshift(F);
F_mag = log(1 + abs(F_shift));

% 傅里叶逆变换（还原回空间域）
img_recon = ifft2(F);

% 逆变换结果可能包含极小虚部，取实部并转换回uint8类型
img_recon = uint8(real(img_recon));

% 比较原图与逆变换图像矩阵是否相同
isSame = isequal(img_gray, img_recon);
if isSame
    fprintf('原图与逆变换图像矩阵是否完全相同: 是\n');
else
    fprintf('原图与逆变换图像矩阵是否完全相同: 否（微小差异由浮点运算精度导致）\n');
end

% 可视化结果
figure('Name', 'FFT 与逆变换实验', 'NumberTitle', 'off');
subplot(1,4,1); imshow(img_rgb); title('原图');
subplot(1,4,2); imshow(img_gray); title('灰度图');
subplot(1,4,3); imshow(F_mag, []); title('二维傅里叶变换');
subplot(1,4,4); imshow(img_recon); title('傅里叶逆变换');

% 读取两张彩色图像
img1_rgb = imread('ALi.jpg'); 
img2_rgb = imread('天火三玄变.jpg'); 

% 转换为灰度图
img1_gray = rgb2gray(img1_rgb);
img2_gray = rgb2gray(img2_rgb);

% 转换为double类型以进行傅里叶变换
img1_double = double(img1_gray);
img2_double = double(img2_gray);

% 二维傅里叶变换
F1 = fft2(img1_double);
F2 = fft2(img2_double);

% 频谱中心化
F1_shift = fftshift(F1);
F2_shift = fftshift(F2);

% 计算幅值谱，并取对数以便显示
F1_mag = log(1 + abs(F1_shift));
F2_mag = log(1 + abs(F2_shift));

% 可视化对比结果
figure('Name', '傅里叶频谱对比', 'NumberTitle', 'off');
subplot(2,2,1); imshow(img1_gray); title('原图灰度图1');
subplot(2,2,2); imshow(img2_gray); title('原图灰度图2');
subplot(2,2,3); imshow(F1_mag, []); title('频谱图1');
subplot(2,2,4); imshow(F2_mag, []); title('频谱图2');

% 读取彩色图像并转换为灰度图
img_rgb = imread('阿狸2.jpg'); 
img_gray = rgb2gray(img_rgb);
img_double = double(img_gray);

%% 方法1：直接使用fft2进行二维傅里叶变换
F_fft2 = fft2(img_double);
F_fft2_shift = fftshift(F_fft2);
F_fft2_mag = log(1 + abs(F_fft2_shift));

%% 方法2：使用fft分两步实现（先对行，再对列）
F_row = fft(img_double, [], 2);  % 对每一行做一维傅里叶变换
F_2step = fft(F_row, [], 1);     % 对结果的每一列做一维傅里叶变换
F_2step_shift = fftshift(F_2step);
F_2step_mag = log(1 + abs(F_2step_shift));

%% 可视化两种方法的结果
figure('Name', '二维DFT分解验证', 'NumberTitle', 'off');
subplot(1,3,1); imshow(img_gray); title('原图');
subplot(1,3,2); imshow(F_fft2_mag, []); title('fft2实现变换');
subplot(1,3,3); imshow(F_2step_mag, []); title('fft实现变换');

%% 验证两种方法结果的一致性（允许微小浮点误差）
diff_matrix = abs(F_fft2 - F_2step);
max_error = max(diff_matrix(:));
fprintf('两种方法结果的最大误差: %.2e\n', max_error);
if max_error < 1e-8
    fprintf('验证结论：两种方法在理论上一致，微小误差由浮点运算精度导致。\n');
else
    fprintf('验证结论：结果差异超出预期，需检查代码。\n');
end

% 配置参数：学号末位数M=7，截止频率D0=M*10=70
student_id_last_digit = 7;
D0 = student_id_last_digit * 10;  % 理想低通滤波器截止频率

% 读取彩色图像并转换为灰度图
img_raw = imread('AALi.jpg'); 
img_gray = rgb2gray(img_raw);

% 为灰度图添加高斯白噪声
img_noise = imnoise(img_gray, 'gaussian');

% 对加噪图像做傅里叶变换并中心化（用于可视化和滤波）
fft_img_noise = fft2(double(img_noise));
fftshift_img_noise = fftshift(fft_img_noise);

% 获取图像尺寸，构造理想低通滤波器
[rows, cols] = size(img_noise);
H = zeros(rows, cols);  % 初始化滤波器掩模
center_row = rows / 2;  % 频域中心行坐标
center_col = cols / 2;  % 频域中心列坐标

% 遍历每个频率点，计算到中心的距离并生成滤波器
for i = 1:rows
    for j = 1:cols
        % 计算频率点(i,j)到中心的欧氏距离
        distance = sqrt((i - center_row)^2 + (j - center_col)^2);
        % 理想低通滤波器：距离≤D0则保留（H=1），否则抑制（H=0）
        if distance <= D0
            H(i, j) = 1;
        else
            H(i, j) = 0;
        end
    end
end

% 频域滤波：滤波器掩模与频域图像逐元素相乘
filtered_fftshift = fftshift_img_noise .* H;

% 逆傅里叶变换，还原回空间域
filtered_img = ifft2(ifftshift(filtered_fftshift));
filtered_img = real(filtered_img);  % 去除微小虚部

% 计算滤波后图像的频谱（用于可视化）
filtered_fft = fftshift(fft2(filtered_img));

% 可视化所有结果
figure('Name', '高斯加噪与理想低通滤波', 'NumberTitle', 'off');
subplot(2, 2, 1); imshow(img_noise); title('加噪后图像');
subplot(2, 2, 2); imshow(log(abs(fftshift_img_noise) + 1), []); title('噪声fft');
subplot(2, 2, 3); imshow(filtered_img, []); title('滤波后噪声图');
subplot(2, 2, 4); imshow(log(abs(filtered_fft) + 1), []); title('滤波后噪声图fft');