MATLAB TV去噪程序
基于全变分(Total Variation, TV)的图像去噪MATLAB实现
%% TV去噪程序
% 功能: 实现多种TV模型的图像去噪
% 模型: ROF, LLT, CD, FGP-TV等
clear; clc; close all;
%% 1. 参数设置
% 图像参数
imagePath = 'lena_gray.png'; % 输入图像路径
noiseLevel = 0.1; % 噪声水平 (0-1)
useCustomImage = true; % 使用自定义图像
% TV模型参数
lambda = 0.07; % 正则化参数 (典型值: 0.01-0.2)
modelType = 'ROF'; % 模型类型: 'ROF', 'LLT', 'CD', 'FGP'
iterations = 100; % 迭代次数
stepSize = 0.125; % 梯度下降步长
epsilon = 1e-3; % 平滑参数 (避免除零)
% 高级参数
adaptLambda = true; % 自适应调整lambda
displayInterval = 10; % 显示间隔
saveResults = true; % 保存结果
%% 2. 图像准备
if useCustomImage
% 加载自定义图像
try
noisyImg = imread(imagePath);
if size(noisyImg, 3) == 3
noisyImg = rgb2gray(noisyImg);
end
noisyImg = im2double(noisyImg);
catch
warning('无法加载图像，使用默认测试图像');
noisyImg = phantom(256);
noisyImg = im2double(noisyImg);
end
else
% 生成测试图像
[x, y] = meshgrid(linspace(-1, 1, 256));
noisyImg = exp(-5*(x.^2 + y.^2));
noisyImg = noisyImg/max(noisyImg(:));
end
% 添加高斯噪声
noisyImg = imnoise(noisyImg, 'gaussian', 0, noiseLevel);
% 显示原始噪声图像
figure('Name', '输入图像', 'NumberTitle', 'off');
subplot(121), imshow(noisyImg), title('含噪图像');
subplot(122), imhist(noisyImg), title('噪声图像直方图');
%% 3. TV去噪核心算法
% 初始化
denoisedImg = noisyImg;
[M, N] = size(noisyImg);
% 自适应调整lambda
if adaptLambda
lambda = estimateOptimalLambda(noisyImg, modelType);
fprintf('自适应调整lambda: %.4f\n', lambda);
end
% 选择模型
switch lower(modelType)
case 'rof'
denoisedImg = ROF_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'llt'
denoisedImg = LLT_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'cd'
denoisedImg = CD_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'fgp'
denoisedImg = FGP_TV_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
otherwise
error('未知模型类型: %s', modelType);
end
% 显示结果
figure('Name', '去噪结果', 'NumberTitle', 'off', 'Position', [100, 100, 1200, 500]);
subplot(131), imshow(noisyImg), title('含噪图像');
subplot(132), imshow(denoisedImg), title(sprintf('%s去噪结果', modelType));
subplot(133), imshowpair(noisyImg, denoisedImg, 'montage'), title('对比');
%% 4. 结果评估
% 计算PSNR
psnrNoisy = psnr(noisyImg, denoisedImg);
psnrClean = psnr(idealImg, denoisedImg); % 需要原始干净图像
% 计算SSIM
ssimNoisy = ssim(denoisedImg, noisyImg);
ssimClean = ssim(denoisedImg, idealImg); % 需要原始干净图像
% 计算RMSE
rmse = sqrt(mean((denoisedImg(:) - noisyImg(:)).^2));
fprintf('===== 评估结果 =====\n');
fprintf('模型: %s\n', modelType);
fprintf('噪声水平: %.2f\n', noiseLevel);
fprintf('正则化参数λ: %.4f\n', lambda);
fprintf('PSNR: %.2f dB\n', psnrNoisy);
fprintf('SSIM: %.4f\n', ssimNoisy);
fprintf('RMSE: %.4f\n', rmse);
%% 5. 参数敏感性分析
% 扫描lambda参数
lambdaRange = linspace(0.01, 0.2, 10);
psnrValues = zeros(size(lambdaRange));
ssimValues = zeros(size(lambdaRange));
for i = 1:length(lambdaRange)
tempImg = ROF_Denoising(noisyImg, lambdaRange(i), iterations, stepSize, epsilon);
psnrValues(i) = psnr(tempImg, noisyImg);
ssimValues(i) = ssim(tempImg, noisyImg);
end
% 绘制参数敏感性曲线
figure('Name', '参数敏感性分析', 'NumberTitle', 'off');
subplot(211), plot(lambdaRange, psnrValues, 'bo-');
xlabel('正则化参数λ');
ylabel('PSNR (dB)');
title('PSNR随λ的变化');
grid on;
subplot(212), plot(lambdaRange, ssimValues, 'rs-');
xlabel('正则化参数λ');
ylabel('SSIM');
title('SSIM随λ的变化');
grid on;
%% 6. 模型比较
models = {'ROF', 'LLT', 'CD', 'FGP'};
psnrComparison = zeros(1, length(models));
ssimComparison = zeros(1, length(models));
for i = 1:length(models)
switch models{i}
case 'ROF'
result = ROF_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'LLT'
result = LLT_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'CD'
result = CD_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
case 'FGP'
result = FGP_TV_Denoising(noisyImg, lambda, iterations, stepSize, epsilon);
end
psnrComparison(i) = psnr(result, noisyImg);
ssimComparison(i) = ssim(result, noisyImg);
end
% 绘制模型比较图
figure('Name', '模型比较', 'NumberTitle', 'off');
subplot(211), bar(psnrComparison);
set(gca, 'XTickLabel', models);
ylabel('PSNR (dB)');
title('不同模型的PSNR比较');
grid on;
subplot(212), bar(ssimComparison);
set(gca, 'XTickLabel', models);
ylabel('SSIM');
title('不同模型的SSIM比较');
grid on;
%% 7. 高级TV模型实现
% ROF模型 (Rudin-Osher-Fatemi)
function u = ROF_Denoising(f, lambda, maxIter, tau, eps)
[M, N] = size(f);
u = f; % 初始化
for iter = 1:maxIter
% 计算梯度
ux = [diff(u, 1, 2), zeros(M, 1)];
uy = [diff(u, 1, 1); zeros(1, N)];
% 计算梯度模
gradMag = sqrt(ux.^2 + uy.^2 + eps);
% 计算散度
px = ux ./ gradMag;
py = uy ./ gradMag;
div_p = [px(:,1), diff(px, 1, 2), -px(:,end)] + ...
[py(1,:); diff(py, 1, 1); -py(end,:)];
% 更新图像
u = u + tau * (lambda * (f - u) + div_p);
% 显示进度
if mod(iter, displayInterval) == 0
fprintf('ROF: Iteration %d/%d\n', iter, maxIter);
end
end
end
% LLT模型 (Local Laplacian Filter)
function u = LLT_Denoising(f, lambda, maxIter, tau, eps)
[M, N] = size(f);
u = f;
for iter = 1:maxIter
% 计算局部拉普拉斯算子
laplacian = del2(u);
% 更新图像
u = u + tau * (lambda * (f - u) - laplacian);
% 显示进度
if mod(iter, displayInterval) == 0
fprintf('LLT: Iteration %d/%d\n', iter, maxIter);
end
end
end
% CD模型 (Curvature Diffusion)
function u = CD_Denoising(f, lambda, maxIter, tau, eps)
[M, N] = size(f);
u = f;
for iter = 1:maxIter
% 计算梯度
ux = [diff(u, 1, 2), zeros(M, 1)];
uy = [diff(u, 1, 1); zeros(1, N)];
% 计算曲率
curvature = (uxx.*uyy - uxy.^2) ./ (ux.^2 + uy.^2 + eps).^(3/2);
% 更新图像
u = u + tau * (lambda * (f - u) + curvature);
% 显示进度
if mod(iter, displayInterval) == 0
fprintf('CD: Iteration %d/%d\n', iter, maxIter);
end
end
end
% FGP-TV模型 (Fast Gradient Projection)
function u = FGP_TV_Denoising(f, lambda, maxIter, tau, eps)
[M, N] = size(f);
u = f;
for iter = 1:maxIter
% 计算梯度
ux = [diff(u, 1, 2), zeros(M, 1)];
uy = [diff(u, 1, 1); zeros(1, N)];
% 计算梯度模
gradMag = sqrt(ux.^2 + uy.^2 + eps);
% 投影梯度
proj = max(0, 1 - 1./(lambda*gradMag));
px = ux .* proj;
py = uy .* proj;
% 计算散度
div_p = [px(:,1), diff(px, 1, 2), -px(:,end)] + ...
[py(1,:); diff(py, 1, 1); -py(end,:)];
% 更新图像
u = u + tau * (div_p - lambda * f);
% 显示进度
if mod(iter, displayInterval) == 0
fprintf('FGP-TV: Iteration %d/%d\n', iter, maxIter);
end
end
end
%% 8. 辅助函数
% 估计最优lambda
function lambda_opt = estimateOptimalLambda(noisyImg, modelType)
% 使用L曲线法估计最优lambda
lambdaRange = logspace(-3, 0, 10);
residualNorm = zeros(size(lambdaRange));
regularizationNorm = zeros(size(lambdaRange));
for i = 1:length(lambdaRange)
tempImg = ROF_Denoising(noisyImg, lambdaRange(i), 50, 0.1, 1e-3);
residualNorm(i) = norm(tempImg(:) - noisyImg(:), 2);
regularizationNorm(i) = norm(gradient(tempImg), 'fro');
end
% 寻找曲率最大点
curvature = computeCurvature(lambdaRange, residualNorm);
[~, idx] = max(curvature);
lambda_opt = lambdaRange(idx);
end
% 计算曲率
function curvature = computeCurvature(x, y)
dx = gradient(x);
dy = gradient(y);
ddx = gradient(dx);
ddy = gradient(dy);
curvature = abs(dx .* ddy - dy .* ddx) ./ (dx.^2 + dy.^2).^(3/2);
end
% 计算SSIM
function ssimValue = ssim(img1, img2)
K = [0.01, 0.03];
L = 1; % 动态范围
C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
C3 = C2/2;
mu1 = filter2(fspecial('average', 11), img1, 'same');
mu2 = filter2(fspecial('average', 11), img2, 'same');
mu1_sq = mu1.^2;
mu2_sq = mu2.^2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(fspecial('average', 11), img1.^2, 'same') - mu1_sq;
sigma2_sq = filter2(fspecial('average', 11), img2.^2, 'same') - mu2_sq;
sigma12 = filter2(fspecial('average', 11), img1.*img2, 'same') - mu1_mu2;
luminance = (2*mu1_mu2 + C1) ./ (mu1_sq + mu2_sq + C1);
contrast = (2*sqrt(sigma1_sq.*sigma2_sq) + C2) ./ (sigma1_sq + sigma2_sq + C2);
structure = (sigma12 + C3) ./ (sqrt(sigma1_sq.*sigma2_sq) + C3);
ssimValue = mean(luminance .* contrast .* structure, 'all');
end
%% 9. 可视化分析
% 残差分析
residual = denoisedImg - noisyImg;
figure('Name', '残差分析', 'NumberTitle', 'off');
subplot(221), imshow(residual, []), title('残差图像');
subplot(222), histogram(residual, 50), title('残差直方图');
subplot(223), imshow(log(abs(fftshift(fft2(residual)))), []), title('残差频谱');
subplot(224), plot(sort(residual(:))), title('残差排序');
% 梯度分析
[Gx, Gy] = gradient(denoisedImg);
gradientMagnitude = sqrt(Gx.^2 + Gy.^2);
figure('Name', '梯度分析', 'NumberTitle', 'off');
subplot(131), imshow(Gx, []), title('水平梯度');
subplot(132), imshow(Gy, []), title('垂直梯度');
subplot(133), imshow(gradientMagnitude, []), title('梯度幅值');
% 边缘保持分析
edgesNoisy = edge(noisyImg, 'canny');
edgesDenoised = edge(denoisedImg, 'canny');
figure('Name', '边缘保持分析', 'NumberTitle', 'off');
subplot(131), imshow(noisyImg), title('含噪图像');
subplot(132), imshow(edgesNoisy), title('噪声图像边缘');
subplot(133), imshow(edgesDenoised), title('去噪图像边缘');
%% 10. 高级应用
% 彩色图像去噪
function denoisedColor = tvDenoiseColor(colorImg, lambda, modelType)
% 转换到Lab颜色空间
labImg = rgb2lab(colorImg);
L = labImg(:,:,1)/100;
a = labImg(:,:,2)/120;
b = labImg(:,:,3)/120;
% 对每个通道去噪
L_denoised = ROF_Denoising(L, lambda, 100, 0.1, 1e-3);
a_denoised = ROF_Denoising(a, lambda*0.8, 100, 0.1, 1e-3);
b_denoised = ROF_Denoising(b, lambda*0.8, 100, 0.1, 1e-3);
% 合并通道
labDenoised = cat(3, L_denoised*100, a_denoised*120, b_denoised*120);
denoisedColor = lab2rgb(labDenoised);
end
% 视频去噪 (帧处理)
function denoisedVideo = tvDenoiseVideo(videoPath, lambda)
videoReader = VideoReader(videoPath);
denoisedVideo = VideoWriter('denoised_video.mp4', 'MPEG-4');
denoisedVideo.FrameRate = videoReader.FrameRate;
open(denoisedVideo);
while hasFrame(videoReader)
frame = readFrame(videoReader);
if size(frame, 3) == 3
denoisedFrame = tvDenoiseColor(frame, lambda, 'ROF');
else
denoisedFrame = ROF_Denoising(im2double(frame), lambda, 50, 0.1, 1e-3);
end
writeVideo(denoisedVideo, denoisedFrame);
end
close(denoisedVideo);
end
%% 11. 保存结果
if saveResults
% 保存去噪图像
imwrite(denoisedImg, sprintf('denoised_%s.png', modelType));
% 保存评估结果
results = struct('Model', modelType, ...
'NoiseLevel', noiseLevel, ...
'Lambda', lambda, ...
'PSNR', psnrNoisy, ...
'SSIM', ssimNoisy, ...
'RMSE', rmse);
save('tv_denoising_results.mat', 'results');
% 保存中间结果
save('intermediate_results.mat', 'noisyImg', 'denoisedImg', 'lambdaRange', ...
'psnrValues', 'ssimValues');
end
%% 12. 性能优化技巧
% 使用GPU加速
function u = ROF_GPU(f, lambda, maxIter, tau, eps)
if gpuDeviceCount > 0
f_gpu = gpuArray(f);
u_gpu = f_gpu;
for iter = 1:maxIter
% GPU加速计算
ux = [diff(u_gpu, 1, 2), zeros(size(u_gpu,1), 1, 'gpuArray')];
uy = [diff(u_gpu, 1, 1); zeros(1, size(u_gpu,2), 'gpuArray')];
gradMag = sqrt(ux.^2 + uy.^2 + eps);
px = ux ./ gradMag;
py = uy ./ gradMag;
div_p = [px(:,1), diff(px, 1, 2), -px(:,end)] + ...
[py(1,:); diff(py, 1, 1); -py(end,:)];
u_gpu = u_gpu + tau * (lambda * (f_gpu - u_gpu) + div_p);
end
u = gather(u_gpu);
else
u = ROF_Denoising(f, lambda, maxIter, tau, eps);
end
end
% 多分辨率处理
function u = multiResTV(f, lambda, levels)
% 金字塔下采样
pyramid = cell(levels, 1);
pyramid{1} = f;
for i = 2:levels
pyramid{i} = imresize(pyramid{i-1}, 0.5);
end
% 从粗到细处理
u = pyramid{levels};
for i = levels:-1:1
if i < levels
u = imresize(u, size(pyramid{i}));
end
u = ROF_Denoising(u, lambda/levels, 20, 0.1, 1e-3);
end
end
TV去噪原理与模型
1. 全变分去噪基本原理
全变分去噪基于变分原理，通过最小化能量泛函来去除噪声：
E(u) = ∫_Ω |∇u| dx + λ/2 ∫_Ω (u - f)^2 dx
其中：
u是去噪后的图像f是含噪图像λ是正则化参数|∇u|是梯度的L1范数（全变分）
2. 主要TV模型比较
| 模型 | 全称 | 特点 | 适用场景 |
|---|---|---|---|
| ROF | Rudin-Osher-Fatemi | 基础TV模型，保留边缘好 | 通用去噪 |
| LLT | Local Laplacian Filter | 局部自适应平滑 | 纹理保留 |
| CD | Curvature Diffusion | 基于曲率的扩散 | 医学图像 |
| FGP-TV | Fast Gradient Projection | 快速投影算法 | 实时处理 |
3. 数值求解方法
- 梯度下降法：简单直观，但收敛慢
- Chambolle对偶法：高效精确，适合ROF模型
- Split-Bregman：处理复杂约束，速度快
- Proximal梯度法：处理非光滑问题
关键功能模块
1. 核心去噪算法
function u = ROF_Denoising(f, lambda, maxIter, tau, eps)
u = f; % 初始化
for iter = 1:maxIter
% 计算梯度
ux = [diff(u, 1, 2), zeros(M, 1)];
uy = [diff(u, 1, 1); zeros(1, N)];
% 计算梯度模
gradMag = sqrt(ux.^2 + uy.^2 + eps);
% 计算散度
px = ux ./ gradMag;
py = uy ./ gradMag;
div_p = computeDivergence(px, py);
% 更新图像
u = u + tau * (lambda * (f - u) + div_p);
end
end
2. 参数优化模块
% 自适应调整lambda
function lambda_opt = estimateOptimalLambda(noisyImg, modelType)
lambdaRange = logspace(-3, 0, 10);
for i = 1:length(lambdaRange)
tempImg = ROF_Denoising(noisyImg, lambdaRange(i), 50, 0.1, 1e-3);
residualNorm(i) = norm(tempImg(:) - noisyImg(:), 2);
regularizationNorm(i) = norm(gradient(tempImg), 'fro');
end
% 寻找L曲线拐点
curvature = computeCurvature(lambdaRange, residualNorm);
[~, idx] = max(curvature);
lambda_opt = lambdaRange(idx);
end
3. 性能评估模块
% 计算PSNR
psnrValue = psnr(denoisedImg, noisyImg);
% 计算SSIM
ssimValue = ssim(denoisedImg, noisyImg);
% 计算RMSE
rmse = sqrt(mean((denoisedImg(:) - noisyImg(:)).^2));
4. 高级处理模块
% 彩色图像去噪
function denoisedColor = tvDenoiseColor(colorImg, lambda, modelType)
labImg = rgb2lab(colorImg);
L = labImg(:,:,1)/100;
a = labImg(:,:,2)/120;
b = labImg(:,:,3)/120;
L_denoised = ROF_Denoising(L, lambda, 100, 0.1, 1e-3);
a_denoised = ROF_Denoising(a, lambda*0.8, 100, 0.1, 1e-3);
b_denoised = ROF_Denoising(b, lambda*0.8, 100, 0.1, 1e-3);
labDenoised = cat(3, L_denoised*100, a_denoised*120, b_denoised*120);
denoisedColor = lab2rgb(labDenoised);
end
参数选择与优化
1. 正则化参数λ
- 过小：去噪不足，残留噪声
- 过大：过度平滑，丢失细节
- 选择策略： 经验值：0.01-0.2 自适应估计：
estimateOptimalLambda()L曲线法：平衡残差和正则化项
2. 迭代次数
- 过少：收敛不充分
- 过多：计算开销大，可能引入伪影
- 建议值：50-200次
3. 步长τ
- 过大：发散不收敛
- 过小：收敛缓慢
- 建议值：0.05-0.25
4. 平滑参数ε
- 过小：数值不稳定
- 过大：边缘模糊
- 建议值：1e-3-1e-2
参考代码 TV去噪程序 www.youwenfan.com/contentcnp/83418.html
扩展功能
1. 非局部TV去噪
function u = NLTV_Denoising(f, lambda, patchSize, searchWindow)
% 非局部均值思想与TV结合
% 1. 计算图像块相似度
% 2. 加权平均TV去噪
end
2. 复合噪声处理
function u = mixedNoiseTV(f, lambdaImpulse, lambdaGaussian)
% 脉冲噪声检测
impulseMask = detectImpulseNoise(f);
% 分别处理
u = f;
u(impulseMask) = medianFilter(f(impulseMask)); % 中值滤波处理脉冲噪声
u = ROF_Denoising(u, lambdaGaussian); % TV处理高斯噪声
end
3. 盲去卷积
function [u, kernel] = blindDeconvTV(f, lambda, maxIter)
% 交替优化图像和模糊核
u = f;
kernel = ones(5,5)/25; % 初始核
for iter = 1:maxIter
% 固定核，更新图像
u = ROF_Denoising(deconvblind(f, kernel), lambda);
% 固定图像，更新核
kernel = updateKernel(f, u);
end
end
4. 深度学习结合
function u = deepTV_Denoising(f)
% 预训练CNN提取特征
features = cnnFeatureExtractor(f);
% TV正则化优化
u = ROF_Denoising(features, 0.05);
% 后处理
u = guidedFilter(u, f, 5, 0.1);
end
浙公网安备 33010602011771号