全变分(Total variation),也称为全变差,是图象复原中常用的一个名词。本文简要介绍全变分的概念以及在图象去噪中的应用。
1 一维信号的全变分和去噪
1.1 一维连续函数的全变分
一维连续实函数f(x)f(x)在区间[a,b]⊂R[a,b]⊂R上的全变分定义为参数曲线x→f(x),x∈[a,b]x→f(x),x∈[a,b]的弧长。其表达式为
Vba(f)=∫ba|f′(x)|dx
Vab(f)=∫ab|f′(x)|dx
说白了,所谓的“变分”就是|f(x+Δx)−f(x)||f(x+Δx)−f(x)|,对于连续函数Δx→0Δx→0。而全变分是对函数定义的区间而言的,就是将“变分”在区间上累加起来。
一维离散信号的全变分
从上面连续实函数的全变分,我们可以很容易想到它的离散形式。给出信号序列{yi},i=1,..,n{yi},i=1,..,n,它的全变分定义为
V(y)=∑i=1n|yi+1−yi|
V(y)=∑i=1n|yi+1−yi|
用一句话来概括,全变分是前后项之差的绝对值之和。
1.2 一维信号去噪
当我们得到观察信号xixi,希望xixi变得平滑,也就是对xixi去噪。一种很直观的想法就是让信号的全变分变小。全变分对应的物理意义就是输入信号的平滑度。设得到的恢复信号为yiyi,它应该满足两个条件:
yiyi与观察信号xixi的差距不大。这个差距的常用数学表达式就是
E(x,y)=12∑i(xi−yi)2
E(x,y)=12∑i(xi−yi)2
yiyi的全变分不大。
将物理约束转化为数学模型,求解yy等价于求解下面这个优化问题:
minyE(x,y)+λV(y)
minyE(x,y)+λV(y)
其中参数λλ是正常数,用于调节两个约束的作用大小。注意到E(x,y)E(x,y)和V(y)V(y)都是凸函数,这是一个无约束凸优化问题,有很多经典方法可以求解。
2 二维离散信号(图象)的全变分和去噪
图象是典型的二维离散信号,Rudin在1992年将其全变分定义为
V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2−−−−−−−−−−−−−−−−−−−−−−−√
V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2
这个函数是各项同性的,但是不可微,也并不是凸函数。非凸函数的优化求解难度、速度和稳定性都无法与凸函数相比。因此二维全变分有另一种常用定义
V(y)=∑i,j|yi+1,j−yi,j|2−−−−−−−−−−−√+|yi,j+1−yi,j|2−−−−−−−−−−−√=∑i,j|yi+1,j−yi,j|+|yi,j+1−yi,j|
V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2=∑i,j|yi+1,j−yi,j|+|yi,j+1−yi,j|
这个函数是凸函数了。
仿照一维信号的去噪,基于全变分的图象去噪可以看成求解优化问题
minyE(x,y)+λV(y)
minyE(x,y)+λV(y)
其中E(x,y)E(x,y)作为数据误差项定义为
E(x,y)=12∑i,j(xi,j−yi,j)2
E(x,y)=12∑i,j(xi,j−yi,j)2
当VV有凸函数形式时,问题变为无约束凸优化问题,从而容易求解。
%% This file demonstrates the Split Bregman method for Total Variation denoising % % SB_ATV.m Split Bregman Anisotropic Total Variation Denoising % SB_ITV.m Split Bregman Isotropic Total Variation Denoising clc; clear all; close all; N = 512; n = N^2; f = double(imread('Lena512','png')); g = f(:) + 0.09*max(f(:))*randn(n,1); mu = 20; g_denoise_atv = SB_ATV(g,mu); g_denoise_itv = SB_ITV(g,mu); fprintf('ATV Rel.Err = %g\n',norm(g_denoise_atv(:) - f(:)) / norm(f(:))); fprintf('ITV Rel.Err = %g\n',norm(g_denoise_itv(:) - f(:)) / norm(f(:))); figure; colormap gray; subplot(221); imagesc(f); axis image; title('Original'); subplot(222); imagesc(reshape(g,N,N)); axis image; title('Noisy'); subplot(223); imagesc(reshape(g_denoise_atv,N,N)); axis image; title('Anisotropic TV denoising'); subplot(224); imagesc(reshape(g_denoise_itv,N,N)); axis image; title('Isotropic TV denoising'); function u = SB_ATV(g,mu) % Split Bregman Anisotropic Total Variation Denoising % % u = arg min_u 1/2||u-g||_2^2 + mu*ATV(u) % % g : noisy image % mu: regularisation parameter % u : denoised image % % Refs: % *Goldstein and Osher, The split Bregman method for L1 regularized problems % SIAM Journal on Imaging Sciences 2(2) 2009 % *Micchelli et al, Proximity algorithms for image models: denoising % Inverse Problems 27(4) 2011 % % Benjamin Trémoulhéac % University College London % b.tremoulheac@cs.ucl.ac.uk % April 2012 g = g(:); n = length(g); [B Bt BtB] = DiffOper(sqrt(n)); d = b; u = g; err = 1;k = 1; tol = 1e-3; lambda = 1; while err > tol fprintf('it. %g ',k); [u,~] = cgs(speye(n)+BtB, g-lambda*Bt*(b-d),1e-5,100); Bub = B*u+b; d = max(abs(Bub)-mu/lambda,0).*sign(Bub); b = Bub-d; err = norm(up-u)/norm(u); k = k+1; end fprintf('Stopped because norm(up-u)/norm(u) <= tol=%.1e\n',tol); end function [B Bt BtB] = DiffOper(N) D = spdiags([-ones(N,1) ones(N,1)], [0 1], N,N+1); D(:,1) = []; D(1,1) = 0; B = [ kron(speye(N),D) ; kron(D,speye(N)) ]; Bt = B'; BtB = Bt*B; end
1 matlab版本
2014a
2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社,2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社,2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社,2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社,2015.