function [handles ] = optimizeFilters(handles)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

    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);

    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);

    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);
    imgTot = zeros(handles.numYpixels,handles.numXpixels);
    
    for k = 1:handles.numFiltersD
        if handles.useSynthTestImg
            Z{k} = imgFunc2D(handles.numXpixels,handles.numYpixels,pi/4);    
        else
            if isempty(handles.hyperReflectanceData)
                load('hyperReflectanceData');
                handles.hyperReflectanceData = hyperReflectanceData;
            end
            Z{k} = imgFromHyperSpecReflectance(handles.numXpixels,handles.numYpixels,handles.hyperReflectanceData,k);
            imgParams.Zfft(:,:,k) = fftshift(fft2(ifftshift(Z{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;
    
    imgMeanAfterMask = mean(mean(imgTot));
    
    imgMean = mean(baseMean);
    imgMult = handles.numPElectrons/imgMeanAfterMask;

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

    imgTotFFTt = fftshift(fft2(ifftshift(imgTot)));
    imgTotFFTEx = [...
                    imgTotFFTt imgTotFFTt imgTotFFTt;...
                    imgTotFFTt imgTotFFTt imgTotFFTt;...
                    imgTotFFTt imgTotFFTt imgTotFFTt
                  ];              
       
    for k = 1:numChannels
        indX = indXT{k};
        indY = indYT{k};    
        indStartX = indX - indRadX;
        indStopX = indX + indRadX;
        indStartY = indY - indRadY;
        indStopY = indY + indRadY;

        tempFFTExFilt = imgTotFFTEx(indStartY:indStopY,indStartX:indStopX);
        FFTblocks(:,:,k) = tempFFTExFilt(1:end-1,1:end-1);
        
%         radius(k) = 1/6;
%         eps(k) = 0.1;
    end
    
    [X,Y] = meshgrid(handles.xiC(1:end-1),handles.etaC(1:end-1));
    
    imgParams.numPElectrons = handles.numPElectrons;
    imgParams.imgTransFactor = imgTransFactor;
    imgParams.baseMax = baseMax;
    imgParams.baseMean = baseMean;
    imgParams.imgMult = imgMult;
    
    NN = 1;
    numGAGens = 16;
%     meanMSE = genMeanMSEResultAWGN( FFTblocks, @planckTaper2D, radius,eps,...
%                                     NN, X, Y, handles.noiseMult,...
%                                     imgParams, newLfInv  );

    genMeanMSEResultAWGNEval = @(filtParams)genMeanMSEResultAWGN( filtParams, FFTblocks, @planckTaper2D, NN, X, Y, handles.noiseMult, imgParams, newLfInv  );
                                
%     gaoptions = optimoptions('ga','MaxGenerations',1,'Display','iter','UseParallel',true);

    popSize = 200;
    initPop(:,1) = linspace(0.05,0.5,popSize);
    initPop(:,2) = 0.05*ones(popSize,1);
    UB = [0.5; 1];
    LB = [0; 0];
    
    gaoptions = gaoptimset('Generations',numGAGens,'Display','iter','UseParallel',true,...
        'MutationFcn',@mutationadaptfeasible,'InitialPopulation',initPop,'PopulationSize',popSize,'OutputFcn',@gaStateFcn);    
    
    [x,fval] = ga(genMeanMSEResultAWGNEval,2,[],[],[],[],LB,UB,[],gaoptions)
    
    for k = 1:numChannels
        handles.radius(k) = x(1);
        handles.eps(k) = x(2);
    end
end

