function [ meanMSE ] = genMeanMSEResultAWGN( FFTblocks, filtFunc, filtParams1,filtParams2,numNoiseRealizations,...
                                                     X,Y,noiseMagVec, imgParams, newLfInv  )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here


%tic();

numChannels = length(FFTblocks);

for k = 1:numChannels
    F{k} = planckTaper2D(filtParams1(k),filtParamsk(k),X,Y);    
end

for l=1:numNoiseRealizations
    if usePoissonNoise
        imgTotN = fPoissrnd( imgTot, 1000 );    
    end

    % snrStr = ['SNR_' num2str(snrIndB,'%5.1f') 'dB'];
    %thermal and readout noise round(sqrt(s2).*randn(size(N)) + s2);
    % imgTot = round(imgTot + fPoissrnd(imgMean*handles.noiseMult.*ones(handles.numYpixels,handles.numXpixels),500));
    if noiseMagVec(ll) > 0
    %     if handles.noiseMult > 0    
%         imgTotN = round(imgTot + fPoissrnd(handles.numPElectrons*handles.noiseMult.*ones(handles.numYpixels,handles.numXpixels),500));    
%         imgTotN = round(imgTot + fPoissrnd(handles.numPElectrons*noiseMagVec(ll).*ones(handles.numYpixels,handles.numXpixels),500));
          randNoise = handles.numPElectrons.*noiseMagVec(ll).*randn(handles.numYpixels,handles.numXpixels);
          imgTotN = round(imgTot + randNoise);
          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};
    %     F = planckTaper2D(handles.radius(k),handles.eps(k),X-chanCenters(k,1),Y-chanCenters(k,2));        
    %     F = planckTaper2D(handles.radius(k),epS,X-chanCenters(k,1),Y-chanCenters(k,2));        
    %     figure(handles.filterTestPlots(1));
    %     imagesc(F);
    %     ax = gca;
    %     ax.YDir = 'normal';
    %     pause(0.5);
        imgTotFFTExFilt = F{k}.*imgTotFFTEx;
    %     imagesc(imag(imgTotFFTExFilt));
    %     ax = gca;
    %     ax.YDir = 'normal';
    %     pause(0.5);
        indStartX = indX - indRadX;
        indStopX = indX + indRadX;
        indStartY = indY - indRadY;
        indStopY = indY + indRadY;

        tempFFTExFilt = imgTotFFTExFilt(indStartY:indStopY,indStartX:indStopX);
    %     tempFFTExFilt = fftshift(fft2(2*real(ifft2(ifftshift(tempFFTExFilt)))));

        chanImg{k} = tempFFTExFilt;
    %     imagesc(real(tempFFTExFilt));
    %     chanImg{k} = fftshift(2*real(ifft2(ifftshift(tempFFTExFilt))));
    %     imagesc(real(fftChanImg{k}));

    end
    

    % for k = 1:chanHalfPlaneIdx
    %     if isEdge(k) == 1
    %         chanImg{k} = chanImg{k} + chanImg{numChannels-k+1};
    % %         chanImg{k} = chanImg{k} + (chanImg{numChannels-k+1});
    %         chanImg{numChannels-k+1} = conj(chanImg{k});
    %         tempFFTExFilt = chanImg{k};
    %         figure(20);
    %         imagesc(real(tempFFTExFilt));
    %         figure(21);
    %         imagesc(imag(tempFFTExFilt));
    %     elseif isEdge(k) == 2
    %         %todo : add corner case, need to add all 4 corners appropriately
    %     else
    %     end
    % end

    % dcImg = zeros(size(Z{1}));
    % imgROI = [300 550;1135 1385];
    % imgROI = [280 540;1130 1390]; %ROI_1
%       imgROI = [600 860;800 1060]; %ROI_2
    imgROI = [151 350;151 350]; %Synth
    % imgROI = [1 1680;1 1680];
%     imgROI = [200 1480;200 1480];
    % imgROI = [1 1400;1 1400];
    roiStr = 'ROI_2';
    % roiStr = '';
    diffCRange = [-0.2 0];

    genTikz = false;
    
    for k = 1:numRows
        tempImg = zeros(size(chanImg{1}));
        for j = 1:numCols
    %         figure(15);
    %         imagesc(real(fftChanImg{j}));
    %         ax = gca;
    %         ax.YDir = 'normal';
            tempImg = tempImg + newLfInv(k,j).*chanImg{j};
    %         figure(16);
    %         imagesc(real(tempImgFFT));
    %         ax = gca;
    %         ax.YDir = 'normal';
        end
    %     figure(16);
    %     imagesc(real(tempImgFFT));

        finalImg{k} = fftshift(ifft2(ifftshift(tempImg)));
    %     finalImg{k} = tempImg;
    %     filtStr = ['_filtRad_' num2str(handles.radius(k),'%1.3f')];
    %     handles.filtStr = filtStr;
        figure(handles.filterTestPlots(k));
        clf;    
    %     sp(1) = subplot(2,2,1);
        dispImg = real(finalImg{k}(1:end-1,1:end-1))./(baseMax(k)*imgTransFactor(k)*imgMult);
%         dispImg = real(finalImg{k}(1:end-1,1:end-1))./imgMax;
        imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
        set(gca, 'XTick', [],'YTick',[]);
        axis square;
    %     title('Reconstructed Img');
    %     caxis([baseMin(k) baseMax(k)]);
%         caxis([imgMin./imgMax 1]);
        if baseMax(k) > baseMin(k)
            caxis([baseMin(k)./baseMax(k) 1]);
        end
    %     colormap(magma(256));
        colormap(gray(256));
        if genTikz
            if handles.useSynthTestImg
    %             matlab2tikz(['./figures/imgReconstructionSynth_' num2str(handles.xSizeD)...
    %                          'x' num2str(handles.ySizeD) '_' roiStr filtStr poissStr snrStr '_filtNum_'...
    %                          num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
                matlab2tikz(['./figures/imgReconstruction_' num2str(handles.xSizeD)...
                             'x' num2str(handles.ySizeD) '_Synth_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            else
                matlab2tikz(['./figures/imgReconstruction_' num2str(handles.xSizeD)...
                             'x' num2str(handles.ySizeD) '_' roiStr '_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            end
        end
    %     colorbar;
    %     sp(2) = subplot(2,2,2);    
    %     dispImg = -abs(real(finalImg{k}(1:end-1,1:end-1))-Z{k})./baseMax(k);
    %     dispImg = -abs(real(finalImg{k}(1:end-1,1:end-1))-Z{k})./imgMax;
%         dispImg = -abs(real(finalImg{k}(1:end-1,1:end-1))./(max(max(real(finalImg{k}))))-Z{k});
          dispImg = -abs(real(finalImg{k}(1:end-1,1:end-1))./(baseMax(k)*imgTransFactor(k)*imgMult)-Z{k});
%           dispImg = dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2));
        imagesc(dispImg(imgROI(1,1):imgROI(1,2),imgROI(2,1):imgROI(2,2)));
        set(gca, 'XTick', [],'YTick',[]);
        axis square;
    %     imagesc(real(finalImg{k}(1:end-1,1:end-1)));
    %     colormap(magma(256));
        colormap(gray(256));
    %     caxis([minZ maxZ])
    %     colorbar;
    %     caxis([baseMin(k)./baseMax(k) 1]);
        caxis(diffCRange);
        MSEtemp(l,k,ll) = sum(sum(dispImg.^2))./((imgROI(1,2)-imgROI(1,1))*(imgROI(2,2)-imgROI(2,1)));
    %     maxDiff = max(max(dispImg))
    %     minDiff = min(min(dispImg))
    %     MSE = sqrt(sum(sum(dispImg.^2)))./(handles.numPElectrons);
    %     mseStr = ['_MSE_' num2str(MSE,'%10.2e')];
    %     title(['Difference Img, MSE : ' num2str(MSE,'%10.2e') ' : ' num2str(minDiff,'%10.2f') ' : ' num2str(maxDiff,'%10.2f')])
        if genTikz
            if handles.useSynthTestImg
    %             matlab2tikz(['./figures/imgDiffSynth_' num2str(handles.xSizeD)...
    %                          'x' num2str(handles.ySizeD) '_' roiStr filtStr poissStr snrStr mseStr '_filtNum_'...
    %                          num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
                matlab2tikz(['./figures/imgDiff_' num2str(handles.xSizeD)...
                             'x' num2str(handles.ySizeD) '_Synth_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            else
    %             matlab2tikz(['./figures/imgDiff_' num2str(handles.xSizeD)...
    %                          'x' num2str(handles.ySizeD) '_' roiStr filtStr poissStr snrStr mseStr '_filtNum_'...
    %                          num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
                matlab2tikz(['./figures/imgDiff_' num2str(handles.xSizeD)...
                             'x' num2str(handles.ySizeD) '_' roiStr '_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            end
        end
%     %     sp(3) = subplot(2,2,3);    
        dispImg = Z{k}./baseMax(k);
        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(magma(256));
        colormap(gray(256));
        if genTikz
            if handles.useSynthTestImg
                matlab2tikz(['./figures/imgOriginal_Synth_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            else
                matlab2tikz(['./figures/imgOriginal_' roiStr '_c'...
                             num2str(k) '.tikz'],'width','\fpITwidth','height','\fpITheight');
            end
        end
    %     colorbar;
    %     linkaxes(sp,'xy');
    %     dcImg = dcImg + Z{k};
    %     title(['Difference Img, MSE : ' num2str(MSE,'%10.2e') ' : ' num2str(minDiff,'%10.2f') ' : ' num2str(maxDiff,'%10.2f')]);
    end    
end

toc()

end

