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