%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       MODULATION TRANSFER FUNCTION FROM EDGE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This does the following functions:
% Find the position of the edge in the image
% Using the edge position, find the distance of each point from the edge
% Plot "distance vs signal" using all the points
% Use a local weighted linear fit approach to average / smooth all the data
% points to get a smooth edge spread function
% Differentiate the results, and apply a little more smoothing to get the
% point spread function
% Perform the FFT
%
% In the source code, various settings can be changed relating to how
% the data points are smoothed, what data ranges we work with etc.
%
% One big question is how to find the equation of the original line
% properly - even a small error in this can lead to problems. 
% This script interpolates the data to 0.1 pixel resolution and finds the
% point in each column where we're closest to the middle of the edge.
% This means we're making good use of the intensity information. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Variables you need to set up beforehand:

% InitialImage
% InitialMask - 256*256 matrix. Any pixels you want masked should be set to 
% 0, others to 1

% SET THE PIXEL RANGE TO WORK WITH, IF NEEDED
% This allows to to crop to the region around the edge

% Planar uses Y 50-200. 3D should use Y 80-200
XMin=10;
XMax=246;
%YMin=50;
%YMax=200;
YMin=80;
YMax=200;

% Note that we should crop the image before processing:

EdgeImage=InitialImage(YMin:YMax,XMin:XMax);
Mask=InitialMask(YMin:YMax,XMin:XMax);

[Ysize,Xsize] = size(EdgeImage);


%%%%%%%%% MASKING BEFORE EDGE FINDING

% For edge finding processes, we'll apply median filtering to the noisy
% pixels so they don't throw off the results.

FullFiltEdgeImage=medfilt2(EdgeImage);

MaskFiltEdgeImage=EdgeImage;
MaskFiltEdgeImage(~Mask)=FullFiltEdgeImage(~Mask);


%%%%%%%%% EDGE FINDING

% Now, we want to find the middle of the edge spread function in each
% column of pixels. As a start point, need to actually find the range of
% values

MeanPixelVal=mean(MaskFiltEdgeImage(:));

% Histogram the pixel values
% Assuming we have 1000s of counts, and at least 20000 pixels in each
% region, then using 200 bins will give us reasonable precision and a
% decent no of counts in the appropriate bins. "Bins" is array of X-values
[ImageHist,bins] = hist(MaskFiltEdgeImage(:),200);

figure;
stairs(bins,ImageHist)

% The histogram should have one peak above the mean and one below it. Find
% them, and get the midrange
HighHist=ImageHist.*(bins>MeanPixelVal);
LowHist=ImageHist.*(bins<MeanPixelVal);
[HighNo,HighBin] = max(HighHist);
HighSignal=bins(HighBin)
[LowNo,LowBin] = max(LowHist);
LowSignal=bins(LowBin)

MidSignal=(HighSignal+LowSignal)/2

% Blur the image used for edge finding to reduce noise effects, if needed,
% using a convolution. Syntax below
% blur=ones(3,3)./9;
% FindEdgeImage = conv2(A,B)

% Interpolate the image to get 0.1 pixel spacing in the Y direction.

XI=1:1:Xsize;
YI=(1:0.1:Ysize).';

FindEdgeImageInterp = interp2(MaskFiltEdgeImage,XI,YI);

% Find the Y-point in each column with the closest value
FindEdgeImage_DiffFromMid = abs(FindEdgeImageInterp-MidSignal);
[EdgeValues,YEdgeIndices] = min(FindEdgeImage_DiffFromMid);

% Given vector of indices, need to grab corresponding points
YEdgePoints=YI(YEdgeIndices);



%%%%%%%%%%%%%%%
% LINEFIT EDGE FINDING

% Since our image has a single line, plus a few outliers, we should be
% able to find it using a robust line fit
% First, we need to turn the "Find edge" result into arrays of X-Y data. We
% do this by getting the row and column indices of all these points
%[YEdgePoints,XEdgePoints] = find(FindEdgeSobel);

% Next, do fit to a "1st order polynomial"
% poly1   Y = p1*x+p2

EdgeLineFitOptions = fitoptions('Method','LinearLeastSquares','Robust','Bisquare');

cEdgeLineFit = fit((XI).',YEdgePoints,'poly1',EdgeLineFitOptions);

% Get the fit parameters - M, C

EdgeLineParams=coeffvalues(cEdgeLineFit);
EdgeLineM=EdgeLineParams(1);
EdgeLineC=EdgeLineParams(2);

% Display the fit
figure;
plot(cEdgeLineFit,XI,YEdgePoints);
title('Edge fit - 3D','FontName','Arial','FontSize',16);
xlabel('X','FontName','Arial','FontSize',14);
ylabel('Y','FontName','Arial','FontSize',14);
axes_handle=gca;
%set(axes_handle,'FontName','Arial')
%set(axes_handle,'FontSize','12')

% Then, get this in polar co-ordinates
% The line is rho = x*cos(theta) + y*sin(theta)
% But also y=mx+c
% So, can convert back
% Note that, unlike before, we get the result in radians

ThetaLineRads=atan(-1/EdgeLineM)
RhoLine=EdgeLineC*sin(ThetaLineRads)
% For comparison, get Theta in degrees
ThetaLineDegrees=ThetaLineRads*180/pi

% These parameters will be in terms of the cropped image, so we can use
% them to find the distances directly provided we continue to work from the
% same cropped image.



%%%%%%%%%%% DISTANCE OF POINTS FROM LINE

% Then, need to find distance of each pixel from the line. Check area of
% image we will work with 


% DistPixels gives distance from line in terms of pixels
% Basic expression is:
% DistPixels(y,x)=x*cos(ThetaLineRads)+y*sin(ThetaLineRads)-RhoLine;
% Due to the repeated use of the same cos and sin terms, pre-calculate

CTheta=cos(ThetaLineRads);
STheta=sin(ThetaLineRads);

for x = 1:Xsize
    for y = 1:Ysize
    DistPixels(y,x)=x*CTheta+y*STheta-RhoLine;
    end
end
% Then, convert to microns
Pixelsize=55;
DistMicrons=DistPixels*Pixelsize;



%%%%%%%%%%%%%% CREATING VECTORS WITH SIGNAL VS DISTANCE


% When working with the signal vs distance results, we simply remove all
% the bad pixels, rather than trying to compensate with the median
% filtering.

% At each point, only include results in output vectors for unmasked pixels
% Also, introduce a range of distances to work with, to reduce no of points

% CHOOSE THE DISTANCE RANGE FROM THE EDGE WE WORK WITH
% A typical edge is spread across 100um or so - choose 500um either side to
% be safe

DistMin=-500;
DistMax=500;

n=1;
for x = 1:Xsize
    for y = 1:Ysize
        if ((Mask(y,x)==1) & (DistMicrons(y,x)>DistMin) & (DistMicrons(y,x)<DistMax))
            Distance(n)=DistMicrons(y,x);
            Signal(n)=EdgeImage(y,x);
            n=n+1;
        end
    end
end

% Sort the vector by distance
[Distance_sorted, i] = sort(Distance);
Signal_sorted = Signal(i);



%%%%%%%%%%%%% SMOOTHING THE DATA

% In tests, to smooth out these clusters, a robust lowess smoothing with a span of 2e-3
% (i.e. 0.2% of the data, a span of about half a pixel) works well. Note that the
% lowess weights the data by distance.

% NOTE - The "span" of the smoothing is expressed as a fraction of the no
% of data points. We need to use it in a form which doesn't get clobbered
% by changes in the range. So, will specify the span in microns, then use
% the DistMin and DistMax to convert it.

% From various tests, the settings below supress high-frequency components (due
% to pixel-to-pixel variations) without affecting the results below 20
% lp/mm that we're interested in.

% CHOOSE THE SETTINGS FOR THE SMOOTHER
% Default 20um
% SpanMicrons=20
% If the data quality is good, less smoothing might work well
SpanMicrons=15;
SpanFraction=SpanMicrons/(DistMax-DistMin)

% For each point, the lowess smoother takes the surrounding data points
% over a range SpanMicrons, does a weighted linear fit to these points
% (with the weight decreasing with distance from the relevant point) and
% then uses this to find the smoothed point value.
% Robust lowess (rlowess) also ignores points more than 6 sigma from the fit
% It's also possible to use loess/rloess, which does a quadratic fit rather
% than a linear one.
Signal_rlowess = smooth(Distance_sorted,Signal_sorted,SpanFraction,'rlowess');

% Plot the edge spread function before and after smoothing
figure;
plot(Distance_sorted,Signal_sorted,'+',Distance_sorted,Signal_rlowess);
title('Edge Spread Function, with smoothing - 3D, THL358','FontName','Arial','FontSize',16)
xlabel('Distance (microns)','FontName','Arial','FontSize',14)
ylabel('Amplitude','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')



%%%%%%%%% INTERPOLATING AND DIFFERENTIATING THE SIGNAL

% Even after the smooth, the data points are still spaced very irregularly,
% which can cause problems with differentiation.
% So, we do linear interpolation of the function before finding the
% derivative. 
% Will use a spacing of 1um - gives Nyquist frequency of 100 lp/mm, which
% is a factor of 10 higher than we expect to have to work with.

% Can set the interpolation spacing as a variable here.

% SET INTERPOLATION SPACING
%Interp_spacing=1
% Once again, with good data smaller interpolation can give better results
Interp_spacing=0.5

Distance_interp = DistMin:Interp_spacing:DistMax;
Signal_interp = interp1(Distance_sorted,Signal_rlowess,Distance_interp,'linear');

% Now we have regular spacings, can just differentiate by finding
% difference between successive points.

DiffSignal_interp = diff(Signal_interp);

% Need to cut out final element of distance array before plotting
Distance_interp_short=Distance_interp;
Distance_interp_short(end) = [];

plot(Distance_interp_short,DiffSignal_interp)



%%%%%%%%%%% SMOOTHING THE LINE SPREAD FUNCTION

% After this, we then apply some robust smoothing, to get rid of outliers
% and reduce random fluctuations. The smoothing settings are pretty low
% here.

% CHOOSE SMOOTHER SETTINGS
DiffSpanMicrons=5;
DiffSpanFraction=DiffSpanMicrons/(DistMax-DistMin);

DiffSignal_rlowess = smooth(Distance_interp_short,DiffSignal_interp,DiffSpanFraction,'rlowess');

% Since these variables are very useful later on, will rename them

LineSpreadFunction=DiffSignal_rlowess;
LineSpreadDistance=Distance_interp_short;

LineSpreadFunction_3D_THL358=LineSpreadFunction;
LineSpreadDistance_3D_THL358=LineSpreadDistance;

% Plot the Line Spread Function
figure
plot(Distance_interp_short,DiffSignal_interp,'+',Distance_interp_short,DiffSignal_rlowess);
title('Line Spread Function - 3D, THL358','FontName','Arial','FontSize',16')
xlabel('Distance (microns)','FontName','Arial','FontSize',14')
ylabel('Amplitude','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([-100 100 0 1]);
axis 'auto y';




%%%%%%%%%%%%%% FAST FOURIER TRANSFORM

% If necessary, we can choose to use a limited range of the data (in um)
% If not, can comment out the next 2 lines of code, and just use FFTInputSignal=DiffSignal_rlowess;

FFTInputEndpoint=60;
FFTInputSignal=DiffSignal_rlowess(1:((FFTInputEndpoint-DistMin)/Interp_spacing));
% Recall that we interpolated the data in 1um steps as below:
% Distance_interp = DistMin:Interp_spacing:DistMax;


% Then, since the FFT algorithm works on vectors with 2^n points, find next
% highest power of 2

FFTLength = 2^nextpow2(length(FFTInputSignal)); 

% Find FFT
FFTOutput = fft(FFTInputSignal,FFTLength);

%The no of elements in FFTOutput will match FFTLength. The highest
%spatial frequency component present will be equal to 1/Spatial freq of
%data. However, the Nyquist frequency will be half of this. So, only want
%to look at results up to this. Note that the first term listed here is the
%Nyquist frequency, given by 0.5/(1um)=0.5/(0.001mm)=500 lp/mm when our
%original function was sampled at 1um intervals.

FNyquist=0.5/(Interp_spacing*0.001);
SpatialF = FNyquist*linspace(0,1,FFTLength/2);

% Normalise the FFT using the DC component

FFTNorm=FFTOutput/FFTOutput(1);

% Plot single-sided amplitude spectrum, up to the Nyquist frequency
figure;
plot(SpatialF,abs(FFTNorm(1:length(FFTNorm)/2))) 
title('3D - THL385 - Modulation Transfer Function','FontName','Arial','FontSize',16);
xlabel('Frequency (lp/mm)','FontName','Arial','FontSize',14)
ylabel('|MTF|','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 25 0 1]);


% Comments only - this is the theoretical curve, calculated on the basis
% that we have a 55um pixel (18.2 lp/mm) with full sensitivity to X-rays
% falling within the pixel, zero to those falling outside, and a perfect 
% edge pattern falling on the detector. 
% TheoreticalCurve=sinc(SpatialFPlanar*pi/18.2);

% Example code for plotting multiple results together
% plot(SpatialFPlanar,abs(FFTNormPlanar(1:length(FFTNormPlanar)/2)),SpatialF3D1,abs(FFTNorm3D1(1:length(FFTNorm3D1)/2)),SpatialFPlanar,TheoreticalCurve)

% Some code for saving variables with new names
%LSF_distance_3D_THL358=Distance_interp_short;
%LineSpreadFunction_3D_THL358=DiffSignal_rlowess;
%save




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FINDING INFORMATION FROM THE LSF
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The variation in the LSF with the threshold gives us info about charge
% sharing, and the LFS as a whole gives us info about the behaviour of the
% pixel

% The following section will automatically find information from the LSF.



% Sort in ascending order
LSFsort= sort(LineSpreadFunction);
PeakWidth=40; % Chosen by user
NPeakSamples=PeakWidth/Interp_spacing;

PeakAmplitude=mean(LSFsort((end-NPeakSamples):end))

% Create separate variables for the LSF below and above distance zero

LSF_neg=LineSpreadFunction(LineSpreadDistance<=0);
LSDist_neg=LineSpreadDistance(LineSpreadDistance<=0);
LSF_pos=LineSpreadFunction(LineSpreadDistance>=0);
LSDist_pos=LineSpreadDistance(LineSpreadDistance>=0);

% Find the positions closest to half maximum
LSF_neg_DiffFromMid=LSF_neg-(PeakAmplitude./2);
[LSF_neg_MidEdgeVal,LSF_neg_MidEdgeIndex]=min(abs(LSF_neg_DiffFromMid));

LSF_neg_MidEdgeVal;
LSF_neg_MidEdgePos=LSDist_neg(LSF_neg_MidEdgeIndex)

LSF_pos_DiffFromMid=LSF_pos-(PeakAmplitude./2);
[LSF_pos_MidEdgeVal,LSF_pos_MidEdgeIndex]=min(abs(LSF_pos_DiffFromMid));

LSF_pos_MidEdgeVal;
LSF_pos_MidEdgePos=LSDist_pos(LSF_pos_MidEdgeIndex)

FullWidthHalfMaximum=LSF_pos_MidEdgePos-LSF_neg_MidEdgePos


% Next, want a measure of the steepness of the rising and falling edges
% Do this by finding the 20% to 80% rise and fall distance.
[LSF_neg_val,LSF_neg_RisingEdgeIndex1]=min(abs(LSF_neg-0.2*PeakAmplitude));
LSF_RisingEdge1=LSDist_neg(LSF_neg_RisingEdgeIndex1)

[LSF_neg_val,LSF_neg_RisingEdgeIndex2]=min(abs(LSF_neg-0.8*PeakAmplitude));
LSF_RisingEdge2=LSDist_neg(LSF_neg_RisingEdgeIndex2)

LSF_RisingEdgeDist=LSF_RisingEdge2-LSF_RisingEdge1

%LSF_pos=LSF_pos(1:160);
[LSF_pos_val,LSF_pos_FallingEdgeIndex1]=min(abs(LSF_pos-0.2*PeakAmplitude));
LSF_FallingEdge1=LSDist_pos(LSF_pos_FallingEdgeIndex1)

[LSF_pos_val,LSF_pos_FallingEdgeIndex2]=min(abs(LSF_pos-0.8*PeakAmplitude));
LSF_FallingEdge2=LSDist_pos(LSF_pos_FallingEdgeIndex2)

LSF_FallingEdgeDist=LSF_FallingEdge1-LSF_FallingEdge2

% Measure of variation in signal due to beam intensity etc
LSFVariationSigma=std(LSF_pos((50/Interp_spacing):(100/Interp_spacing)))