%function [YSpatialF,NPSNormalised_1D,NPSNormalised]=NoisePowerSpectrum(FilenamePrefix,StartIndex,StopIndex, Step)

% Set pixel size info
Pixel_size_mm = 0.055;

% Test data
FilenamePrefix='S95Planar_THL387_tube_flatfield__r'
StartIndex=0;
StopIndex=99;
Step=1;

% When finding the NPS, we need to find the "individual" NPS for a large no
% of images, then take the average. So, the main routine is a loop over all
% the images

% Have to decide how to treat dead and noisy pixels
% I will just insert the mean value for now, since the no of noisy pixels
% is small.

% Should have already calculated the mean of all the images, the flat-field
% coefficients, and the mask.



% Correct the mean image with the flat-field correction
MeanImageCorrected=MeanImage.*FFCoeff;
% Get the mean of the unmasked pixels
OverallMeanPixel=mean(MeanImageCorrected(Mask))
MeanImageCounts=sum(MeanImageCorrected(:));

% Initialisation of NPS parts
NImages=0;
NPSSum=zeros(256,256);
% Loop over all the images we've got
for k=StartIndex:Step:StopIndex
    % For each k value,
    % grab the corresponding file
    if ((k<10)&(StopIndex>=10))
        Filename=[FilenamePrefix '0' num2str(k) '_.txt'];
    else
        Filename=[FilenamePrefix num2str(k) '_.txt'];
    end

    Image=load(Filename);
    
    % Apply the flat-field correction
    ImageFF=Image.*FFCoeff;
    
    % When doing this, we want to reject the images with significantly
    % different count rates, since I messed up a couple of acquisitions.
    
    ImageCounts=sum(ImageFF(:));
    
    if ((abs(ImageCounts-MeanImageCounts))<(MeanImageCounts/10))
        % Do the processing
      
        % Find the deviation from the mean. This is set to zero for all masked
        % pixels, to stop them screwing up the results
        ImageDev=ImageFF-OverallMeanPixel;
        ImageDev(~Mask)=0;

        % At this point, we might still have bad or noisy pixels. We can
        % compensate for this by finding pixels more than 6 std. devs from
        % zero, and setting these to zero. Note this is optional - for
        % perfectly good sensor, shouldn't be necessary.

        ImageDevCorrected=ImageDev;
        %{
        ImageDevStdDev=std(ImageDev(:));
        ImageDevBadPixels=(abs(ImageDev)>(ImageDevStdDev*5));
        ImageDevCorrected(ImageDevBadPixels)=0;
        %}
        % Take 2D FFT of this

        FFTImage=fft2(ImageDevCorrected);

        % Find square of abs value, and then divide by the no of pixels, to find
        % the effective NPS for the particular image
        NPSImage=(abs(FFTImage)).^2./(256^2);

        NPSSum=NPSSum+NPSImage;
        % While doing this, count the images loaded
        NImages=NImages+1;    
        
    else
        % Reject the image
        Message=['The file ' Filename ' was rejected'];
        disp(Message);
    end
    

end

% After looping through all this and summing the results, want to divide by
% no of images used and then normalise by dividing by mean intensity.

NPSRaw=NPSSum./(NImages);

% Then, have to decide how to work with it. Current version is NPS, in
% terms of counts per pixel in each image, without normalisation. Better to
% normalise it by dividing by the mean no of counts per pixel, getting a
% dimensionless result, for use in the DQE.

NPSNormalised=NPSRaw./OverallMeanPixel;

NPS_per_mm2=NPSRaw./(Pixel_size_mm.^2);

% Generate x- and y-vectors with the corresponding spatial frequencies

XSpatialF=(0:1:255)./(Pixel_size_mm*256);
YSpatialF=(0:1:255)./(Pixel_size_mm*256);

% Then, get the 1-D NPS, along the Y-direction (perp to the edge). Get this
% by taking the column y(:), x=1.

NPSRaw_1D=NPSRaw(:,1);
NPSNormalised_1D=NPSNormalised(:,1);

% Display this
figure;
plot(YSpatialF,NPSRaw_1D);

figure;
plot(YSpatialF,NPSNormalised_1D);

% Try using a projection with more points
NPSNormalised_1D_alt=mean(NPSNormalised(1:8,:));

% Then, get a bit of smoothing
SmoothSpan_lppmm=1.5;
SmoothSpanFraction=SmoothSpan_lppmm*Pixel_size_mm;

NPSNormalised_1D_smoothed = smooth(YSpatialF,NPSNormalised_1D_alt,SmoothSpanFraction,'rlowess');


figure;
plot(YSpatialF,NPSNormalised_1D_alt,YSpatialF,NPSNormalised_1D_smoothed);
title('Planar - THL387 - Normalised Noise Power Spectrum','FontName','Arial','FontSize',16);
xlabel('Frequency (lp/mm)','FontName','Arial','FontSize',14)
ylabel('Noise Power / hits detected','FontName','Arial','FontSize',14)
fig_handle=gcf;
set(fig_handle,'Units','centimeters')
set(fig_handle,'Position',[3,3,20,15])
% Position vector is [left, bottom, width, height]
axes_handle=gca;
set(axes_handle,'FontSize',12)
set(axes_handle,'FontName','Arial')
axis([0 18 0 2]);