发帖
雷达卡
分数阶微积分在图像去噪中表现良好，因为它能够同时保持图像的边缘和纹理细节，同时有效地抑制噪声。以下是一个完整的分数阶微积分图像去噪MATLAB实现，涵盖了多种方法和对比分析。
1. 基于分数阶微分掩模的去噪方法
% 基于分数阶微分掩模的图像去噪
function fractional_order_denoising()
% 清除环境
clear; clc; close all;
% 读取并添加噪声
I = im2double(rgb2gray(imread('lena.png')));
I_noisy = imnoise(I, 'gaussian', 0, 0.02);
% 分数阶微分去噪参数
v = 0.8; % 分数阶次 (0 < v < 1)
k = 3; % 掩模尺寸
% 构建分数阶微分掩模
mask = construct_frac_diff_mask(v, k);
% 应用分数阶微分去噪
I_denoised = imfilter(I_noisy, mask, 'replicate');
% 展示结果
figure;
subplot(1,3,1); imshow(I); title('原始图像');
subplot(1,3,2); imshow(I_noisy); title('含噪图像 (PSNR='+string(psnr(I_noisy,I))+'dB)');
subplot(1,3,3); imshow(I_denoised); title(['分数阶去噪(v=',num2str(v),', PSNR=',...
num2str(psnr(I_denoised,I)),'dB)']);
% 计算SSIM
disp(['含噪图像SSIM: ', num2str(ssim(I_noisy, I))]);
disp(['去噪后SSIM: ', num2str(ssim(I_denoised, I))]);
end
% 构建分数阶微分掩模函数
function mask = construct_frac_diff_mask(v, k)
mask = zeros(2*k+1);
center = k+1;
% 根据Grünwald-Letnikov定义构建掩模
for i = 1:2*k+1
for j = 1:2*k+1
if i == center && j == center
mask(i,j) = gamma(2-v)/(gamma(1-v)^2);
else
r = sqrt((i-center)^2 + (j-center)^2);
mask(i,j) = (-1)^(i+j) * gamma(1-v) / ...
(gamma(i-center-v+1)*gamma(j-center-v+1));
end
end
end
% 归一化处理
mask = mask / sum(abs(mask(:)));
end
2. 分数阶全变分(TV)去噪方法
% 分数阶全变分去噪方法
function fractional_TV_denoising()
% 清除环境
clear; clc; close all;
% 读取并添加噪声
I = im2double(imread('lena.tif'));
I_noisy = imnoise(I, 'salt & pepper', 0.05);
% 算法参数
alpha = 1.2; % 分数阶次 (1 < α ≤ 2)
lambda = 0.1; % 正则化参数
iter = 100; % 迭代次数
dt = 0.05; % 时间步长
% 初始化
u = I_noisy;
% 迭代求解
for n = 1:iter
% 计算分数阶梯度
[Gx, Gy] = fractional_gradient(u, alpha);
% 计算扩散系数
c = 1 ./ (sqrt(Gx.^2 + Gy.^2) + eps);
% 计算分数阶散度
div = fractional_divergence(c.*Gx, c.*Gy, alpha);
% 更新图像
u = u + dt * (div - lambda*(u - I_noisy));
% 展示中间结果
if mod(n,20) == 0
figure(1);
subplot(1,2,1); imshow(I_noisy); title('含噪图像');
subplot(1,2,2); imshow(u);
title(['迭代 ',num2str(n),', α=',num2str(alpha),...
', PSNR=',num2str(psnr(u,I)),'dB']);
drawnow;
end
end
% 展示最终结果
figure;
subplot(1,3,1); imshow(I); title('初始图像');
subplot(1,3,2); imshow(I_noisy); title('含噪图像');
subplot(1,3,3); imshow(u); title(['分数阶TV去噪(α=',num2str(alpha),')']);
% 评估标准
disp(['含噪图像PSNR: ', num2str(psnr(I_noisy,I))]);
disp(['去噪后PSNR: ', num2str(psnr(u,I))]);
disp(['含噪图像SSIM: ', num2str(ssim(I_noisy,I))]);
disp(['去噪后SSIM: ', num2str(ssim(u,I))]);
end
% 分数阶导数计算
function [Gx, Gy] = fractional_gradient(u, alpha)
[m,n] = size(u);
Gx = zeros(m,n);
Gy = zeros(m,n);
% x方向分数阶差分 (基于Grünwald-Letnikov定义)
for i = 1:m
for j = 2:n
Gx(i,j) = u(i,j) - u(i,j-1);
for k = 1:j-2
Gx(i,j) = Gx(i,j) + (-1)^k * gamma(alpha+1) / ...
(gamma(k+1)*gamma(alpha-k+1)) * u(i,j-k-1);
end
end
end
% y方向分数阶差分
for j = 1:n
for i = 2:m
Gy(i,j) = u(i,j) - u(i-1,j);
for k = 1:i-2
Gy(i,j) = Gy(i,j) + (-1)^k * gamma(alpha+1) / ...
(gamma(k+1)*gamma(alpha-k+1)) * u(i-k-1,j);
end
end
end
end
% 分数阶散度计算
function div = fractional_divergence(u1, u2, alpha)
[m,n] = size(u1);
div = zeros(m,n);
% x方向
for i = 1:m
for j = 1:n-1
div(i,j) = div(i,j) + u1(i,j+1) - u1(i,j);
for k = 1:n-j-1
div(i,j) = div(i,j) + (-1)^k * gamma(alpha+1) / ...
(gamma(k+1)*gamma(alpha-k+1)) * u1(i,j+k+1);
end
end
end
% y方向
for j = 1:n
for i = 1:m-1
div(i,j) = div(i,j) + u2(i+1,j) - u2(i,j);
for k = 1:m-i-1
div(i,j) = div(i,j) + (-1)^k * gamma(alpha+1) / ...
(gamma(k+1)*gamma(alpha-k+1)) * u2(i+k+1,j);
end
end
end
end
3. 分数阶非线性扩散去噪
% 分数阶非线性扩散去噪
function fractional_diffusion_denoising()
% 清除环境
clear; clc; close all;
% 读取并加入噪声
I = im2double(rgb2gray(imread('lena.png')));
I_noisy = imnoise(I, 'gaussian', 0, 0.03);
% 算法参数
alpha = 0.6; % 分数阶次 (0 < α < 1)
K = 0.05; % 扩散系数阈值
iter = 50; % 迭代次数
dt = 0.1; % 时间步长
% 初始化
u = I_noisy;
% 迭代过程
for n = 1:iter
% 计算分数阶导数
[Gx, Gy] = fractional_gradient(u, alpha);
% 计算梯度幅度
grad_mag = sqrt(Gx.^2 + Gy.^2);
% 计算扩散系数 (Perona-Malik模型)
c = 1 ./ (1 + (grad_mag/K).^2);
% 计算分数阶散度
div = fractional_divergence(c.*Gx, c.*Gy, alpha);
% 更新图像
u = u + dt * div;
% 显示中间结果
if mod(n,10) == 0
figure(1);
subplot(1,2,1); imshow(I_noisy); title('含噪图像');
subplot(1,2,2); imshow(u);
title(['迭代 ',num2str(n),', α=',num2str(alpha),...
', PSNR=',num2str(psnr(u,I)),'dB']);
drawnow;
end
end
% 展示最终成果
figure;
subplot(1,3,1); imshow(I); title('初始图像');
subplot(1,3,2); imshow(I_noisy); title('含噪图像');
subplot(1,3,3); imshow(u); title(['分数阶扩散去噪(α=',num2str(alpha),')']);
% 评价指标
disp(['含噪图像PSNR: ', num2str(psnr(I_noisy,I))]);
disp(['去噪后PSNR: ', num2str(psnr(u,I))]);
disp(['含噪图像SSIM: ', num2str(ssim(I_noisy,I))]);
disp(['去噪后SSIM: ', num2str(ssim(u,I))]);
end
4. 分数阶小波转换去噪
% 分数阶小波转换去噪
function fractional_wavelet_denoising()
% 清除环境
clear; clc; close all;
% 加载并添加噪声
I = im2double(rgb2gray(imread('lena.png')));
I_noisy = imnoise(I, 'gaussian', 0, 0.04);
% 分数阶小波参数
alpha = 0.5; % 分数阶
wavelet = 'db4'; % 小波基
level = 3; % 分解层级
threshold = 0.05; % 阈值
% 分数阶小波变换
[C, S] = frac_wavedec2(I_noisy, level, wavelet, alpha);
% 阈值处理
for l = 1:level
% 水平细节
start_idx = S(1,1)*S(1,2) + 3*sum(S(2:l,1).*S(2:l,2)) + 1;
end_idx = start_idx + S(l+1,1)*S(l+1,2) - 1;
C(start_idx:end_idx) = wthresh(C(start_idx:end_idx), 's', threshold);
% 垂直细节
start_idx = end_idx + 1;
end_idx = start_idx + S(l+1,1)*S(l+1,2) - 1;
C(start_idx:end_idx) = wthresh(C(start_idx:end_idx), 's', threshold);
% 对角细节
start_idx = end_idx + 1;
end_idx = start_idx + S(l+1,1)*S(l+1,2) - 1;
C(start_idx:end_idx) = wthresh(C(start_idx:end_idx), 's', threshold);
end
% 分数阶小波重构
I_denoised = frac_waverec2(C, S, wavelet, alpha);
% 展示成果
figure;
subplot(1,3,1); imshow(I); title('初始图像');
subplot(1,3,2); imshow(I_noisy); title('含噪图像');
subplot(1,3,3); imshow(I_denoised); title(['分数阶小波去噪(α=',num2str(alpha),')']);
% 评价指标
disp(['含噪图像PSNR: ', num2str(psnr(I_noisy,I))]);
disp(['去噪后PSNR: ', num2str(psnr(I_denoised,I))]);
disp(['含噪图像SSIM: ', num2str(ssim(I_noisy,I))]);
disp(['去噪后SSIM: ', num2str(ssim(I_denoised,I))]);
end
% 分数阶小波分解
function [C, S] = frac_wavedec2(x, level, wavelet, alpha)
[m,n] = size(x);
S = zeros(level+2,2);
S(1,:) = [m,n];
% 初始化系数向量
C = [];
% 逐层分解
for k = 1:level
% 分数阶微分预处理
x_frac = fractional_diff(x, alpha);
% 小波分解
[x, h, v, d] = dwt2(x_frac, wavelet);
% 存储细节系数
C = [h(:); v(:); d(:); C];
S(k+1,:) = size(h);
end
% 保存近似系数
C = [x(:); C];
S(level+2,:) = size(x);
end
% 分数阶小波重构
function x = frac_waverec2(C, S, wavelet, alpha)
level = size(S,1)-2;
% 获取近似系数
a_len = prod(S(end,:));
a = reshape(C(1:a_len), S(end,1), S(end,2));
% 逐层重建
for k = level:-1:1
% 获取细节系数
h = reshape(C(a_len+1:a_len+prod(S(k+1,:))), S(k+1,1), S(k+1,2));
v = reshape(C(a_len+prod(S(k+1,:))+1:a_len+2*prod(S(k+1,:))), S(k+1,1), S(k+1,2));
d = reshape(C(a_len+2*prod(S(k+1,:))+1:a_len+3*prod(S(k+1,:))), S(k+1,1), S(k+1,2));
a_len = a_len + 3*prod(S(k+1,:));
% 小波重建
a = idwt2(a, h, v, d, wavelet, S(k,:));
% 分数阶积分后处理
a = fractional_int(a, alpha);
end
x = a;
end
% 分数阶微分
function y = fractional_diff(x, alpha)
[m,n] = size(x);
y = zeros(m,n);
% 基于Grünwald-Letnikov定义
for i = 1:m
for j = 2:n
y(i,j) = x(i,j) - alpha*x(i,j-1);
for k = 2:j-1
y(i,j) = y(i,j) + (-1)^k * gamma(alpha+1) / ...
(gamma(k+1)*gamma(alpha-k+1)) * x(i,j-k);
end
end
end
end
% 分数阶积分
function y = fractional_int(x, alpha)
[m,n] = size(x);
y = zeros(m,n);
% 基于Riemann-Liouville定义
for i = 1:m
for j = 1:n
for k = 0:j-1
y(i,j) = y(i,j) + gamma(alpha) / gamma(alpha+k+1) * x(i,j-k);
end
end
end
end
方法对比与选择指南
分数阶微分掩模法
：
优点：实现简便，运算迅速
缺点：对强烈噪声效果不佳
适用：轻微噪声，需快速处理的情况
分数阶全变分(TV)法
：
优点：边缘保持良好
缺点：可能引发阶梯效应
适用：包含丰富边缘结构的影像
分数阶非线性扩散法
：
优点：能够自动调节扩散程度
缺点：参数选择敏感
适用：纹理复杂的噪声影像
分数阶小波变换法
：
优点：多尺度分析能力出色
缺点：实现较为复杂
适用：同时存在高频和低频噪声的影像
参数选择建议
分数阶次
α\alpha
α
的选择：
针对微分操作：通常选取 0.5-0.8
针对积分操作：通常选取 0.2-0.5
迭代次数：
一般 50-200 次，可根据PSNR变化来确定
正则化参数
λ\lambda
λ
：
通常 0.01-0.2，噪声较强时取较大值
实际运用中，建议通过实验确定最优参数组合，也可使用自适应方法动态调整参数。
京公网安备 11010802022788号
