Robust estimation approach for blind denoising implementation issues.

I've been trying to implement a gaussian denoising script following a paper called 'Robust Estimation Approach for Blind Denoising' by Tamer Rabie. I've been having some issues with the results and I can't find what is the problem, I will focus on the possible causes and maybe if someone is experienced in image processing they can help me:

Script: nothing interesting here.

clc
close all
clear
%profile on
f = imread('peppers.gif');
g = imnoise(f,'gaussian',0,0.016);
N = 4;      %Window maximum size is (2*N+1)*(2*N+1).
C = 0.3;    %Smoothing factor/threshold.
mu = 0.4;   %Over-relaxation constant.
[Image,r] = RobustFilter(g,N,C,mu);
PSNR = 10*log10(255^2/((sum(sum((double(Image)-double(f)).^2)))/(size(f,1)*size(f,2))));
PSNRg = 10*log10(255^2/((sum(sum((double(g)-double(f)).^2)))/(size(f,1)*size(f,2))));
figure
subplot(2,2,1)
imshow(f)
subplot(2,2,2)
imshow(g)
subplot(2,2,3)
imshow(r)
subplot(2,2,4)
imshow(Image)
%profile viewer

Robust Filter: the problem is most certainly here.

%Robust noise filtering.
%g: intensity image.
%N: maximum distance from center pixel.
%C: weighting constant that affects the smoothing amount.
%mu: convergence control.
%Image: denoised image.
%r: found edges of the input image.
function [Image,r] = RobustFilter(g,N,C,mu)
  Sigma = EstimateNoise(g);   %Standard deviation of noise estimation.
  V = LocalVariance(g);       %Local variance for each pixel.
  o = C*sqrt(V)/sqrt(2);      %Outlier rejection parameters.
  [H,W] = size(g);
  Image = uint8(zeros(H,W));
  r = uint8(zeros(H,W));
  for i = 1:H
      for j = 1:W
          Wr = N;
          Wa = N;
          Wl = N;
          Wb = N;
          if(i-Wa < 1); Wa = i-1; end
          if(i+Wb > H); Wb = H-i; end
          if(j-Wl < 1); Wl = j-1; end
          if(j+Wr > W); Wr = W-j; end
          T = V(i,j)-Sigma^2; %Signal noise activity for current pixel.
          if(T > C*Sigma^2)   %If the condition is met, an edge was found.
              Wr = 0;
              Wa = 0;
              Wl = 0;
              Wb = 0;         
              S = 0;          %Obtain optimal size for window.
              while(Wr < N && T > S && (j+Wr) < W)
                  Wr = Wr+1;
                  S = V(i,j+Wr)-Sigma^2;
                  if(T < S)
                      Wr = Wr-1;
                  end
              end
              S = 0;
              while(Wa < N && T > S && (i-Wa) > 1)
                  Wa = Wa+1;
                  S = V(i-Wa,j)-Sigma^2;
                  if(T < S)
                      Wa = Wa-1;
                  end
              end
              S = 0;
              while(Wl < N && T > S && (j-Wl) > 1)
                  Wl = Wl+1;
                  S = V(i,j-Wl)-Sigma^2;
                  if(T < S)
                      Wl = Wl-1;
                  end
              end
              S = 0;
              while(Wb < N && T > S && (i+Wb) < H)
                  Wb = Wb+1;
                  S = V(i+Wb,j)-Sigma^2;
                  if(T < S)
                      Wb = Wb-1;
                  end
              end
              r(i,j) = 255;
          end
          Wn = double(g(i-Wa:i+Wb,j-Wl:j+Wr));    %Gradient based technique
          p = double(g(i,j));                     %to minimize cost function.
          f = p-mu*o(i,j)^2*sum(sum(2*(Wn-p)./(2*o(i,j)^2+(Wn-p).^2)));
          while (abs(f-p) > o(i,j))
              p = f;
              f = p-mu*o(i,j)^2*sum(sum(2*(Wn-p)./(2*o(i,j)^2+(Wn-p).^2)));
          end
          Image(i,j) = f;
      end
  end
end

EstimateNoise: uses the Immerkaer method to estimate gaussian noise standard deviation, no problems here.

%Immerkaer AWGN deviation estimate.
%I: image with AWGN.
%Sigma: noise standard deviation of input image.
function Sigma = EstimateNoise(I)
  [H,W] = size(I);
  I = double(I);
  M = [1 -2 1;-2 4 -2;1 -2 1];
  Sigma = sum(sum(abs(conv2(I,M))));
  Sigma = Sigma*sqrt(0.5*pi)./(6*(W-2)*(H-2));
end

LocalVariance: obtains variance for each pixel in the image, I'm confident in this method but the result is different form squaring every result in stdfilt(), why?

%Returns local variance of matrix.
%M: input matrix. V: local variance for every value in the matrix.
function V = LocalVariance(M)
  M = double(M);
  N = 3;
  W = ones(N)/N.^2;
  u = conv2(M,W,'same');
  U = conv2(M.^2,W,'same');
  V = U-u.^2;
end

Focusing in RobsutFilter.m where the problem probably is:

_ The threshold checking if(T > C*Sigma^2) functions as a local variance edge detector and works pretty well in my opinion, inside there optimal sized windows are calculated for the specified pixel.

_ This filter uses optimal size windowing for each pixel, it looks for edges and limits the window size to reduce statistical errors in the estimation of the real value of the pixel as explained in 'Adaptive hybrid mean and median filtering of high-ISO long-exposure sensor noise for digital photography' by Tamer Rabie. Not using optimal sized windows would result in poor performance but noise should be partially removed, commenting the if(T > C*Sigma^2) section would still reduce noise.

_Then comes the gradient method for finding a global minimum in an error function which uses the Lorentzian influence function. I can't see how to make this one work, if a corrupted pixel is 0 and the values are uint8() the resulting estimation will clip at 0, if the values are double() it will go to negative values which are meaningless in a gray scale image. I understand the gradient method but in this case, where values are limited to [0,255] I can't make this work properly.

Any ideas or suggestions? The results in this paper are amazing. Without the optimal size windows it should still work, so I guess the main issue is at the pixel estimation part. Thanks for your time, I would really appreciate your help.