% updated March 6th, 2015 by Heng YAO(contact@yaoheng.info)


function [Filter_NLFR, Filter_NLFG, Filter_NLFB, Xmean,Xstd,blen, Valid,minfval]=NoiseEstimation(im)
warning off all

[ny, nx, nc]=size(im);

% divide the image into two non-overlapped regions:  edge region&  non-edge region
display('start edge detection and dilation...');
imEDGE=~edge(rgb2gray(im),'canny');   % edge detecion using Canny operators
se=ones(3,3);
DilatedEDGE=~imdilate(~imEDGE,se);   % edge dilation by three pixels
display('end of edge extraction');

% imgfile = 'denoisedimage.jpg';
% DNim=loadimg(imgfile);

DNim=waveletdenoise(im);   % DNim: denoised image

Nim=double(im)-double(DNim);   % Nim: noise map
% maxStdR=std2(Nim(:,:,1)./255);
% maxStdG=std2(Nim(:,:,2)./255);
% maxStdB=std2(Nim(:,:,3)./255);
maxStdR=0.15;
maxStdG=0.15;
maxStdB=0.15;

blksize=256;
my=ceil(ny/blksize);
mx=ceil(nx/blksize);
display('start image segmentation using Mean-shift Algorithm...');

% hbar = waitbar(0,'Mean Shift...');
MeanImg=im;
T1=150;

Tmean=[]; Tstd=[];  Tnum=[];
  for i=1:my
    for j=1:mx
         s1=(i-1)*blksize+1;
         s2=(j-1)*blksize+1;
        subDNIm=double(DNim(s1:min(s1+blksize-1,ny), s2:min(s2+blksize-1,nx),:));
        subIm=double(im(s1:min(s1+blksize-1,ny), s2:min(s2+blksize-1,nx),:));
        subNim=double(Nim(s1:min(s1+blksize-1,ny), s2:min(s2+blksize-1,nx),:));
        subDilatedEDGE=double(DilatedEDGE(s1:min(s1+blksize-1,ny), s2:min(s2+blksize-1,nx)));
        hs=10;%   hs - spatial bandwith for mean shift analysis
        hr=20;  %   hr - range bandwidth for mean shift analysis
        M=30; %   M  - minimum size of final output regions
        [SegImg Label] = msseg(subDNIm,hs,hr,M);
        [Seg, bmean, bstd, bnum] =MeanStd (subIm, subNim,subDilatedEDGE, Label, T1,maxStdR, maxStdG, maxStdB );
        Tmean=[Tmean bmean];
        Tstd=[Tstd bstd];
        Tnum=[Tnum bnum];
        MeanImg(s1:min(s1+blksize-1,ny), s2:min(s2+blksize-1,nx),:)=Seg;
%         waitbar(((i-1)*mx+j)/(my*mx));
    end
  end
% close(hbar)
display('End of segmentation');

%searching the minimum variance in the same interval
nn=length(Tnum);
h=4; % uniform interval
Tmean1=round(Tmean/h)*h;
[Tmean2, indexmean]=sort(Tmean1,2);
Tstd2=[Tstd(1,indexmean(1,:));Tstd(2,indexmean(2,:)); Tstd(3,indexmean(3,:))];
Tnum2=[Tnum(indexmean(1,:));Tnum(indexmean(2,:)); Tnum(indexmean(3,:))];
Tmean3=Tmean2(:,1);
Tstd3=Tstd2(:,1);
Tnum3=Tnum2(:,1);
blen=zeros(1,3);  % the available sample number in each color component
Tmean4=[Tmean(1,indexmean(1,:));Tmean(2,indexmean(2,:)); Tmean(3,indexmean(3,:))];
Tmean5=Tmean4(:,1);
% 
for bb=1:3
    tt=1;
   for aa=2:nn
          if Tmean2(bb,aa)==Tmean3(bb,tt)
              if Tstd2(bb,aa)<Tstd3(bb,tt)
          Tstd3(bb,tt)=Tstd2(bb,aa);
          Tnum3(bb,tt)=Tnum2(bb,aa);
          Tmean5(bb,tt)=Tmean4(bb,aa);
              end
          else
              tt=tt+1;
              Tmean3(bb,tt)=Tmean2(bb,aa);
              Tstd3(bb,tt)=Tstd2(bb,aa);
              Tnum3(bb,tt)=Tnum2(bb,aa);
              Tmean5(bb,tt)=Tmean4(bb,aa);
          end
   end
   blen(bb)=tt;
end

Xmean=Tmean5;
Xstd=Tstd3;
Xnum=Tnum3;
blen=blen;


% obtain the edge information D2
threshold1= 0.001; % control the color variance
threshold2= 0.4;  % control the color difference
threshold3=0.4;  % control the edge color in the range of two non-edge colors
display('analyzing valid patches for constraint the shape of NLF...');
[Valid ,ColorS1, ColorS2, ColorEdge] = PatchSelect(im, threshold1, threshold2, threshold3);  % searching valid patches
display('ending analyzing...');

save Result.mat ColorS1 ColorS2 ColorEdge Xmean Xstd Xnum blen Tmean Tstd maxStdR maxStdG maxStdB DNim
display('Minilize the objective function...');
load VARC  % obtained for construct the prior model of invCRF coeffiencts C

C1=zeros(9,16);
lb=[MINC 0.005  0 0 0]';
ub=[MAXC  0.04 0.03 0.02 0.02]';
tt=0;
for tt1=0:3
    for tt2=1
    tt=tt+1; 
    display(['the minimum searching for no.',num2str(tt),' group of initial values...']);
    C0=[[0.8 0.4 0.2 0.1 0.05]*tt1  0.018*tt2 0.01*tt2 0.01*tt2 0.01*tt2];
    C1(:,tt) = LMFsolve ('YaoBayesianMap', C0,'MaxIter',200);
    C1(6:9,tt) = abs( C1(6:9,tt) );
    if isreal(C1(:,tt)) & C1(:,tt)>lb & C1(:,tt)<ub                       
    fval(tt)=RefineNLF(C1(:,tt));
    else fval(tt)=1000;
    end
    end
end
[minfval,indexfval]=min(fval);
C=C1(:,indexfval);


EstG=g0+C(1:5)'*H;
DerivativeEstG=diff(EstG)*(length(EstG)-1);
DerivativeEstG=[DerivativeEstG  DerivativeEstG(length(EstG)-1)];

EstSigma1=C(6);
EstSigma2R=C(7);
EstSigma2G=C(8);
EstSigma2B=C(9);
EstNLFR=sqrt(EstSigma1^2*EstG+EstSigma2R^2)./(DerivativeEstG+eps);
EstNLFG=sqrt(EstSigma1^2*EstG+EstSigma2G^2)./(DerivativeEstG+eps);
EstNLFB=sqrt(EstSigma1^2*EstG+EstSigma2B^2)./(DerivativeEstG+eps);
windowsize=64;
Filter_NLFR=filter(ones(1,windowsize)/windowsize,1,EstNLFR);
Filter_NLFG=filter(ones(1,windowsize)/windowsize,1,EstNLFG);
Filter_NLFB=filter(ones(1,windowsize)/windowsize,1,EstNLFB);



display('end optimization');



function NI= waveletdenoise (I)

%  reference:  M K Mihcak, etc. 'SPATIALLY ADAPTIVE STATISTICAL MODELING OF
%  WAVELET IMAGE COEFFICIENTS AND ITS APPLICATION TO DENOISING' 2009 IASSP
%  written by H Yao (contact@yaoheng.info)
% date: 2011-11-14
%===================================================


% I=double(imread('Img.TIF'));

display('start wavelet based image denoising....' );

[mm,nn,dd]=size(I);
bmm=ceil(mm/512);
bnn=ceil(nn/512);
NI=I;
sigma0=5;

for k=1:dd
    for i=1:bmm
        for j=1:bnn
% k=1;
% i=1;
% j=1;
            s1=(i-1)*512+1;
            s2=(j-1)*512+1;
            subI=double(I(s1:min(s1+511,mm), s2:min(s2+511,nn),k));
            [C, S]=wavedec2(subI,4,'db8');
            cA4=appcoef2(C,S,'db8');
            [cH4,cV4,cD4]=detcoef2('all', C, S, 4);
            [cH3,cV3,cD3]=detcoef2('all', C, S, 3);
            [cH2,cV2,cD2]=detcoef2('all', C, S, 2);
            [cH1,cV1,cD1]=detcoef2('all', C, S, 1);
           
            DNcH1=EstimateLocalVar(cH1,sigma0);
            DNcV1=EstimateLocalVar(cV1,sigma0);
            DNcD1=EstimateLocalVar(cD1,sigma0);
            DNcH2=EstimateLocalVar(cH2,sigma0);
            DNcV2=EstimateLocalVar(cV2,sigma0);
            DNcD2=EstimateLocalVar(cD2,sigma0);
            DNcH3=EstimateLocalVar(cH3,sigma0);
            DNcV3=EstimateLocalVar(cV3,sigma0);
            DNcD3=EstimateLocalVar(cD3,sigma0);
            DNcH4=EstimateLocalVar(cH4,sigma0);
            DNcV4=EstimateLocalVar(cV4,sigma0);
            DNcD4=EstimateLocalVar(cD4,sigma0);
            
            % reconstruct the image
            K4=[cA4 DNcH4 DNcV4 DNcD4];
            K3=[DNcH3 DNcV3 DNcD3];
            K2=[DNcH2 DNcV2 DNcD2];
            K1=[DNcH1 DNcV1 DNcD1];
     
            
             DNC=[reshape(K4, 1, S(1,1)*S(1,2)*4),  reshape(K3, 1, S(3,1)*S(3,2)*3),  reshape(K2, 1, S(4,1)*S(4,2)*3),    reshape(K1, 1, S(5,1)*S(5,2)*3) ];
            NI(s1:min(s1+511,mm), s2:min(s2+511,nn),k)=uint8(waverec2(DNC,S, 'db8'));
            
        end
    end
end          

display('End of denoising');

function [S L] = msseg(I,hs,hr,M)
  gray = 0;
  if(size(I,3)==1)
    gray = 1;
    I = repmat(I,[1 1 3]);
  end
  
  if(nargin < 4)
    hs = 10; hr = 7; M = 30;
  end
    
  [fimg labels modes regsize grad conf] = edison_wrapper(I,@RGB2Luv,...
      'SpatialBandWidth',hs,'RangeBandWidth',hr,...
      'MinimumRegionArea',M,'speedup',3);
  L=labels;
  S = Luv2RGB(fimg);

  if(gray == 1)
    S = rgb2gray(S);
  end

function  [Seg, bmean, bstd, bnum] = MeanStd(im, NoiseMap, DilatedEDGE, l, T1,maxStdR, maxStdG, maxStdB) 
% for each segments i, geting samples mean and standard deviation {bmean(i)  bstd(i)}
% input   im: color image
%            NoiseMap: extracted noise map
%            DilatedEDGE: dilated edge map
%            I  : labels segmentated by meanshift algorithm
%            T1: threlshold for filtering the small sample regions.
%output  Seg: image after segmentation
%             bmean: the mean value of each segment
%             bstd: the standard deviation of each segment
%             bnum: the total number of each segment
% written by HengYao(conatct@yaoheng.info)
% date: 2011-10-16
%===============================================================


lmax=max(max(l));
labels = unique( l(:) );
rmean=zeros(3, length(labels) );  % r,g,b 3 channels
rsize=zeros(1, length(labels) );
rstd=zeros(3,length(labels));
% invlabels=zeros(1, lmax, 'int32');
NonNoiseMap=im-NoiseMap;
NoiseMapR=double(NoiseMap(:,:,1));
NoiseMapG=double(NoiseMap(:,:,2));
NoiseMapB=double(NoiseMap(:,:,3));
NonNoiseMapR=double(NonNoiseMap(:,:,1));
NonNoiseMapG=double(NonNoiseMap(:,:,2));
NonNoiseMapB=double(NonNoiseMap(:,:,3));
bmean=zeros(3,length(labels));
bstd=zeros(3,length(labels));
bnum=zeros(1, length(labels));



Seg=zeros(size(im));
segR=Seg(:,:,1);
segG=Seg(:,:,2);
segB=Seg(:,:,3);

% hbar = waitbar(0,'Noise Deviation Extraction...');
t=0;
for i=1:length(labels)
    j=labels(i);
%     invlabels(j)=i;
    b=(l==j);
    rsize(i)=sum(sum(b));
    rmean(1,i)=mean(NonNoiseMapR(b));
    rmean(2,i)=mean(NonNoiseMapG(b));
    rmean(3,i)=mean(NonNoiseMapB(b));
    b1=b & DilatedEDGE;  % filtering the edges
    rsize1=sum(sum(b1));
 
    if rsize1>T1
        rmean(1,i)=mean(NonNoiseMapR(b1));
        rmean(2,i)=mean(NonNoiseMapG(b1));
        rmean(3,i)=mean(NonNoiseMapB(b1));
        rstd(1,i)=std(NoiseMapR(b1)./255);
        rstd(2,i)=std(NoiseMapG(b1)./255);
        rstd(3,i)=std(NoiseMapB(b1)./255);
        if rstd(1,i)<maxStdR && rstd(2,i)<maxStdG && rstd(3,i)<maxStdB 
           t=t+1;
           bmean(:,t)=rmean(:,i);   %sortting the useful regions
           bstd(:,t)=rstd(:,i);
           bnum(t)=rsize1;
        end
    end

    segR=segR+double(b)*rmean(1,i);
    segG=segG+double(b)*rmean(2,i);
    segB=segB+double(b)*rmean(3,i);
    
%     waitbar(i/length(labels))
    
end
% close (hbar)

bmean=bmean(:,1:t);
bstd=bstd(:,1:t);
bnum=bnum(1,1:t);

Seg(:,:,1)=segR;
Seg(:,:,2)=segG;
Seg(:,:,3)=segB;  

function  [Valid ,ColorS1, ColorS2, ColorEdge] = PatchSelect(Img, threshold1, threshold2, threshold3)
% input: Img: an color image
%        threshold1= 0.001; % control the color variance
%        threshold2= 0.4;  % control the color difference
%        threshold3=0.4;  % control the edge color in the range of two non-edge colors
% output: the total number of observed set
% =========================================================================


[mm,nn,cc]=size(Img);
% display('analyzing valid patches for inverse response estimation...');
if cc==1
    display ('pls load a color image');
    return
end

% edge detecting using Canny edge detector 
imEDGE=edge(rgb2gray(Img),'canny');  %0: non-edge pixels  1: edeg pixels 
% dilating the edge path by three pixels
% BW1=bwmorph(imEDGE,'clean');
se=ones(3,3);
dilatedEDGE=imdilate(imEDGE,se);
Img=double(Img)/255;
Valid=zeros(1,2);
% ValidR=zeros(1,2);   % valid patch cordinates to contain monotonic color edges between uniform region colors
% ValidG=zeros(1,2);
% ValidB=zeros(1,2);
T=0;
Color=zeros(1,3,3);
VarP=zeros(1,2,3);
% ColorR=zeros(1,3);   % 1st value: S1's mean color, 2nd value: S2's mean color, 3rd value: edge color
% ColorG=zeros(1,3);
% ColorB=zeros(1,3);
% VarR=zeros(1,2);
% VarG=zeros(1,2);
% VarB=zeros(1,2);

windowsize=15;
for i=8:mm-7
    for j=8:nn-7
        if imEDGE(i,j)==1
            window=dilatedEDGE(i-7:i+7, j-7:j+7);
            windowR=Img(i-7:i+7, j-7:j+7,1);
            windowG=Img(i-7:i+7, j-7:j+7,2);
            windowB=Img(i-7:i+7, j-7:j+7,3);
            [L,num]=bwlabel(~window,4);  %find two partitioned non-edge 
            if num==2
                indexS1=find(L==1);
                indexS2=find(L==2);
                numS1=size(indexS1,1);
                numS2=size(indexS2,1);
                if (numS1>windowsize^2*0.35) && (numS2>windowsize^2*0.35)     % to control the areas of the two partitioned non-edge regions shoud be as close to each other as possilble
                    varS1R=var(windowR(indexS1));
                    varS1G=var(windowG(indexS1));
                    varS1B=var(windowB(indexS1));
                    varS2R=var(windowR(indexS2));
                    varS2G=var(windowG(indexS2));
                    varS2B=var(windowB(indexS2));
                    meanS1R=mean(windowR(indexS1));
                    meanS1G=mean(windowG(indexS1));
                    meanS1B=mean(windowB(indexS1));
                    meanS2R=mean(windowR(indexS2));
                    meanS2G=mean(windowG(indexS2));
                    meanS2B=mean(windowB(indexS2));
                    % every color variance should below a threshold
                    varconstr=(varS1R<threshold1) & (varS1G<threshold1) & (varS1B<threshold1) & (varS2R<threshold1)& (varS2G<threshold1) &(varS2B<threshold1) ;
                   % every edge color must lie within the range of both region
                   if varconstr == 1
                       p1=Img(i,j,1)-min(meanS1R,meanS2R);
                       p2=max(meanS1R,meanS2R)-min(meanS1R,meanS2R);
                       p3=Img(i,j,2)-min(meanS1G,meanS2G);
                       p4=max(meanS1G,meanS2G)-min(meanS1G,meanS2G);
                       p5=Img(i,j,3)-min(meanS1B,meanS2B);
                       p6=max(meanS1B,meanS2B)-min(meanS1B,meanS2B);
                       egdeconstr= ((p1/p2)>threshold3) & ((p1/p2)<(1-threshold3))& ((p3/p4)>threshold3) & ((p3/p4)<(1-threshold3))& ((p5/p6)>threshold3) & ((p5/p6)<(1-threshold3));
                      if egdeconstr ==1 
                    % difference of mean color should above a threshold
                       colordistanceconstr= (abs(meanS1R-meanS2R)>threshold2) &(abs(meanS1G-meanS2G)>threshold2) & (abs(meanS1B-meanS2B)>threshold2);
                       if colordistanceconstr==1
                           X1=abs([Valid(:,1)-i, Valid(:,2)-j]);
                           X2=X1<windowsize;
                           if sum( X2(:,1) & X2(:,2))== 0  % every valid patches should not be overlapped
                               T=T+1;
                               Valid(T,:)=[i, j];
                               Color(T,:,1)=[meanS1R meanS2R Img(i,j,1)];
                               Color(T,:,2)=[meanS1G meanS2G Img(i,j,2)];
                               Color(T,:,3)=[meanS1B meanS2B Img(i,j,3)];
                               VarP(T,:,1)=[varS1R  varS2R];
                               VarP(T,:,2)=[varS1G  varS2G];
                               VarP(T,:,3)=[varS1B  varS2B];
                           end    
                       end
                    
                      end  
                   end
                   
                end
            end
        end
    end
end

if T>100  % constrain the maximum valid patches to 100
    T=100;
    Y1=sum(VarP(:,1,:)+VarP(:,2,:));
    [Z1, indexZ1]=sort(Y1);
    ValidR=Valid(indexZ1(1:100),:,:);
    ColorR=Color(indexZ1(1:100),:,:);
end

ColorS1=[Color(:,1,1) Color(:,1,2) Color(:,1,3)];   % T*{R G B} 
ColorS2=[Color(:,2,1) Color(:,2,2) Color(:,2,3)];
ColorEdge=[Color(:,3,1) Color(:,3,2) Color(:,3,3)];



function  DC=RefineNLF(C)
alpha=C(1:5);
sigma1=C(6);
sigma2r=C(7);
sigma2g=C(8);
sigma2b=C(9);
% [sigma1r sigma1g sigma1b sigma2r sigma2g sigma2c]=C(6:11);
sigma1r=sigma1;
sigma1g=sigma1;
sigma1b=sigma1;


load VARC
load Result

g=g0+alpha'*H;
interval=[0:1/(length(g)-1):1];

lamda1=blen(1);  
lamda2=size(ColorS1,1);
DerivativeG=diff(g)*(length(g)-1);
DerivativeG=[DerivativeG  DerivativeG(length(g)-1)];

NLFR=sqrt(sigma1^2*g+sigma2r^2)./(DerivativeG+eps);
NLFG=sqrt(sigma1^2*g+sigma2g^2)./(DerivativeG+eps);
NLFB=sqrt(sigma1^2*g+sigma2b^2)./(DerivativeG+eps);

factor=1.5;   % parameter 3 $ 1.5
if sum(NLFR>max(Xstd(1,:))*factor) | sum(NLFR<0) | sum(NLFG>max(Xstd(2,:))*factor) |  sum(NLFG<0)  | sum(NLFB>max(Xstd(3,:))*factor) |  sum(NLFB<0)
    DC=1000;
else
    indexXmean=round(Xmean*(length(g)-1)./255);   % size: 3*blen
    sum1=(1/blen(1))*sum((sqrt(g(indexXmean(1,1:blen(1)))*sigma1r^2+sigma2r^2)./(DerivativeG(indexXmean(1,1:blen(1))))-Xstd(1,1:blen(1))).^2+eps);
    sum2=(1/blen(2))*sum((sqrt(g(indexXmean(2,1:blen(2)))*sigma1g^2+sigma2g^2)./(DerivativeG(indexXmean(2,1:blen(2))))-Xstd(2,1:blen(2))).^2+eps);
    sum3=(1/blen(3))*sum((sqrt(g(indexXmean(3,1:blen(3)))*sigma1b^2+sigma2b^2)./(DerivativeG(indexXmean(3,1:blen(3))))-Xstd(3,1:blen(3))).^2+eps);
    DC=lamda1*(sum1+sum2+sum3);
end




function Hden=EstimateLocalVar(H,sigma0)
    W=[3,5,7,9];
    [mmH,nnH]=size(H);
    
    
  Hden=zeros(mmH,nnH);  
 for i=1:mmH
      for j=1:nnH 
          sigmaw=zeros(1,4);
          for k=1:4
            neighbor=W(k);
     
                borderupmm=max(1,i-neighbor);
                borderdownmm=min(mmH,i+neighbor);
                borderupnn=max(1,j-neighbor);
                borderdownnn=min(nnH,j+neighbor);
                AA=H(borderupmm:borderdownmm,borderupnn:borderdownnn);
               [mmAA,nnAA]=size(AA);
               %                  sigmaw(k)=0;
               sigmaw(k)=max(0,sum(sum(AA.^2))/(mmAA*nnAA)-sigma0^2);

          end
            Hden(i,j)=H(i,j)*min(sigmaw)/(min(sigmaw)+sigma0^2);
        end
 end
    

