How the following code works in image denoising?

function [y] = EADTV(x,lam,alp,costheta,sintheta,varargin)
% Input parameters
% x : noisy image
% lam : lambda parameter in the functional
% alp : length of the major axis of the ellipse
% costheta, sintheta : 
% varargin : number of iterations (default = 100)
% Output : y : minimizer of the functional above

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

R = (costheta + 1i*sintheta);%%%===&&&(((
kappa = alp.^(-2);
kappa = kappa/(8*lam^2);

vx = zeros(size(x)); % vx and vy hold the vector fields 
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,calculating the EADTV...');
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 alpha
    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);

[David Young]

Please format your code by removing the blank lines, selecting the code, and pressing the "{} Code" button. It will make it more readable.

[D Joseph]

Thanks @David. But the answer is not clear for me. I'vl provide the full code. Please consider this question once again

[David Young]

I've edited the answer.

Answer 1

[David Young]

To understand this properly, I think you may need to consult the author of the FEX submission from which it is taken. If you haven't already, you could try messaging him or her. Also, you should read the paper that the submission is based on.

The operation implemented by conv2 is an image differencing operation for finding the gradient of the image. For each pixel, it computes the difference in value between that pixel and one of its neighbours (say, the neighbour to the right). This is finding directional image gradients, as used in the Canny edge detector, but without the smoothing.

I can't say - without doing much more analysis of the code - exactly how this fits in with the A matrix and the function in general. It's fairly complex, and you may need to check carefully that the FEX submission is in fact an accurate implementation of the paper.

I can try to tell what the convolution operation does.

Since g is [-1 1], this line

ux = conv2(real(u),g,'same');

is equivalent to

uReal = real(u);
ux(row, col) = uReal(row, col) - uReal(row, col+1);

for every value of the indices row and col in u. This is like shifting uReal one position to the left, and then subtracting it from the original.

To deal with the problem that we refer to uReal(row, size(u,2)+1) - that is, we need to subtract something from the rightmost column - conv2 pads uReal with a column of zeros on the right.

[D Joseph]

Thanks@ David