%Chapter 8,  Fig 8.8
%Eric Dubois, updated 2018-12-14
%Create CFA image, convert to LCC basis, and apply lowpass filtering.
%uses function convert_color_basis and loads LCC_filters.mat
clear all , close all
RGBin = im2double(imread('05.tif'));
figure, imshow(RGBin);
N1 = size(RGBin,1);
N2 = size(RGBin,2);
Z = 0.0;  %set missing components to Z
%downsample the RGB components according to the Bayer pattern
%Upper left pixel is green, pixel to right of it is red, pixel below it is
%blue
%
% G R
% B G
%
g1 = 4; g2=2; g3=4;
RGB_bayer = Z*ones(size(RGBin));
%Write four cosets of RGBin to RGB_bayer
RGB_bayer(1:2:N1,2:2:N2,1) = g1*RGBin(1:2:N1,2:2:N2,1); %red
RGB_bayer(1:2:N1,1:2:N2,2) = g2*RGBin(1:2:N1,1:2:N2,2); %green
RGB_bayer(2:2:N1,2:2:N2,2) = g2*RGBin(2:2:N1,2:2:N2,2); %green
RGB_bayer(2:2:N1,1:2:N2,3) = g3*RGBin(2:2:N1,1:2:N2,3); %blue
%Convert to LCC basis
A_RGB2LCC = [1 2 1; -1 2 -1; -1 0 1]/4;
LCC_bayer = convert_color_basis(RGB_bayer,A_RGB2LCC);
load('LCC_filters');
%Filter the three components of LLL_bayer
LCC_bayer_filt(:,:,1) = imfilter(LCC_bayer(:,:,1),h_11,'symmetric','same');
LCC_bayer_filt(:,:,2) = imfilter(LCC_bayer(:,:,2),h_22,'symmetric','same');
LCC_bayer_filt(:,:,3) = imfilter(LCC_bayer(:,:,3),h_33,'symmetric','same');
%Convert back to sRGB
RGB_bayer_filt= convert_color_basis(LCC_bayer_filt,inv(A_RGB2LCC));
BND = 10;
%compute CMSE and CPSNR
err1 = RGBin(BND+1:N1-BND,BND+1:N2-BND,:) - RGB_bayer_filt(BND+1:N1-BND,BND+1:N2-BND,:);
RMSE1 = sqrt(sum(sum(sum(err1.^2)))/(size(err1,1)*size(err1,2)));
PSNR1 = 20*log10(1/RMSE1);
figure, imshow(RGB_bayer_filt);
%compare with bilinear
RGB_bilinear = im2double(imread('data/05_bilinear.tif'));
err2 = RGBin(BND+1:N1-BND,BND+1:N2-BND,:) - RGB_bilinear(BND+1:N1-BND,BND+1:N2-BND,:);
RMSE2 = sqrt(sum(sum(sum(err2.^2)))/(size(err2,1)*size(err2,2)));
PSNR2 = 20*log10(1/RMSE2);
PSNRdiff = PSNR1-PSNR2;