function mtfcurveA(varargin)  % version 04/30/2002
% Visualize the effects of MTF (modulation transfer function)
% in film, lenses, scanners, and sharpening using simulation.  Format:
% mtfcurveA f50   % f50 is the 50% MTF density of film in line pairs/mm, typically 20 - 60.
% mtfcurveA f50 fboost   % fboost is the frequency shift of the 2nd order equalizer.
% mtfcurveA f50 fboost flens  % flens is the 50% MTF density of the lens.
% mtfcurveA f50 fboost flens lord  % lord is lens rolloff order (default = 2).
% mtfcurveA f50 fboost flens lord dscan sincpwr  % dscan is scanner pixels per mm (dpmm).
%   sincpwr is the power of the sinc function describing the scanner or digital sensor roloff.
% mtfcurveA f50 fboost flens lord dscan sincpwr ksharp  % ksharp is sharpening boost.
%
% mtfcurveA 45 13 61 2 4000/25.4 3  % test example
%
% Parameters default to 0 if not entered, except lord=2.
% If f50<0, Create printer test targets. |f50| is the 50% film density.
% If fboost>0, a boost is applied by shifting the peak frequency of the rolloff.
%    normal: 1/(1+(f/f50)^2)
%    with fboost: k1/(1+((f-fboost)/f50a)^2) ,  where k1 normalizes response to 1
%    at f=0 and f50a is set so f50 remains the 50% MTF density. fboost < f50/2.
% If flens>0, a simplified lens MTF is added to the transfer, using 1/(1+(f/flens)^2)
% If abs(dscan)>0, a scanner with dscan dpmm (pixels/mm) is added.
% sincpwr is the exponent of the sinc functhin describing scanner MTF: |sinc|^sincpwr
%    3 is appropriate for a film scanner; 2 for a Bayer sensor digital camera;
%    1 for a Foveon X3 digital camera; 4 or 5 for an inexpensive flatbed scanner.
% If abs(ksharp)>0, sharpening equivalent to a cosine boost is applied, where ksharp is 
%    the amount of boost. Sharpening radius rsharp = 1 if ksharp>0; 2 if ksharp<0.
%    Equation: (1+|ksharp|*cos(2*pi*f*rsharp/dscan))/(1+|ksharp|)   %%%%%%%%%%%%%%
%
% Copyright (C) 2001-2004 by Norman Koren, http://www.normankoren.com  normkoren@earthlink.net
% Every conceivable disclaimer applies; the author assumes no legal liability whatsoever.
% You may distribute this program freely, but please maintain the original copyright.

fboost = 0;  flens = 0; lord = 2;  dscan=0;  ksharp=0;  % Defaults.
if (nargin==0 | nargin>7)
   error(['Entered ' num2str(nargin) ' arguments. Must enter 1 - 6.']);  end
if (nargin>=1)
   f50 =    abs(str2num(varargin{1}));  % Density (line pairs per mm) of 50% MTF.
   descript = ['f50 = ' num2str(f50,4) ' lp/mm'];  end 
   prtarg = ((str2num(varargin{1}))<0);  % Create printer target if first argument <0.
if (nargin>=2)        % Boost frequency from cosine eql or peak center shift.
   fboost = str2num(varargin{2});  % Boost by shifting peak
   if (fboost>eps) descript = [descript '; fboost = ' num2str(fboost,4) ' lp/mm']; end, end
if (nargin>=3)        % Lens f50.
   flens =  str2num(varargin{3});  % Lens 50% density.
   descript = [descript '; flens = ' num2str(flens,4) ' lp/mm'];  end
if (nargin>=4)        % Lens rolloff order (default 2).
   lord =  str2num(varargin{4});  % Lens 50% density.
   if (lord<eps) lord = 2;  % Default value if 0 is entered.
   else descript = [descript '; lord = ' num2str(lord,3)];  end, end
if (nargin>=5)        % Scanner pixels per mm.
   dscan =   abs(str2num(varargin{5}));  % Scanner pixels per mm.
   sincpwr = abs(str2num(varargin{6}));  % Sinc power: 3 for film scanner, 2 for std digital camera.
   end;      % sincpwr < 2 for Foveon sensor; >=4 for inexpensive flatbed scanners.
if (nargin>=7)        % Slimming.
   ksharp = str2num(varargin{7});  % Boost constant: 1-ksharp*cos(pi*f/fboost)
   if (abs(ksharp)>eps) descript = [descript '; ksharp = ' num2str(ksharp,3) '      '];
      rsharp = 1;  if (ksharp<0) rsharp = 2;  end;  end;  % rsharp is sharpening radius.
   ksharp = abs(ksharp);
end
if (fboost>.45*f50)   % Boost by shifting peak.
   error(['fboost must be < .45*f50 = ' num2str(.45*f50,4)]);  end
close all
xmax = 0.5;    % Maximum of x1 in mm. The physical length of the image.
ymax = 200;   % Maximum line pairs per mm (spatial frequency) the plots.
ncyc = log10(100);      % Cycles of frequency to display. Argument is frequency ratio.
ymin = ymax/10^ncyc;  ymid = ymax/10^(ncyc/2);  % Minimum, middle spatial frequencies.
lfft = 2048;   % Total array length;
ldisp = 1920;  % Length to display;
lgr = lfft-ldisp;  % Length of gray needed to avoid cyclic problems.
ld1 = ldisp-1;  lf2 = lfft/2;
x1 = ([0:ld1]/ld1)*xmax;      % Distance in mm along axis to display.
xinc = x1(2)-x1(1);           % x1 (spatial) increment in mm.
fdel = 1/(lfft*xinc);         % Spatial frequency increment in 1/mm.
fq1 = 10.^(ncyc*[0:ld1]/ld1);  fq1 = ymax*fq1/max(fq1);  % Spatial frequency 1/mm.

xint = 1/ld1;  % Substitute for xinc
% figure;  semilogy([0:ld1]*xint,fq1);  grid on;  zoom on;
% ax = axis;  ax(2) = ld1*xint;  ax(4) = ymax;  axis(ax);
% title('Spatial frequency as a function of distance');
% xlabel('Normalized distance');  ylabel('Spatial frequency in line pairs per mm (1/mm)');
% p1 = get(gcf,'Position');  p1(4) = .8*p1(4);  p1(3) = .8*p1(3);  set(gcf,'Position',p1);

u1 = 2.*pi*cumsum(fq1)*xinc;  % Equation for varying frequency is tricky!
ycos = .5*(1+cos(u1));  % Original cosine pattern before rounding.
temp = find(round(ycos(1:ldisp-1))~=round(ycos(2:ldisp)));  % Anti-aliasing calculation.
temp1 = (.5-ycos(temp))./(ycos(temp+1)-ycos(temp));  % Fraction for crossing.
ybw = round(ycos);  % {0,1} Basic black/white pattern to analyze.
ybw(temp+round(temp1)) = (1-abs(temp1-.5)).*ybw(temp+  round(temp1)) + ...
                         (abs(temp1-.5))  .*ybw(temp+1-round(temp1)); clear temp,temp1;
% figure;  plot(x1,ycos,x1,ybw);  zoom on; % For test.  %%%%%%%%%%%%%%%%%%%
half = .5*ones(1,lgr);  ybw = [ybw half];  ycos = [ycos half];  % Last lgr points have a value of 1/2.


if (prtarg)  %              COUNTERPARTS FOR PRINTER
   xpmax = 10;     % Max x1 mm for printer test image.
   lpfft = 16384; lpdsp = 15360; lpgr = lpfft-lpdsp;  % Counterparts of lfft... for printer.
   lpd1 = lpdsp-1;  lpf2 = lpfft/2;
   xp1 = ([0:lpd1]/lpd1)*xpmax;  % Distance in mm along axis to display.
   xpinc = xpmax/lpd1;           % Spatial increment in mm.
   fpel = 1/(lpfft*xpinc);
   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.
   temp = find(round(ypw(1:lpdsp-1))~=round(ypw(2:lpdsp)));  % Anti-aliasing calculation.
   temp1 = (.5-ypw(temp))./(ypw(temp+1)-ypw(temp));  % Fraction for crossing.
   ypw = round(ypw);  % {0,1} basic black/white pattern to analyze. Comment out for testing.
   ypw(temp+round(temp1)) = (1-abs(temp1-.5)).*ypw(temp+  round(temp1)) + ...
                            (abs(temp1-.5))  .*ypw(temp+1-round(temp1));
   clear temp, temp1;
   % figure;  plot(1:lpdsp,ypw,'bx-');  zoom on;     %%%%%%%
   % ax = axis;  ax(3) = -1;  ax(4) = 2;  axis(ax);  %%%%%%%
   phalf = .5*ones(1,lpgr);  ypw = [ypw phalf];  % Last lgr points have a value of 1/2.
end  %                     END PRINTER FOR NOW.

cmap = zeros(256,3);  % Color map: gray scale. Apparently gamma corrected.
gamma = 1;  cmap(1:256,1) = ([0:255]'/255).^(1/gamma);
cmap(1:256,2) = cmap(1:256,1);  cmap(1:256,3) = cmap(1:256,1);

lvon = 40;  lv2 = ceil(.4*lvon);  % Lvon is pattern height. lv2 is text position.
vones = ones([lvon 1]);  % Start sequence display.
yp2 = real(ybw(1:ldisp));  yp3 = vones*yp2;  % Original bands. yp3 is a 2D matrix.
yp3 = yp3-min(min(yp3));  yp3 = yp3/max(max(yp3));  % Normalize plot.
ypp = yp3;  % Original bands

f1 = fft(ybw);  f1cos = fft(ycos);  fq2 = [0:lf2, -lf2+1:-1];
% figure; plot(fq2,real(f1),fq2,imag(f1));  grid on;  zoom on;
tfrfn = 1./(1+(fq2*fdel/f50).^2);  % MTF (without cos boost).
if (prtarg)  fp1 = fft(ypw); fp2 = [0:lpf2,-lpf2+1:-1];   % For printer.
   pfrfn = 1./(1+(fp2*fpel/f50).^2);  end  % MTF (without cos boost).
% Add cosine boost, but not at high frequencies. Use tricks.
if (fboost>eps)  % MTF with shifted zero point.
   tfrfn = 1./(1+((abs(fq2*fdel)-fboost)/sqrt(f50^2-2*f50*fboost-fboost^2)).^2);
   tfrfn = tfrfn/tfrfn(1);  % normalize tfrfn to lowest frequency.
   if (prtarg)
      pfrfn = 1./(1+((abs(fp2*fpel)-fboost)/sqrt(f50^2-2*f50*fboost-fboost^2)).^2); % Printer
      pfrfn = pfrfn/pfrfn(1);  end  % normalize tfrfn to lowest frequency.
end

if (flens>eps)                           % LENS
   tlens = 1./(1+abs(fq2*fdel/flens).^lord);  % Lens MTF. Order (lord) defaults to 2.
   if (prtarg)  plens = 1./(1+abs(fp2*fpel/flens).^lord);  end  % Lens MTF for printer.
   f2 = ifft(tlens.*f1);  % Lens only.
   yp2 = real(f2(1:ldisp));  yp3 = vones*yp2;  % Lens only. yp3 is a 2D matrix.
   yp3 = yp3-min(min(yp3));  yp3 = yp3/max(max(yp3));  % Normalize plot.
   ypp  = [ypp ; yp3];       % LENS
   f2 = ifft(tfrfn.*f1);
   if (f50>eps & f50<999)        % Omit film.
      yp2 = real(f2(1:ldisp));  yp3 = vones*yp2;  % Lens only. yp3 is a 2D matrix.
      yp3 = yp3-min(min(yp3));  yp3 = yp3/max(max(yp3));  % Normalize plot.
      ypp  = [ypp ; yp3];    % FILM
   end
   tfrfn = tfrfn.*tlens;     % Make tfrfn lens + film.
   if (prtarg) pfrfn = pfrfn.*plens;  end
end  % Multiply tfrfn by lens MTF.

[t50,i50] = min((tfrfn(1:lf2)-.5).^2);  fq50 = fq2(i50)*fdel;  % Find 50% point.
if (tfrfn(i50)>.5) i51 = i50+1;  else i51 = i50-1;  end
f51 = abs((tfrfn(i50)-.5)/(tfrfn(i50)-tfrfn(i51)));  % Fraction i51 for interpolation.
fq50 = fdel*((1-f51)*fq2(i50)+f51*fq2(i51)); % 50% point by interpolation.
[t10,i10] = min((tfrfn(1:lf2)-.1).^2);  fq10 = fq2(i10)*fdel;  % Find 10% point.
if (tfrfn(i10)>.1) i11 = i10+1;  else i11 = i10-1;  end
f11 = abs((tfrfn(i10)-.1)/(tfrfn(i10)-tfrfn(i11)));  % Fraction i11 for interpolation.
fq10 = fdel*((1-f11)*fq2(i10)+f11*fq2(i11)); % 10% point by interpolation.

% figure; loglog(fq2(2:lf2)*fdel,100*tfrfn(2:lf2));  grid on;  zoom on;
% ax = axis;  ax(1) = ymin;  ax(2) = ymax;  ax(3) = 1;
% if (max(tfrfn)<2 & ax(4)>200) ax(4) = 200;  end;  axis(ax);
xlblm = ['Line pairs per mm;  MTF = 50%,10% @ ' num2str(fq50,3) ', ' num2str(fq10,3) '/mm'];
if (fboost>eps) xlblm = [xlblm ';  Max MTF = ',num2str(max(tfrfn),3)];  end
% xlabel(xlblm);  title(descript);

f2 = ifft(tfrfn.*f1);  f2cos = ifft(tfrfn.*f1cos);
yp2 = real(f2(1:ldisp));  yp3 = vones*yp2;  % yp3 is a 2D matrix.  Lens+film
yp3 = yp3-min([min(yp3) 0]);  yp3 = yp3/max([max(yp3) 1]);  % Normalize plot.
if (f50>eps & f50<999) ypp  = [ypp ; yp3];  end       % LENS + FILM

if (prtarg)fp2 = ifft(pfrfn.*fp1);  %%%%%%%%%%%%%%% FOR PRINTER %%%%%%%%%%%%
   figure; plot([0:lpd1]*xpinc,real(fp2(1:lpdsp)),[0:lpd1]*xpinc,imag(fp2(1:lpdsp)));  %%%%%%
   ax = axis;  ax(2) = lpd1*xpinc;  axis(ax);  %%%%%%%
   title('For printer');  grid on;  zoom on;   %%%%%%%
end

y6 = real(f2);  y6cos = real(f2cos(1:ldisp));
y6 = y6-min([min(y6(1:ldisp)) 0]);  y6 = y6/max([max(y6(1:ldisp)) 1]);  % Normalize plot.
y6cos = y6cos-min([min(y6cos(1:ldisp)) 0]);  y6cos = y6cos/max([max(y6cos(1:ldisp)) 1]);
y7 = ones(12,1)*ycos(1:ldisp);
y7 = [y7; ones(38,1)*y6cos];  % y7(13:50,:)
y7 = [y7; ones(12,1)*ybw(1:ldisp)];    % y7(13:50,:)
y7 = [y7; ones(38,1)*y6(1:ldisp)];     % y7(63:100,:)
imwrite(y7,'MTFcurve_image_1.tif','tif');

% COMBINED PLOT
figure; p1 = get(gcf,'Position');
p1(2) = p1(2)-60;   p1(4) = p1(4)+60; p1(3) = p1(3)-80;  set(gcf,'Position',p1);
p2 = get(gcf,'PaperPosition');
p2(2) = p2(2)-1;  p2(4) = p2(4)+2;  set(gcf,'PaperPosition',p2);

subplot('Position',[.10 .4 .88 .55]);  % [l b w h]
image(255*y7);  colormap(cmap);  % get(gcf) gets properties.
h = gca;  set(h,'YTickLabel',{' '});  set(h,'XTickLabel',{' '});
title(descript);  ylabel('Bar pattern (bottom); Sine pattern (top)');
htxt = text(ldisp/20, 6,['Sine pattern: Original']);
set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
htxt = text(ldisp/20,56,['Bar pattern: Original']);
set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
if (flens<eps)
   htxt = text(ldisp/20,69,['Bar pattern: Film only']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   htxt = text(ldisp/20,19,['Sine pattern: Film only']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
elseif (f50<eps | f50>5000)
   htxt = text(ldisp/20,69,['Bar pattern: Lens only']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   htxt = text(ldisp/20,19,['Sine pattern: Lens only']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
else
   htxt = text(ldisp/20,69,['Bar pattern: Film & lens']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   htxt = text(ldisp/20,19,['Sine pattern: Film & lens']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
end

subplot('Position',[.10 .09 .88 .31]);   % [l b w h]
loglog(fq2(2:lf2)*fdel,100*tfrfn(2:lf2),'b-', 'LineWidth',1.5);  grid on;  zoom on;  hold on;
loglog(fq1,10.^(2*y6(1:ldisp)),'r-', 'LineWidth',1.5);
if (flens>eps) loglog(fq2(2:lf2)*fdel,100*tlens(2:lf2),'b:');  end
plot([2 200],[10 10],'r--');  plot([2 200],[50 50],'b:');  hold off;
ax = axis;  ax(1) = ymin;  ax(2) = ymax;  % ax(3) = 1; % old
ax(3) = min(1,  min(10.^(2*y6(1:ldisp))));
ax(4) = max([100 max(10.^(2*y6(1:ldisp))) max(100*tfrfn(2:lf2))]);  axis(ax);
xlabel(xlblm);  ylabel('MTF %, amplitude');

if (dscan>2)  % Results for a scanner with dscan pixels/mm (dscan*25.4 pixels/inch)
   tfrtot = tfrfn;  % Total transfer function with scanner.
   sscan = 1/(dscan*xinc);  % Scan steps in array points (non-integer).
   sstep = [1:sscan:lfft];  % Scan boundaries for averaging.
   if (sstep(length(sstep))<lfft) sstep = [sstep lfft];  end
   lscan = floor(sscan);
   if (sincpwr<2.5) lspat = 2*lscan+1;  % sinc^2: Length of sensitivity pattern.
      lsphf = lscan;  % Half the sensitivity pattern.
      senspat = zeros(1,lspat);  % Sensitivity pattern (Windowing = triangle = sinc^2)
      senspat(1:lscan+1) = [1:lscan+1]/(lscan+1);  % 1st half of triangle.
      senspat(lscan+2:lspat) = [lscan:-1:1]/(lscan+1);  % 2nd half.
   else  % dscan entered as a negative number: |sinc|^3 for film scanner.
      sens1 = ones(1,lscan);
      senspat = conv(sens1,sens1);  senspat = conv(senspat,sens1);  % 2 convolutions.
      lspat = length(senspat);  lsphf = floor(lspat/2)-1;
   end
   fstep = zeros([1 length(sstep)-1]);  lfs = length(fstep);
   lton = 60;  % Length of tilted plots.
   fscan = ones([lton 1])*y6;  % The scanned image.
   fsharp = fscan;  % The sharpened image.
   
   for (i2=1:lton)  % Tilted edge inside this loop (the slow one).
      sshft = sscan*(i2-1)/(lton-1);
      for (i1=1:lfs)  % Use windowing function for convolution.
         xs2 = sstep(i1)+sscan/2-sshft;  % Midpoint.
         xs2 = min(xs2,lfft);  xs2 = max(xs2,1);
         xs1 = max(xs2-lsphf,1);  xs3 = min(xs2+lsphf,lfft);
         % Convolve with y6.
         fstep(i1) = mean(senspat(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf)) ...
            .*y6(round(xs1):round(xs3)))/ ...
            mean(senspat(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf))); % Correct at ends.
         ix0 = max(sstep(i1)-sshft,1);  % Non-integer left boundary for anti-aliasing.
         ix1 = ceil(ix0);  ix2 = max(floor(sstep(i1+1)-sshft),1);
         fscan(i2,ix1:ix2) = fstep(i1);
         if (i1>1 & ix1>1)  % fscan MUST be antialiased to avoid the jaggies!
            fx0 = ix1-ix0;  % Weighting fraction.
            fscan(i2,ix1-1) = fx0*fstep(i1)+(1-fx0)*fstep(i1-1);
         end
      end
      if (ksharp>eps)  % Sharpen the image using a cosine boost. rsharp = radius.
         fshp = fstep;
         if (rsharp==1)
            fshp(2:lfs-1) = (fstep(2:lfs-1)-.5*ksharp*(fstep(1:lfs-2)+fstep(3:lfs)))/(1-ksharp);
         else   % rsharp == 2.
            fshp(3:lfs-2) = (fstep(3:lfs-2)-.5*ksharp*(fstep(1:lfs-4)+fstep(5:lfs)))/(1-ksharp);
         end
         for (i1=1:lfs)
            ix0 = max(sstep(i1)-sshft,1);  % Non-integer left boundary for anti-aliasing.
            ix1 = ceil(ix0);  ix2 = max(floor(sstep(i1+1)-sshft),1);
            fsharp(i2,ix1:ix2) = fshp(i1);
            if (i1>1 & ix1>1)  % fscan MUST be antialiased to avoid the jaggies!
               fx0 = ix1-ix0;  % Weighting fraction.
               fsharp(i2,ix1-1) = fx0*fshp(i1)+(1-fx0)*fshp(i1-1);
            end
         end
      end  % End cosine boost for sharpening
   end
   
   if (prtarg)  %                 PRINTER TARGET
   pscan =1/(dscan*xpinc);   % Scan steps in points in array fp2 (non-integer).
   pstep = [1:pscan:lpdsp];  % Scan boundaries in array fp2 for averaging.
   lps = length(pstep);      % Number of points (more than target length, typ. 10mm)
   if (pstep(length(pstep))<lpdsp) pstep = [pstep lpdsp];  end
   lpst = length(pstep);
   lpscan = floor(pscan);
      if (sincpwr<2.5) lsprt = 2*lpscan+1;  % sinc^2: Length of sensitivity pattern.
      lsphf = lpscan;  % Half the sensitivity pattern.
      sensprt = zeros(1,lsprt);  % Sensitivity pattern (Windowing = triangle = sinc^2)
      sensprt(1:lpscan+1) = [1:lpscan+1]/(lpscan+1);  % 1st half of triangle.
      sensprt(lpscan+2:lsprt) = [lpscan:-1:1]/(lpscan+1);  % 2nd half.
   else  % dscan entered as a negative number: |sinc|^3 for film scanner.
      sens1 = ones(1,lpscan);
      sensprt = conv(sens1,sens1);  sensprt = conv(sensprt,sens1);  % 2 convolutions.
      lsprt = length(sensprt);  lsphf = floor(lsprt/2)-1;
   end
   fpstep = zeros([1 lpst-1]);  % Sampled scanner array.
   fporig = zeros([1 lpst-1]);  % Original scanner array (no film, lens).
   y6 = real(fp2);  % Recycle y6; was real(f2).
   for (i1=1:lps)  % Use windowing function for convolution.
      xs2 = pstep(i1)+pscan/2;  % Midpoint.
      xs2 = min(xs2,lpdsp);  xs2 = max(xs2,1);
      xs1 = max(xs2-lsphf,1);  xs3 = min(xs2+lsphf,lpdsp);
      fpstep(i1) = mean(sensprt(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf)) ...
         .*y6(round(xs1):round(xs3)))/ ...
         mean(sensprt(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf))); % Correct at ends.
      fporig(i1) = mean(sensprt(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf)) ...
         .*ypw(round(xs1):round(xs3)))/ ...
         mean(sensprt(round(xs1+1-xs2+lsphf):round(xs3+1-xs2+lsphf))); % Correct at ends.
      ix0 = max(pstep(i1),1);  % Non-integer left boundary for anti-aliasing.
      ix1 = ceil(ix0);  ix2 = max(floor(pstep(i1+1)),1);
   end
   if (ksharp>eps)  % Sharpen here using a cosine boost.
      if (rsharp==1)
         fpstep(2:lps-1) = (fpstep(2:lps-1)-.5*ksharp*(fpstep(1:lps-2)+fpstep(3:lps)))/(1-ksharp);
      else   % rsharp == 2.
         fpstep(3:lps-2) = (fpstep(3:lps-2)-.5*ksharp*(fpstep(1:lps-4)+fpstep(5:lps)))/(1-ksharp);
      end
   end  % End cosine boost for sharpening
   htarg = round(lpst/16);  ht1 = round(lpst/72);  ht2 = round(lpst/160);
   fptarg = ones(6*htarg,lpst-1);
   m1 = min(fpstep);  m2 = max(fpstep);
   for (i1=1:2*htarg-1)  fptarg(i1,:) = ...
         ((fpstep-m1)/(m2-m1)*(i1-1)+.5*(2*htarg-i1))/(2*htarg-1);  end
   for (i1=2*htarg+7:3*htarg-1)    fptarg(i1,:) = (fpstep-m1)/(m2-m1);  end
   for (i1=3*htarg+7:4*htarg)      fptarg(i1,:) = fpstep;  end
   grayw = max(round(htarg/8),6);  % Width of gray line.
   for (i1=2*htarg:2*htarg+grayw)  fptarg(i1,:) = .5*ones(1,lpst-1);  end
   for (i1=3*htarg:3*htarg+grayw)  fptarg(i1,:) = .5*ones(1,lpst-1);  end
   % Create tics for 2 to 200 lp/mm logarithmic scale.
   xt1 = log10(2);  xt3 = (lpst-2)/(log10(200)-xt1);
   tic1 = [2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 200];
   tic2 = [2.5 3.5 15 25 35 150];
   ltic1 = round(1+xt3*(log10(tic1)-xt1));
   ltic2 = round(1+xt3*(log10(tic2)-xt1));
   for (i1=4*htarg+1:4*htarg+1+ht1)    fptarg(i1,ltic1) = 0;  end
   for (i1=4*htarg+1:4*htarg+1+ht2)    fptarg(i1,ltic2) = 0;  end
   for (i1=5*htarg-ht1:5*htarg)        fptarg(i1,ltic1) = 0;  end
   for (i1=5*htarg-ht2:5*htarg)        fptarg(i1,ltic2) = 0;  end
   for (i1=5*htarg+1:6*htarg)  fptarg(i1,:) = fporig;  end
   imwrite(fptarg,'printarget.tif','tif');  % Write the target.
   fprintf(' Printer target has been written to printarget.tif.\n');
   figure;  subplot('Position',[.10 .44 .88 .52]);  % [l b w h]
   image(255*fptarg);  colormap(cmap);
   subplot('Position',[.10 .10 .88 .34]);  plot(fpstep);  zoom on;  grid on;
   ax = axis;  ax(2) = length(fpstep);  axis(ax);
   %          Create a super target for photographing, etc.
   if (0)
      clear fptarg  % Get some memory. This target is very long.
      fptarg = zeros(200,lpdsp);
      for (i1=1:200)  fptarg(i1,:) = y6(1:lpdsp);  end
      imwrite(fptarg,'idealtgt.tif','tif');  % Write the target.
   end
   end   %               END PRINTER TARGET
   
   if (ksharp<eps)  % No sharpening.
      yp3 = fscan(:,1:ldisp);
      yp3 = yp3-min(min(yp3));  yp3 = yp3/max(max(yp3));  % Normalize plot.
      ypp  = [ypp ; yp3];
   else             % Sharpen the image using a cosine boost with appropriate radius, rsharp.
      yp3 = [fscan(:,1:ldisp); fsharp(:,1:ldisp)];
      yp3 = yp3-min(min(yp3));  yp3 = yp3/max(max(yp3));  % Normalize plot.
      ypp  = [ypp ; yp3];
      if (1)  % Illustration of sharpening. if (1) turn on; if (0) turn off.
         temp1 = real(ifft(tfrtot));  temp1 = temp1./max(temp1);
         lt = length(temp1);  sr = round(sscan*rsharp);  lt2 = lt-2*sr;
         lt3 = lt2/2-8*sr;  lt4 = lt2/2+8*sr; xtemp = [lt3-lt2/2:lt4-lt2/2]/sscan;
         temp1 = [temp1(lt/2+1:lt) temp1(1:lt/2)];  % Put peak in center.
         temp2 = -.5*ksharp*temp1(1+2*sr:lt);  temp3 = -.5*ksharp*temp1(1:lt-2*sr);
         temp4 = (temp1(1+sr:lt-sr)+temp2+temp3)/(1-ksharp);  % Slimmed pulse.
         temp5 = cumsum(temp1(1+sr:lt-sr));  temp5 = temp5/temp5(length(temp5));  % Edge normlzd to 1.
         temp6 = cumsum(temp4);  temp6 = temp6/temp6(length(temp6));  % Edge normlzd to 1.
         figure;   subplot('Position',[.2 .2 .6 .6]);  % [l b w h]

         plot(xtemp,temp1(lt3+sr:lt4+sr),'k-', xtemp,temp4(lt3:lt4)/max(temp4),'r-', ...
            xtemp,temp2(lt3:lt4),'b:', xtemp,temp3(lt3:lt4),'b:');  hold on;
         plot(xtemp,temp5(lt3:lt4),'k', xtemp,temp6(lt3:lt4),'r', 'LineWidth',2);  hold off;
         xlabel('Spacing in pixels'); ylabel('Relative luminance');  axis tight;  zoom on;
         ax = axis; ax(1) = -5; ax(2) = 5; ax(3) = ax(3)-.05; ax(4) = ax(4)+.1; axis(ax);
      end  % End sharpening illustration.
      tsharp= (1-ksharp*cos(2*pi*rsharp*fdel*fq2/dscan))./(1-ksharp);  % Sharpening MTF.
      tfrtot= tfrtot.*tsharp;  lfs = length(fstep);  fscan = fsharp;
   end  % End cosine boost for sharpening
   
   % y7 is the composite image to plot.
   y7  =  ones(20,1)*ycos(1:ldisp);         % Original sine.
   y7 = [y7; ones(20,1)*ybw(1:ldisp)];      % Original bar.
   fscann = real(fscan(:,1:ldisp));    % Scanned bar.
   y7mn = min(min(fscann)); y7mn = min(y7mn,0);  fscann = fscann-y7mn; 
   y7mx = max(max(fscann)); y7mx = max(y7mx,1);  fscann = fscann/y7mx;  % Normalize the plot.
   y7 = [y7; ones(20,1)*fscann(1,:)];  % Scanned bar.
   y7 = [y7; fscann];  % Combine tilted and gradient plots. Height = 60.
   imwrite(y7,'MTFcurve_image_2.tif','tif');
   
   ffscan = fft(fscan(1,:));  tfrscan = abs(ffscan./f1);  % MTF of scan.
   tfrsc = abs(sinc(fdel*fq2/dscan)).^sincpwr;  % Scanner transfer function (sin(x)/x)^2.
   tfrscp = max(abs(tfrsc),.0001);  % |Scanner tfr fn| to plot.
   tfrtot = tfrtot.*tfrscp;         % Total transfer function to plot.

   [t50,i50] = min((tfrtot(1:lf2)-.5).^2);  fq50 = fq2(i50)*fdel;  % Find 50% point.
   if (tfrtot(i50)>.5) i51 = i50+1;  else i51 = i50-1;  end
   f51 = abs((tfrtot(i50)-.5)/(tfrtot(i50)-tfrtot(i51)));  % Fraction i51 for interpolation.
   fq50 = fdel*((1-f51)*fq2(i50)+f51*fq2(i51)); % 50% point by interpolation.
   [t10,i10] = min((tfrtot(1:lf2)-.1).^2);  fq10 = fq2(i10)*fdel;  % Find 10% point.
   if (tfrtot(i10)>.1) i11 = i10+1;  else i11 = i10-1;  end
   f11 = abs((tfrtot(i10)-.1)/(tfrtot(i10)-tfrtot(i11)));  % Fraction i11 for interpolation.
   fq10 = fdel*((1-f11)*fq2(i10)+f11*fq2(i11)); % 10% point by interpolation.
   xlblm = ['Line pairs per mm;  MTF = 50%,10% @ ' num2str(fq50,3) ', ' num2str(fq10,3) ' lp/mm'];
   if (fboost>eps) xlblm = [xlblm ';  Max MTF = ',num2str(max(tfrfn),3)];  end
   
   % figure;  loglog(fq2(2:lf2)*fdel,100*tfrscan(2:lf2));  grid on;  zoom on;
   % hold on; loglog(fq2(2:lf2)*fdel,100*abs(tfrsc(2:lf2)),'r-'); hold off;
   % ax(1) = ymin;  ax(2) = ymax;  ax(3) = .1; ax(4) = 1000.;
   % if (max(tfrscan)<2 & ax(4)>200) ax(4) = 200;  end;  axis(ax);
   % xlabel(xlblm);  title(descript);

   % COMBINED PLOT
   figure; p1 = get(gcf,'Position');
   p1(2) = p1(2)-60;   p1(4) = p1(4)+60; p1(3) = p1(3)-80;  set(gcf,'Position',p1);
   p2 = get(gcf,'PaperPosition');
   p2(2) = p2(2)-1;  p2(4) = p2(4)+2;  set(gcf,'PaperPosition',p2);
   
   subplot('Position',[.10 .4 .88 .55]);  % [l b w h]
   image(255*y7);  colormap(cmap);  % get(gcf) gets properties.
   h = gca;  set(h,'YTickLabel',{' '});  set(h,'XTickLabel',{' '});
   htxt = text(ldisp/24, 9,['Sine pattern: Original']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   htxt = text(ldisp/24,29,['Bar pattern: Original']);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   htxt = text(ldisp/24,49,['Digitizer resolution = ' num2str(dscan,4) ' pixels/mm = ' ...
      num2str(dscan*25.4,4) '/in;  sinc\^' num2str(sincpwr,2)]);
   set(htxt,'Color',[1 0 0],'Fontsize',10,'FontWeight','bold')
   title(descript);  ylabel('Slanted, original bar pattern');
   
   subplot('Position',[.10 .09 .88 .31]);  % [l b w h]
   loglog(fq2(2:lf2)*fdel,100*tfrfn(2:lf2), 'b-');  hold on;
   loglog(fq2(2:lf2)*fdel,100*tfrtot(2:lf2),'k-','LineWidth',1.5);
   grid on;  zoom on;
   % tfrscp is scanner tfr fn; tsharp is sharpening tfr fn; tfrtot is total transfer fn.
   if (ksharp>eps) tsn = max(tfrtot(2:lf2))/max(tfrscp(2:lf2).*tsharp(2:lf2)); % Normalize.
      loglog(fq2(2:lf2)*fdel,100*tfrscp(2:lf2).*tsharp(2:lf2)*tsn,'b:','LineWidth',1.5);
   else loglog(fq2(2:lf2)*fdel,100*tfrscp(2:lf2),'b--', 'LineWidth',1.5);  end
   % Note f2 (lens+film MTF) and fscan nominally vary between 0 and 1.
   % fmax = max([.5 max(abs(f2(1:ldisp)-.5)) max(abs(fscan(1,1:ldisp)-.5))]);  % Norml. to 0.5.
   fmax = max([.5 max(abs(f2(1:ldisp)-.5))]);  % Norml. to 0.5.
   loglog(fq1,10.^((real(f2(1:ldisp))   -.5)/fmax+1),'m-');
   fmax = max([.5 max(abs(fscan(1,1:ldisp)-.5))]);  % Norml. to 0.5.
   loglog(fq1,10.^((real(fscan(1,1:ldisp))-.5)/fmax+1),'r-', 'LineWidth',1.5);
   plot([2 200],[10 10],'r--');  plot([2 200],[50 50],'b:');  hold off;
   ax = axis;  ax(1) = ymin;  ax(2) = ymax;  % ax(3) = 1; % old
   ax(3) = 1;  ax(4) = max([100 max(100*tfrfn) max(100*tfrtot)]);
   axis(ax);  xlabel(xlblm);  ylabel('MTF %, amplitude');
end

%                                      SEQUENCE PLOT
figure;  subplot('Position',[.01 .115 .98 .87]);
image(255*ypp ); axis off;  colormap(cmap);
imwrite(ypp,'MTFcurve_image_3.tif','tif');

ypt = lv2;  htxt = text(ldisp/30,ypt,['Original bands']);
set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold')

if (nargin>=3)
   ypt = ypt+lvon;  htxt = text(ldisp/30,ypt,['Lens only:  flens = ' num2str(flens,3) ...
      ';  lord = ' num2str(lord,2)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
   if (f50>eps & f50<999) 
      ypt = ypt+lvon;  htxt = text(ldisp/30,ypt,['Film only:  f50 = ' num2str(f50,3) ...
         ';  fboost = ' num2str(fboost,3)]);
      set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
      ypt = ypt+lvon;  htxt = text(ldisp/30,ypt,'Film + lens');
      set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
   end
elseif (f50>eps & f50<999)       % Film only; no lens.
   ypt = ypt+lvon;  htxt = text(ldisp/30,ypt,['Film:  f50 = ' num2str(f50,3) ...
      ';  fboost = ' num2str(fboost,3)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
end

if (nargin>=5) ypt = ypt+lvon;  htxt = text(ldisp/30,ypt,['Digitizer:  ' num2str(dscan,4) ...
      ' dpmm = ' num2str(dscan*25.4,4) ' dpi;  sinc\^' num2str(sincpwr,2)]);
set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');  end

if (nargin>=7) ypt = ypt+lton;  htxt = text(ldisp/30,ypt,['Sharpen:  ' num2str(ksharp,3)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');  end

subplot('Position',[.01 .11 .98 .005]);
loglog(fq1,ones(size(fq1)),'k-');
ax = axis;  ax(1) = ymin;  ax(2) = ymax;  axis(ax);
xlabel(xlblm);

fprintf('%s',7);  % Ring bell.

function y=sinc(x)
%SINC Sin(pi*x)/(pi*x) function.
y=ones(size(x));
i=find(x);
y(i)=sin(pi*x(i))./(pi*x(i));