function mtfscan(varargin)
if (nargin==0 | nargin>6)
   infile = input('Enter the input file name','s');
elseif (nargin>=1)
   infile = varargin{1};
end
infilep = strrep(infile, '\','\\');  % To print correctly in plots.
infilep = strrep(infilep,'_','\_');  % To print correctly in plots.
close all

fileext = infile(findstr(infile(3:length(infile)),'.')+3:length(infile));  % File type
inf = imfinfo(infile,fileext);  % Fileext should be a graphic file extension.
coltyp = inf.ColorType;  fsiz = inf.FileSize;
fprintf('Loading and displaying a %g byte %s image; may take a while...\n',fsiz,coltyp);

testimg = imread(infile);  szim = size(testimg);  sx = szim(2);  sy = szim(1);
if (strcmp(inf.ColorType,'truecolor'))  % RGB color file: must combine!!!
   imgtemp = mean(testimg,3);
   figure;  image(imgtemp);  colormap(gray(256));
   clear testimg;  testimg = imgtemp;  clear imgtemp;
else
   testimg = double(testimg);
   figure;  image(testimg);  colormap(gray(256));
end

fprintf('Image size is %g x %g.\n',szim(1),szim(2));
y1 = input('Choose first y-level from the graph for analysis: ');
y2 = input('Choose 2nd y-level from the graph for analysis (or none): ');
pmagn = input('Enter the print magnification: ');
if (isempty(pmagn))  pmagn = 1;  end
scandpi = sx*2.54/pmagn;  disp(['Scan dpi = ' num2str(scandpi,5)]);

% Code from MTFcurve for pattern.
ymax = 200; ymin = 2;   % Max, min lines per mm (spatial frequency) the plots.
ncyc = log10(ymax/ymin);      % Cycles of frequency to display. Argument is frequency ratio.
xpmax = 10;    % Max x1 mm for printer test image.
sdpi = floor(sx*25.4/xpmax);  % Scan dpi,
lpfft = 16384; lpdsp = 15360; % Counterparts of lfft... for printer.
lpd1 = szim(2)-1;  %  lpd1 = lpdsp-1;
xp1 = ([0:lpd1]/lpd1)*xpmax;  % Distance in mm along axis to display.
xpinc = xpmax/lpd1;           % Spatial increment in mm.
fpq1 = 10.^(ncyc*[0:lpd1]/lpd1);  fpq1 = ymax*fpq1/max(fpq1);  % Spatial frequency 1/mm.
u1 = 2.*pi*cumsum(fpq1)*xpinc;  % Equation for varying frequency is tricky!
ypw = .5*(1+cos(u1));  % Original cosine pattern before rounding.

xlbl = ['lp/mm on ' num2str(xpmax,3) ' mm target; (*' num2str(25.4/pmagn,4) ...
   ' for print lp/in at ' num2str(pmagn,4) 'X)'];
if (pmagn<1.01) xlbl = [xlbl '; ' num2str(sdpi) ' dpi'];  end


if (~isempty(y2))
   xofs = -1;  figure;
   semilogx(fpq1,testimg(y1,:)/255,fpq1,testimg(y2,:)/255+xofs);
   ylabel(['   (Density-1) for y = ' num2str(y2) '         Density for y = ' num2str(y1)]);
else
   figure; semilogx(fpq1,testimg(y1,:)/255); xofs = 0;
   ylabel(['    Density for y = ' num2str(y1)]);
   p1 = get(gcf,'Position');  p1(4) = .52*p1(4);  %  p1(3) = .8*p1(3);
   set(gcf,'Position',p1);
end
hold on;  semilogx(fpq1,.08*ypw+xofs-.09,'r-');
ax = axis;  % ax(1) = ymin;  ax(2) = ymax;
ax(1) = ymin;  ax(2) = 100;  ax(3) = xofs-.1;  ax(4) = 1;  axis(ax)
title(['Print resolution for ' num2str(pmagn,4) 'X magnification:  ' infilep]);
xlabel([xlbl ';  scan dpi = ' num2str(scandpi,5)]);  grid on;  zoom on;  hold off;

% Repeat the figure, smaller, for web page.

if (~isempty(y2))
   xofs = -1;  figure;
   semilogx(fpq1,testimg(y1,:)/255,fpq1,testimg(y2,:)/255+xofs);
   ylabel(['   (Density-1) for y = ' num2str(y2) '         Density for y = ' num2str(y1)]);
else
   figure; semilogx(fpq1,testimg(y1,:)/255); xofs = 0;
   ylabel(['    Density for y = ' num2str(y1)]);
   p1 = get(gcf,'Position');  p1(4) = .52*p1(4);  %  p1(3) = .8*p1(3);
   set(gcf,'Position',p1);
end
p1 = get(gcf,'Position');  p1(4) = .88*p1(4);  p1(3) = .72*p1(3);
set(gcf,'Position',p1);
hold on;  semilogx(fpq1,.08*ypw+xofs-.09,'r-');
ax = axis;  % ax(1) = ymin;  ax(2) = ymax;
ax(1) = ymin;  ax(2) = 100;  ax(3) = xofs-.1;  ax(4) = 1;  axis(ax)
title(['Print resolution for ' num2str(pmagn,4) 'X magnification:  ' infilep]);
xlabel(xlbl);  grid on;  zoom on;  hold off;