Code covered by the BSD License  

Highlights from
An edge adaptive directional total variation model

image thumbnail

An edge adaptive directional total variation model

by

 

An edge adaptive directional total variation model is proposed for image noise removal

DirectionalDenoise(x,lam,alp,theta,varargin)
function [y] = DirectionalDenoise(x,lam,alp,theta,varargin)
% minimizes wrt z, 0.5*||x - z||_2^2 + lam*TV_{alpha,theta}(z) (see eq7 in the manuscript)
% Input parameters
% x : noisy image
% lam : lambda parameter in the functional
% alp : length of the major axis of the ellipse
% theta : direction of the ellipse
% varargin : number of iterations (default = 100)
% Output : y : minimizer of the functional above
%
% Ilker Bayram and Mustafa Kamasak
% ibayram@itu.edu.tr, kamasak@itu.edu.tr
% Istanbul Technical University
% March, 2012
% web.itu.edu.tr/ibayram/DTV

% filters for realizing the Delta operator and its transpose
h=[1 -1];
g = h(end:-1:1);% [-1 1]

R = exp(1i*theta);
kappa = alp.^(-2);
kappa = kappa/(8*lam^2);

vx = zeros(size(x)); % vx and vy hold the vector fields described in the manuscript
vy = zeros(size(x));

if isempty(varargin), % number of default iterations
    MAX_ITER = 100;
else
    MAX_ITER = varargin{1};
end

wb = waitbar(0,'Please wait,DTV...');
for iter = 1:MAX_ITER,
    waitbar(iter/MAX_ITER,wb)
    % apply 'A' to v (see eq 13 for the definition of A)
    ux = alp*vx; % apply Lambda
    uy = vy;
    u = R*(ux + 1i*uy); % apply R
    ux = conv2(real(u),g); ux = ux(:,1:end-1); % apply Delta
    uy = conv2(imag(u),g'); uy = uy(1:end-1,:); 
    % subtract from x
    u = x - lam*(ux + uy);

    % now apply A'
    ux = conv2(u,h); ux = ux(:,2:end);
    uy = conv2(u,h'); uy = uy(2:end,:);
    u = lam*conj(R)*(ux + 1i*uy);
    vx = vx + kappa*real(u)*alp;
    vy = vy + kappa*imag(u);
    
    % now apply the threshold
    m = abs(vx + 1i*vy);
    ind = (m > 10^(-10));
    m = max(1,m);
%    m = max(10^(-10),m); % for stability
    vx(ind) = vx(ind)./m(ind);
    vy(ind) = vy(ind)./m(ind);
end
close(wb);

% compute A*v
vx = alp*vx; % apply Lambda
u = R*(vx + 1i*vy); % apply R
ux = conv2(real(u),g); ux = ux(:,1:end-1); % apply Delta
uy = conv2(imag(u),g');  uy = uy(1:end-1,:);
% subtract from x
y = x - lam*(ux + uy);

Contact us