【影象去噪】基於matlab全變分演算法影象去噪【含Matlab原始碼 626期】
阿新 • • 發佈：2021-06-24
一、簡介
全變分（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有凸函式形式時，問題變為無約束凸優化問題，從而容易求解。
二、原始碼
function [img_estimated,energyi,ISNRi,energyo,ISNRo]=tvmm_debluring(img_noisy,h,lambda,varargin)
% function [img_estimated]=tvmm_debluring(img_noisy,h,lamdba,...
% 'optional_parameter_name1',value1,'optional_parameter_name2', value2,...);
%
% Total Variation-based image deconvolution with
% a majorization-minimization approach.
%
% Written by: Joao Oliveira, Jose Bioucas-Dias, Mario Figueiredo
% email: [email protected]
% SITE: www.lx.it.pt/~jpaos/tvmm
% Date: 21/11/2005
%
%
% ========================== INPUT PARAMETERS (required) =================
% Parameter Values
% name and description
% ========================================================================
% img_noisy (double) Noisy blured image of size ny.
% h (double) Blur kernel.
% lambda Regularization parameter (which is multiplied by the
% TV penalty).
%
% ======================== OPTIONAL INPUT PARAMETERS ====================
% Parameter Values
% name and description
% =======================================================================
% boa_iter (double) The number of outer loop iterations.
% Default: 20
% cg_iter (double) The max number of inner CG iterations.
% Default: 200
% soim_iter (double) The max number of inner SOIM iterations.
% Default: 20
% soim_cg_iter (double) The max number of CG iterations in inner SOIM.
% Default: 6
% cg_thrld (double) Conjugate Gradient threshold.
% Default: 1e-5;
% image (double) Original image.
% displayIm ({'yes','no'}) if 'yes' is passed the restored imaged
% is displayed along the CG iterations
% Default: 'no'
% info_energyi ({'yes','no'}) if 'yes' is passed the isnr and energy of
% the objective function is displayed inside the inner loop.
% Default: 'no'
% info_energyo ({'yes','no'}) if 'yes' is passed the energy of
% the objective function is displayed at each outter loop.
% Default: 'no'
% info_ISNRi ({'yes','no'}) if 'yes' is passed the isnr is displayed.
% Default: 'no'
% Requires: image
% info_ISNRo ({'yes','no'}) if 'yes' is passed the isnr is displayed.
% Default: 'no'
% Requires: image
% x_0 (double) initial image iteration.
% Default: Wiener filter
% method ({'cg','soim'}) Selects the method used to solve the linear
% system: cg=conjugate gradient; soim=second order iterative
% method.
% Default: 'cg'
% l_min (double) Minimum eigenvalue of the system for the second
% order iterative method.
% Default: 1e-5
% l_max (double) Maximum eigenvalue of the system for the second
% order iterative method.
% Default: 1
%
% ===================================================================
% The following functions can be provided in order to overwrite the
% internal ones. Recall that, by default, all calculations are
% performed with circular convolutions. (see Technical Report)
% ===================================================================
%
% mult_H (function_handle) Function handle to the function that
% performes H*x. This function must accept two parameters:
% x : image to apply the convolution kernel h
% h : Blur kernel
% mult_Ht (function_handle) Function handle to the function that
% performes Ht*x (Ht = H transpose). This function must
% accept two parameters:
% x : image to apply the convolution kernel h
% h : Blur kernel
% diffh (function_handle) Function handle to the function that
% computes the horizontal differences of an image x.
% diffv (function_handle) Function handle to the function that
% computes the vertical differences of an image x.
% diffht (function_handle) Function handle to the function that
% computes the horizontal differences (transposed) of an
% image x.
% diffh (function_handle) Function handle to the function that
% computes the vertical differences (transposed) of an
% image x.
%
%
% ====================== Output parameters ===============================
% img_estimated Estimated image
% energyi Energy of inner loop iterations.
% Required input arguments: 'info_energyi' and 'image'.
% ISNRi Improvement Signal-to-noise ratio of inner loop iterations.
% Required input arguments: 'info_ISNRi' and 'image'.
% energyo Energy of outer loop iterations.
% Required input arguments: 'info_energyo' and 'image'.
% ISNRo Improvement Signal-to-noise ratio of outer loop iterations.
% Required input arguments: 'info_ISNRo' and 'image'.
%
% ============= EXAMPLES =========================================
% Normal use:
% [img_estimated]=tvmm_debluring(img_noisy,lambda,sigma)
%
% Display restored images along CG iterations:
% [img_estimated]=tvmm_debluring(img_noisy,lambda,sigma,...
% 'displayIm','yes')
%
% Perform ISNR calculations of the outer loop iterations:
% [img_estimated,energyi,ISNRi,energyo,ISNRo]=tvmm_debluring(img_noisy,...
% lambda,sigma,'info_SNRo','yes','image',image);
%
% ======================= DEFAULT PARAMETERS ===================
global mH mHt diffh diffv diffht diffvt
boa_iter = 20; % Number of outer loop iterations
cg_iter = 200; % Number of inner CG iterations
soim_iter = 20; % Number of inner SOIM iterations
soim_cg_iter = 6; % Number of CG iterations inside SOIM
cg_thrld = 1e-5; % CG threshold
displayIm = 0;
image=[];
info_energyi='no';
info_energyo='no';
info_ISNRi='no';
info_ISNRo='no';
info_int = 0;
info_ext = 0;
energyi=0;
ISNRi=0;
energyo=0;
ISNRo=0;
mH=@conv2c;
mHt=@conv2c;
diffh=@f_diffh;
diffv=@f_diffv;
diffht=@f_diffht;
diffvt=@f_diffvt;
method='cg';
x_0=[];
l_min=1e-5;
l_max=1;
% ====================== INPUT PARAMETERS ========================
% Test for number of required parametres
if (nargin-length(varargin)) ~= 3
error('Wrong number of required parameters');
end
% Read the optional parameters
if (rem(length(varargin),2)==1)
error('Optional parameters should always go by pairs');
else
for i=1:2:(length(varargin)-1)
% change the value of parameter
switch varargin{i}
case 'boa_iter' % Outer loop iterations
boa_iter = varargin{i+1};
case 'cg_iter' % Inner loop iterations
cg_iter = varargin{i+1};
case 'soim_iter' % Inner loop iterations
soim_iter = varargin{i+1};
case 'soim_cg_iter' % CG iterations in Inner loop
soim_cg_iter = varargin{i+1};
case 'displayIm' % display sucessive xe estimates
if (isequal(varargin{i+1},'yes'))
displayIm=1;
end
case 'cg_thrld' % CG threshold
cg_thrld = varargin{i+1};
case 'image' % Original image
image = varargin{i+1};
case 'info_ISNRi' % display ISNR inside inner loop
info_ISNRi = varargin{i+1};
case 'info_ISNRo' % display ISNR on outer loop
info_ISNRo = varargin{i+1};
case 'info_energyi' % display energy inside inner loop
info_energyi = varargin{i+1};
case 'info_energyo' % display energy on outer loop
info_energyo = varargin{i+1};
case 'x_0' % Initial image iteration
x_0 = varargin{i+1};
case 'mult_H' % Function handle to perform H*x
mH = varargin{i+1};
case 'mult_Ht' % Function handle to perform Ht*x
mHt = varargin{i+1};
case 'diffh' % External function to perform
diffh = varargin{i+1}; % horizontal differences
case 'diffv' % External function to perform
diffv = varargin{i+1}; % vertical differences
case 'diffht' % External function to perform
diffht = varargin{i+1}; % horizontal differences (transposed)
case 'diffvt' % External function to perform
diffvt = varargin{i+1}; % vertical differences (transposed)
case 'method' % External function to perform
method = varargin{i+1}; % vertical differences (transposed)
case 'l_min' % Minimum eigenvalue
l_min = varargin{i+1};
case 'l_max' % Maximum eigenvalue
l_max = varargin{i+1};
otherwise
% Hmmm, something wrong with the parameter string
error(['Unrecognized parameter: ''' varargin{i} '''']);
end;
end;
end
三、執行結果
四、備註
版本：2014a
完整程式碼或代寫加1564658423