%mtf50 = [800 2300 1650] / 8736;     %80mm at f/2 on GFX100
%generate 3 rectangles with those mtf50 and theta=4.5 deg
wbr = 2;                                        %wb multiplier
wbb = 3/2;

r = single(imread('red.png')) / wbr;
g = single(imread('green.png'));
b = single(imread('blue.png')) / wbb;

rggb = zeros(size(g),'single');                 %raw before white balance
rggb(1:2:end,1:2:end) = r(1:2:end,1:2:end);
rggb(1:2:end,2:2:end) = g(1:2:end,2:2:end);
rggb(2:2:end,1:2:end) = g(2:2:end,1:2:end);
rggb(2:2:end,2:2:end) = b(2:2:end,2:2:end);

g = imcrop(g,[135 21 230 176]);                 %raw before white balance
r = imcrop(r,[135 21 230 176]);
b = imcrop(b,[135 21 230 176]);
rggb = imcrop(rggb,[135 21 230 176]);           

L = r/4 *wbr + g/2 + b/4 *wbb;                  %wb, Hardware demosaiced
Y = .250*r *wbr + 1.057*g - .307*b *wbb;        %wb, XYZ D65 demosaiced
                        
%% Get MTF curves for each plane
ow = EdgeCropMTF(rggb,'theta',4.5,'ch',4);      %rggb
oL = EdgeCropMTF(L,'theta',4.5);                %mono
oY = EdgeCropMTF(Y,'theta',4.5);                %mono

or = EdgeCropMTF(r,'theta',4.5);                %mono
og = EdgeCropMTF(g,'theta',4.5);                %mono
ob = EdgeCropMTF(b,'theta',4.5);                %mono

%% Output
figure
plot(ow.mtf.f(1:50),ow.mtf.d(1:50))
grid; hold on;
plot(oL.mtf.f(1:50),oL.mtf.d(1:50),'o')
plot((0:49)/64,( 0.25*or.img.edge.sfr(1:50) + 0.5*og.img.edge.sfr(1:50) +...
    0.25*ob.img.edge.sfr(1:50) )/1,'xk')
title('Hardware System MTF Comparison Bayer CFA')
xlabel('c/p'); ylabel('System MTF'); xlim([0 0.8])
legend('WB Raw Image, Measured','MTF_L = MTF( L = 0.25r + 0.5g + 0.25b )',...
    'MTF_{Sys} = 0.25MTF_r + 0.5MTF_g + 0.25MTF_b')

figure
plot(oY.mtf.f(1:50),oY.mtf.d(1:50),'o'); grid; hold on
plot((0:49)/64,( 0.25*or.img.edge.sfr(1:50) + 1.057*og.img.edge.sfr(1:50) -...
    0.307*ob.img.edge.sfr(1:50) )/1,'xk')
title('Luminance System MTF Comparison Bayer CFA')
xlabel('c/p'); ylabel('System MTF')
legend('MTF_Y = MTF( Y = 0.250r + 1.057g - 0.307b )',...
    'MTF_{Sys} = 0.250MTF_r + 1.057MTF_g - 0.307MTF_b')
xlim([0 0.8])

