function [ handles ] = genTestImages( handles )
% Generates test images and their reconstructions from either the 
% hyperReflectanceData.mat file or synthetic resolution targets
%

chanCenters = handles.filterCentersVal(handles.fourierFiltInd,:);

chanXiStart = min(chanCenters(:,1));
chanXiStop = max(chanCenters(:,1));

chanEtaStart = min(chanCenters(:,2));
chanEtaStop = max(chanCenters(:,2));

numChannels = size(chanCenters,1);
chanHalfPlaneIdx = round(numChannels/2);

%Create edge check vector.
isXiEdge = ((chanCenters(:,1) == -0.5) | (chanCenters(:,1) == 0.5));
isEtaEdge = ((chanCenters(:,2) == -0.5) | (chanCenters(:,2) == 0.5));
isEdge = isXiEdge + isEtaEdge;
xiIndStart = find(handles.filterCentersVal(:,1) == chanXiStart);
xiIndStop = find(handles.filterCentersVal(:,1) == chanXiStop);
etaIndStart = find(handles.filterCentersVal(:,2) == chanEtaStart);
etaIndStop = find(handles.filterCentersVal(:,2) == chanEtaStop);

% The below doesn't work correctly yet for odd pixel sizes, must be even.
% Will be fixed in a future release
xiT = [(handles.xiC(1:end-1) - 1) handles.xiC (handles.xiC(2:end-1)+1)];
etaT = [(handles.etaC(1:end-1) - 1) handles.etaC (handles.etaC(2:end-1)+1)];

for k = 1:numChannels
    indXT{k} = find(abs(xiT - chanCenters(k,1)) < 10^-10);
    indYT{k} = find(abs(etaT - chanCenters(k,2)) < 10^-10);
end

indRadX = round(handles.numXpixels/2);
indRadY = round(handles.numYpixels/2);

[X,Y] = meshgrid(xiT,etaT);

%theta = linspace(0,pi/2,handles.numFiltersD);

imgTot = zeros(handles.numYpixels,handles.numXpixels);

for k = 1:handles.numFiltersD
    if handles.useSynthTestImg
%         Z{k} = imgFunc2D(handles.numXpixels,handles.numYpixels,theta(k));
        Z{k} = imgFunc2D(handles.numXpixels,handles.numYpixels,pi/4 + k*pi/64);
    else
        if isempty(handles.hyperReflectanceData)
            load('hyperReflectanceData');
            handles.hyperReflectanceData = hyperReflectanceData;
        end
        %Crops or expands (doesn't scale) the hyperspectral images to be the same as the sensor size
        Z{k} = imgFromHyperSpecReflectance(handles.numXpixels,handles.numYpixels,handles.hyperReflectanceData,k);
    end
    baseImgMean = mean(mean(Z{k}));
    baseImgMax = max(max(Z{k}));
    baseImgMin = min(min(Z{k}));
    if baseImgMean > 0
        baseMax(k) = baseImgMax;
        baseMin(k) = baseImgMin;
        baseMean(k) = baseImgMean;
    else
        baseMax(k) = 1;
        baseMin(k) = 1;
        baseMean(k) = 1;
    end
    imgTot = imgTot + handles.maskCube(:,:,k).*Z{k};
    imgTransFactor(k) = mean(mean(handles.maskCube(:,:,k))); %this should always be < 1 for passive filters.     
end

%set baseband to actual baseband passed.
newLf = handles.Lf(:,~all(abs(handles.Lf)==0));    
newLf = newLf./sum(newLf(:,chanHalfPlaneIdx));
newLf = newLf.*sum(imgTransFactor);
newLfSize = size(newLf);

newLfInv = pinv(newLf');
newLfInv(abs(newLfInv) < 10^-12) = 0;
[numRows,numCols] = size(newLfInv);

imgMeanAfterMask = mean(mean(imgTot));
imgMean = mean(baseMean);
% imgMult = 1;
imgMult = handles.numPElectrons/imgMeanAfterMask;

imgTot = imgMult.*imgTot;
imgMax = max(max(imgTot));
imgMin = min(min(imgTot));
    
imgTotN = imgTot;

% Generate Fourier domain filters at the proper locations
for k = 1:numChannels
    F{k} = planckTaper2D(handles.radius(k),handles.eps(k),X-chanCenters(k,1),Y-chanCenters(k,2));    
end

if handles.shotNoise
    imgTotN = fPoissrnd( imgTot, 1000 );
    snrIndB = 10*log10(handles.numPElectrons/(sqrt((handles.numPElectrons*handles.noiseMult)^2 + handles.numPElectrons)))       
else    
    snrIndB = 10*log10(1/(handles.noiseMult));   
end
snrStr = ['SNR_' num2str(snrIndB,'%5.1f') 'dB'];

%thermal and readout noise round(sqrt(s2).*randn(size(N)) + s2);

if handles.noiseMult > 0
      randNoise = handles.numPElectrons.*handles.noiseMult.*randn(handles.numYpixels,handles.numXpixels);
      imgTotN = round(imgTot + randNoise);
      %we know that at a low number of p-electrons this causes de-whitening
      %of the noise.
      imgTotN(imgTotN < 0) = 0;
%           imgTotN(imgTotN > 2*handles.numPElectrons) = handles.numPElectrons;
end

imgTotFFTt = fftshift(fft2(ifftshift(imgTotN)));

% [sizeY,sizeX] = size(imgTotFFTt);

imgTotFFTEx = [...
                imgTotFFTt imgTotFFTt imgTotFFTt;...
                imgTotFFTt imgTotFFTt imgTotFFTt;...
                imgTotFFTt imgTotFFTt imgTotFFTt
              ];
for k = 1:numChannels
    indX = indXT{k};
    indY = indYT{k};
    imgTotFFTExFilt = F{k}.*imgTotFFTEx;
    indStartX = indX - indRadX;
    indStopX = indX + indRadX;
    indStartY = indY - indRadY;
    indStopY = indY + indRadY;

    tempFFTExFilt = imgTotFFTExFilt(indStartY:indStopY,indStartX:indStopX);

    chanImg{k} = tempFFTExFilt;
end
    
dcImg = zeros(size(Z{1}));

% The user can enter the ROIs here if wanted. Defaults to the entire
% image

% imgROI = [300 550;1135 1385];
% imgROI = [280 540;1130 1390]; %ROI_1
% imgROI = [600 860;800 1060]; %ROI_2 n1
% imgROI = [400 800;950 1350]; %ROI_2 n2
% imgROI = [151 350;151 350]; %Synth
% imgROI = [1 200;1 200]; %Synth2
imgROI = [1 handles.numYpixels;1 handles.numXpixels];
% imgROI = [200 1480;200 1480];
% imgROI = [1 1680;1 1680];
% imgROI = [426 536;326 426]; %Hirakawa and Condat detail for lighthouse
% roiStr = 'ROI_2';
% roiStr = 'ROI_1';
% roiStr = '';

idx = 1;

tempColImgOg = zeros([size(Z{1}),3]);
tempColImgR = tempColImgOg;
for k = 1:numRows
    tempImg = zeros(size(chanImg{1}));

    %Apply Q matrix
    for j = 1:numCols    
        tempImg = tempImg + newLfInv(k,j).*chanImg{j};    
    end

    % IFT
    finalImg{k} = fftshift(ifft2(ifftshift(tempImg)));    
    figure(handles.filterTestPlots(k));
    set(handles.filterTestPlots(k),'Name',['Test Img : Filter ' num2str(k)],'NumberTitle','off');
    clf;    
    sp(1) = subplot(2,2,1);
    dispImg = real(finalImg{k}(1:end-1,1:end-1))./(imgMult);
    if idx < 4
        tempColImgR(:,:,idx) = dispImg;        
    end
    imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
    set(gca, 'XTick', [],'YTick',[]);
    axis square;
    title('Reconstructed Img');        
    colormap(gray(256));        
    colorbar;
    sp(2) = subplot(2,2,2);    
    dispImg = -abs(real(finalImg{k}(1:end-1,1:end-1))./(imgMult)-Z{k});
    imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
    set(gca, 'XTick', [],'YTick',[]);
    axis square;    
    colormap(gray(256));    
    colorbar;
    % Note, MSE is for ROI
    MSE = sum(sum(dispImg.^2))./((imgROI(1,2)-imgROI(1,1))*(imgROI(2,2)-imgROI(2,1)));
    mseStr = ['_MSE_' num2str(MSE,'%10.2e')];
    title(['Difference Img, MSE : ' num2str(MSE,'%10.2e')])        
    sp(3) = subplot(2,2,3);    
%     dispImg = Z{k}./baseMax(k);
    dispImg = Z{k};
    if idx < 4
        tempColImgOg(:,:,idx) = dispImg;        
        idx = idx + 1;
    end
    imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
    set(gca, 'XTick', [],'YTick',[]);
    axis square;
    title('Original Img');
%     if baseMax(k) > baseMin(k)
%         caxis([baseMin(k)./baseMax(k) 1]);
%     end    
    colormap(gray(256));        
    colorbar;
    linkaxes(sp,'xy');
    dcImg = dcImg + Z{k}.*imgTransFactor(k);    
end    

%DC Channel
dcBaseMax = max(max(dcImg));
dcBaseMin = min(min(dcImg));

cIdx = floor(numCols/2)+1;
finalImg{k+1} = fftshift(ifft2(ifftshift(chanImg{cIdx})))./imgMult;
figure(handles.filterTestPlots(k + 1));
set(handles.filterTestPlots(k+1),'Name','Test Img : Base Band','NumberTitle','off');
clf;
sp(1) = subplot(2,2,1);   
dispImg = real(finalImg{k+1}(1:end-1,1:end-1));
imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
set(gca, 'XTick', [],'YTick',[]);
axis square;
caxis([dcBaseMin dcBaseMax]);
colormap(gray(256));
title('Reconstructed DC Img');
colorbar;
sp(2) = subplot(2,2,2);    
dispImg = -abs(dispImg-dcImg);
imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
set(gca, 'XTick', [],'YTick',[]);
axis square;
colorbar;
colormap(gray(256));

MSE = sum(sum(dispImg.^2))./((imgROI(1,2)-imgROI(1,1))*(imgROI(2,2)-imgROI(2,1)));
mseStr = ['_MSE_' num2str(MSE,'%10.2e')];
title(['Difference Img, MSE : ' num2str(MSE,'%10.2e')])
sp(3) = subplot(2,2,3);    
imagesc(dcImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
set(gca, 'XTick', [],'YTick',[]);
axis square;
caxis([dcBaseMin dcBaseMax]);
title('Original DC Img');
colormap(gray(256));
colorbar;
linkaxes(sp,'xy');
figure(handles.filterTestPlots(k + 2));
sp(1) = subplot(1,2,1);   
imshow(tempColImgOg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2),:));
sp(2) = subplot(1,2,2);   
imshow(tempColImgR(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2),:));
linkaxes(sp,'xy');
handles.imgTotFFTt = imgTotFFTt;
