全变分模型(TV)

网上讲全变分模型的博客不少，但如果把技术演变过程讲出，能更好地理解全变分模型。希望对大家有用。

图像复原问题

全变分模型(TV)

传统的图像复原方法多建立线性系统的基础上。图像复原中另一个重要的经典方法称为有约束的最小二乘法。其出发点是：输出图像u在满足条件

算法推导过程

总结：TV项用在图像复原和去噪中的作用就是保持图像的光滑性，消除图像复原可能带来的伪影。缺陷在于会使得复原的图像过于光滑，在一些细节比较多的图像中会使得复原后的图像丢失细节。

实例演示

主函数

 %% Restore an degraded Image (noising and blurring) by TV method. 
 %% The numerical implementation is based on Roe upwind schema
 
clear all;
close all;
clc;

Img = imread('lenna.bmp');
Img = double(Img);
figure(1); imshow(uint8(Img)); title('Origin');
[nrow, ncol] = size(Img);
% blurring image
Img = gauss(Img,7,3);
% add noise to the blurred image
I0=Img+10*randn([nrow, ncol]);
figure(2); imshow(uint8(I0));title('Noise');


lamda=0.2;        % This is a important parameter to get a better result
timestep=0.01;

I_temp=I0;        % Initialization

%Iteration begin
for n=1:600
   
       Ix = 0.5*(I_temp(:,[2:ncol,ncol])-I_temp(:,[1,1:ncol-1]));
       Iy = 0.5*(I_temp([2:nrow,nrow],:)-I_temp([1,1:nrow-1],:));
       Ix_back = I_temp-I_temp(:,[1,1:ncol-1]);
       Ix_forw = I_temp(:,[2:ncol,ncol])-I_temp;
       Iy_back = I_temp-I_temp([1,1:nrow-1],:);
       Iy_forw = I_temp([2:nrow,nrow],:)-I_temp;
       grad = Ix.^2+Iy.^2+eps;
       
       Ixx = I_temp(:,[2:ncol,ncol])+I_temp(:,[1,1:ncol-1])-2*I_temp;
       Iyy = I_temp([2:nrow,nrow],:)+I_temp([1,1:nrow-1],:)-2*I_temp;
       Ixy = 0.25*(I_temp([2:nrow,nrow],[2:ncol,ncol])+I_temp([1,1:nrow-1],[1,1:ncol-1])-...
           I_temp([2:nrow,nrow],[1,1:ncol-1])-I_temp([1,1:nrow-1],[2:ncol,ncol]));
       term1 = (Ixx.*Iy.^2-2*Ix.*Iy.*Ixy+Iyy.*Ix.^2)./grad;
       
       D=gauss(I_temp,7,3)-I0;
       DD = gauss(D,7,3);
       upwind_x=Ix_back;
       upwind_x(Ix.*DD<0)=Ix_forw(Ix.*DD<0);
       upwind_y=Iy_back;
       upwind_y(Iy.*DD<0)=Iy_forw(Iy.*DD<0);
       term2 = lamda*sqrt(upwind_x.^2+upwind_y.^2).*DD;
       I_t = term1-term2;

       I_temp=I_temp+timestep*I_t;
end
figure(3);imshow(uint8(I_temp));title('TV');

gauss.m

function Ig=gauss(I,ks,sigma2)
%%% 子函数: gauss()实现高斯平滑滤波
%%% 参数说明：
%%% I       - 待平滑图像
%%% ks      - 高斯核大小 (奇数)
%%% sigma2  - 高斯函数的方差
%%% Ig      - 返回的高斯平滑后的图象

[Ny,Nx]=size(I);
hks=(ks-1)/2;  % 高斯核的一半
%%- 一维卷积
if (Ny<ks)
    x=(-hks:hks);
    flt=exp(-(x.^2)/(2*sigma2));       % 一维高斯函数
    flt=flt/sum(sum(flt));             % 归一化
    %%- 扩展
    x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx));
    eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)];
    Ig=conv(eI,flt);
    Ig=Ig(:,ks+hks+1:Nx+ks+hks);       % 截尾
else
    %%- 二维卷积
    x=ones(ks,1)*(-hks:hks); y=x';
    flt=exp(-(x.^2+y.^2)/(2*sigma2));  % 二维高斯函数
    flt=flt/sum(sum(flt));             % 归一化
    %%- 扩展
    if (hks>1)
        xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')';
    else
        xL=I(:,1); xR=I(:,Nx);
    end
    eI=[xL*ones(1,hks) I xR*ones(1,hks)];
    if (hks>1)
        xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:));
    else
        xU=eI(1,:); xD=eI(Ny,:);
    end
    eI=[ones(hks,1)*xU; eI; ones(hks,1)*xD];
    Ig=conv2(eI,flt,'valid');
end

运行结果