
% The code can only be used for research purpose.

% Please cite the following paper when you use it:
    %L. Zhang, X. Wu and D. Zhang 
    %“Color Reproduction from Noisy CFA Data of Single Sensor Digital Cameras,” 
    %IEEE Trans. Image Processing, vol. 16, no. 9, pp. 2184-2197, Sept. 2007.

%Note:
%1. You need Matlab Wavelet Toolbox to run the code.
%2. The code is not optimized and may have bugs. Your revision and improvement are welcome!
%3. In the joint denoising&demosaicking procedure, some boundary areas will not be processed. You may extend the image boundary (e.g. using
%   symmetrical extension) to process the full image. 


%%%1. Simulate the noisy CFA image%%%%%%
clear all;
I=imread('kodak_fence','tif');
I=double(I);
[n,m,ch]=size(I);

figure(1);clf;
imshow(I/255);

%%add noise
noi=randn(n,m);
vr=12;%%noise level in red channel
vb=12;%%noise level in blue channel
vg=12;%%noise level in green channel

In(:,:,3)=I(:,:,3)+vb*noi;
In(:,:,1)=I(:,:,1)+vr*noi;
In(:,:,2)=I(:,:,2)+vg*noi;

%%downsampling to CFA (Bayer pattern) image 
mI(1:n,1:m)=In(:,:,2);
mI(1:2:n,2:2:m)=In(1:2:n,2:2:m,1);
mI(2:2:n,1:2:m)=In(2:2:n,1:2:m,3);%mI is the simulated mosaic noisy CFA image


%%2. Joint denoising and demosaicking%%%%%

Id1=ndmsc(mI,vr,vg,vb);
Id2=den_bayernoise(Id1,vr,vg,vb);

snrf=csnr(Id2,I,25,25) %PSNR of the recovered full color image
figure(2);clf;
imshow(Id2/255);



