clear; clc;
I = imread('slanted.jpg'); % gray image!
G = rgb2gray(I);
pixelSize = 2.5e-3; % mm
fNyquist = 1/(2*pixelSize);

dim = size(I);
up  = dim(1);
integmax = -1.0e22;
for cut = 90:90
D   = double(I(:,cut));      % Get a portion of image and convert to double
LSF = cda(D,pixelSize);      % Line Spread Function
MTF = abs(fft(LSF)).^2;      % MTF = FFT of LSF           
MTF = MTF/max(MTF);          % Nomalized MTF
MTF = MTF';
N   = length(LSF);           % Number of points 
f   = (1/pixelSize)*(0:N-1)/N; % Frequency = cycles/mm
x   = (pixelSize)  *(0:N-1); % Position mm

MTFdetector = sin(pi*f*pixelSize) ./ (pi*f*pixelSize);
MTFdetector(1) = 1;

figure(1);
integ = trapz(f(1:10),MTF(1:10));
if integ > integmax
    integmax = integ
    cut
end
str = fprintf('Original image, cut = %d, integral = %f / %f\n',cut,integ, integmax);
imshow(I);
title(str);
rectangle('Position',[cut 1 1 up],'LineWidth',0.5,'EdgeColor','g')

figure(2);
subplot(2,2,1); imshow(uint8(D'));
                title('Cropped image')
subplot(2,2,3); plot(x,D(2:end-1));
                xlabel('Position (mm)');                
                ylabel('Intensity (arb.u.)');
                grid on
subplot(2,2,2); plot(x,LSF);
                xlabel('Position (mm)');
                ylabel('LSF (intensity/mm)');
                grid on
subplot(2,2,4); upper = 20;
                plot(f(1:upper),MTF(1:upper),'k--'); hold on
                plot(f(1:upper),MTFdetector(1:upper),'g');
                plot(f(1:upper),MTF(1:upper)./MTFdetector(1:upper),'b');
                xlabel('Frequency (lp/mm)');
                ylabel('MTF (%)');
                grid on
                hold off
pause(0.01);
end

% first derivative using central differrence approximation method.
function der = cda(y,h)
  j = 1;
  n = length(y);
  der = zeros(n-2);
  for i = 2:n-1
     der(j) = (y(i+1)-y(i-1))/(2*h);
     j = j + 1;
  end
end

