function mtfcurve(varargin)  % version 04/30/2002
% Visualize the effects of MTF (modulation transfer function)
% in film, lenses, scanners, and sharpening using simulation.  Format:
% mtfcurve f50   % f50 is the 50% MTF density of film in lines/mm, typically 20 - 60.
% mtfcurve f50 fboost   % fboost is the frequency shift of the 2nd order equalizer.
% mtfcurve f50 fboost flens  % flens is the 50% MTF density of the lens.
% mtfcurve f50 fboost flens lord  % lord is lens rolloff order (default = 2).
% mtfcurve f50 fboost flens lord dscan  % dscan is scanner pixels per mm (dpmm).
% mtfcurve f50 fboost flens lord dscan ksharp  % ksharp is sharpening boost.
%
% mtfcurve 45 13 61 2  % 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.
% If dscan > 0, scanner MTF follows |sinc|^3 (appropriate for a film scanner).
% If dscan < 0, scanner MTF follows |sinc|^2 (appropriate for digital camera);
% If ksharp>0, sharpening equivalent to a cosine boost is applied, where ksharp is 
%    the amount of boost.
%    Equation: (1+ksharp*cos(2*pi*f/dscan))  %%%%%%%%%%%%%%
%
% Copyright (C) 2001 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>6)
   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) '/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) '/mm']; end, end
if (nargin>=3)        % Lens f50.
   flens =  str2num(varargin{3});  % Lens 50% density.
   descript = [descript '; flens = ' num2str(flens,4) '/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 = 3;       % Sinc power: 3 if input>0; 2 if <0.
   if ((str2num(varargin{5}))<0) sincpwr = 2;  end;  end;
if (nargin>=6)        % Slimming.
   ksharp = str2num(varargin{6});  % Boost constant: 1-ksharp*cos(pi*f/fboost)
   if (ksharp>eps) descript = [descript '; ksharp = ' num2str(ksharp,3) '      ']; end; 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 lp/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.

figure;  semilogy([0:ld1]*xinc,fq1);  grid on;  zoom on;
ax = axis;  ax(2) = ld1*xinc;  ax(4) = ymax;  axis(ax);
title('Spatial frequency as a function of distance');
xlabel('Distance in mm.');  ylabel('Spatial frequency in lp/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!
ybw = .5*(1+cos(u1));  % Original cosine pattern before rounding.
temp = find(round(ybw(1:ldisp-1))~=round(ybw(2:ldisp)));  % Anti-aliasing calculation.
temp1 = (.5-ybw(temp))./(ybw(temp+1)-ybw(temp));  % Fraction for crossing.
ybw = round(ybw);  % {0,1} basic black/white pattern to analyze. Comment out for testing.
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,ybw);  zoom on; % For test.
half = .5*ones(1,lgr); ybw = [ybw 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);

y4 = [0:.01:1];  y5 = y4'*ybw(1:ldisp) + (1-y4)'*.5*ones(size(ybw(1:ldisp))); % y4 is gradient.
figure;  image(255*y5);  colormap(cmap);
title('Original bands');
xlabel([num2str(ymin,4) ' line pairs per mm at left;  ' num2str(ymid,4) ' in middle;  ' ...
   num2str(ymax,4) ' at right']);
ylabel('Maximum contrast at bottom [0,1]; none at top [all 0.5]');
lvon = 40;  lv2 = ceil(.4*lvon);  vones = ones([lvon 1]);  % Start sequence display.
yp2 = real(ybw(1:ldisp));  yp3 = vones*yp2;  % Original bands. yp3 is a 2D matrix.
lyp = length(yp2);  yp4 = yp2(1:5:lyp);  yx1 = round(.76*lyp); % Insert reduced box.
yp3(6:25,yx1+1:yx1+length(yp4)+10) = vones(1:20)*[zeros(1,5) yp4 zeros(1,5)];
yp3(5:20:26,yx1+1:yx1+length(yp4)+10) = vones(1:2)*[zeros(1,length(yp4)+10)];
yp = yp3;

f1 = fft(ybw);  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)
   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.
   lyp = length(yp2);  yp4 = yp2(1:5:lyp);  yx1 = round(.76*lyp); % Insert reduced box.
   yp3(6:25,yx1+1:yx1+length(yp4)+10) = vones(1:20)*[zeros(1,5) yp4 zeros(1,5)];
   yp3(5:20:26,yx1+1:yx1+length(yp4)+10) = vones(1:2)*[zeros(1,length(yp4)+10)];
   yp = [yp; yp3];
   f2 = ifft(tfrfn.*f1);
   yp2 = real(f2(1:ldisp));  yp3 = vones*yp2;  % Lens only. yp3 is a 2D matrix.
   lyp = length(yp2);  yp4 = yp2(1:5:lyp);  yx1 = round(.76*lyp); % Insert reduced box.
   yp3(6:25,yx1+1:yx1+length(yp4)+10) = vones(1:20)*[zeros(1,5) yp4 zeros(1,5)];
   yp3(5:20:26,yx1+1:yx1+length(yp4)+10) = vones(1:2)*[zeros(1,length(yp4)+10)];
   yp = [yp; yp3];
   tfrfn = tfrfn.*tlens;
   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);
yp2 = real(f2(1:ldisp));  yp3 = vones*yp2;  % yp3 is a 2D matrix.  Lens+film
lyp = length(yp2);  yp4 = yp2(1:5:lyp);  yx1 = round(.76*lyp); % Insert reduced box.
yp3(6:25,yx1+1:yx1+length(yp4)+10) = vones(1:20)*[zeros(1,5) yp4 zeros(1,5)];
yp3(5:20:26,yx1+1:yx1+length(yp4)+10) = vones(1:2)*[zeros(1,length(yp4)+10)];
yp = [yp; yp3];
figure; plot([0:ld1]*xinc,real(f2(1:ldisp)),[0:ld1]*xinc,imag(f2(1:ldisp)));

title(descript);  ylabel('Band amplitude');
xlabel(['Distance in mm;  ' num2str(ymin,4) ' line pairs per mm at left;  ' num2str(ymid,4) ...
   ' in middle;  ' num2str(ymax,4) ' at right']);
ax = axis;  ax(2) = ld1*xinc;  axis(ax);
grid on;  zoom on;

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);  y7 = y4'*y6(1:ldisp) + (1-y4)'*.5*ones(size(y6(1:ldisp)));
y7 = min(y7,1);  y7 = max(y7,0);  % Limit the array to a range of {0,1}.
% 1/5 size box.
sy7 = size(y7);  y8 = y7(1:5:101, 1:5:sy7(2));  sy8 = size(y8);
y9 = zeros([sy8(1)+2 sy8(2)+20]); y9(2:sy8(1)+1,11:sy8(2)+10) = y8;  sy9 = size(y9);
yy1 = 11;  yx1 = round(.76*sy7(2));
y7(yy1+1:yy1+sy9(1),yx1+1:yx1+sy9(2)) = y9;

% 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]);  % subplot(2,1,1);
image(255*y7);  colormap(cmap);  % get(gcf) gets properties.
title(descript);  ylabel('Max contrast at bottom [0,1]; none at top [0.5]     ');
if (f50>10*ymax)
   htxt = text(ldisp/2,7,['Original image']);
elseif (flens<eps)
   htxt = text(ldisp/2,7,['Film only']);
else
   htxt = text(ldisp/2,7,['Film & lens']);
end   
set(htxt,'Color',[1 1 .8],'Fontsize',10,'HorizontalAlignment','center','FontWeight','bold')
subplot('Position',[.10 .09 .88 .31]);  % subplot(2,1,1);
loglog(fq2(2:lf2)*fdel,100*tfrfn(2:lf2),'b-');  grid on;  zoom on;  hold on;
if (flens>eps) loglog(fq2(2:lf2)*fdel,100*tlens(2:lf2),'b:');  end
loglog(fq1,10.^(2*y6(1:ldisp)),'r-');
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))));  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.
         fshp = fstep;
         fshp(2:lfs-1) = (fstep(2:lfs-1)-ksharp*(fstep(1:lfs-2)+fstep(3:lfs)))/(1-2*ksharp);
         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.
      fpstep(2:lps-1) = (fpstep(2:lps-1)-ksharp*(fpstep(1:lps-2)+ ...
         fpstep(3:lps)))/(1-2*ksharp);
   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]);
   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
   
   yp3 = fscan(:,1:ldisp);  syp3 = size(yp3);
   yp4 = yp3(1:2:syp3(1),1:5:syp3(2));  % Needs anti-aliasing.
   syp4 = size(yp4);  yx1 = round(.76*syp3(2));  % Reduced box.
   yp4(:,2:syp4(2)-1) = (yp3(1:2:syp3(1),6:5:syp3(2)-5)+yp3(1:2:syp3(1),5:5:syp3(2)-6) ...
      +yp3(1:2:syp3(1),4:5:syp3(2)-7)+yp3(1:2:syp3(1),7:5:syp3(2)-4) ...
      +yp3(1:2:syp3(1),8:5:syp3(2)-3))/5;  % Anti aliasing!
   yp3(6:5+syp4(1),yx1+1:yx1+syp4(2)) = yp4;
   yp3(6:35,yx1+1:yx1+5) =                    vones(1:30)*[zeros(1,5)];  % Box outline left
   yp3(6:35,yx1+1+syp4(2)+5:yx1+10+syp4(2)) = vones(1:30)*[zeros(1,5)];
   yp3(5:30:36,yx1+1:yx1+syp4(2)+10) = vones(1:2)*[zeros(1,syp4(2)+10)]; % Box outline top,bot
   yp = [yp; yp3];
   
   if (ksharp>eps)  % Sharpen the image using a cosine boost.
      yp3 = fsharp(:,1:ldisp);   yp4 = yp3(1:2:syp3(1),1:5:syp3(2));  % Needs anti-aliasing.
      % syp4 = size(yp4);  yx1 = round(.76*syp3(2));  % Reduced box.
      yp4(:,2:syp4(2)-1) = (yp3(1:2:syp3(1),6:5:syp3(2)-5)+yp3(1:2:syp3(1),5:5:syp3(2)-6) ...
         +yp3(1:2:syp3(1),4:5:syp3(2)-7)+yp3(1:2:syp3(1),7:5:syp3(2)-4) ...
         +yp3(1:2:syp3(1),8:5:syp3(2)-3))/5;  % Anti aliasing!
      yp3(6:5+syp4(1),yx1+1:yx1+syp4(2)) = yp4;
      yp3(6:35,yx1+1:yx1+5) =                    vones(1:30)*[zeros(1,5)];  % Box outline left
      yp3(6:35,yx1+1+syp4(2)+5:yx1+10+syp4(2)) = vones(1:30)*[zeros(1,5)];
      yp3(5:30:36,yx1+1:yx1+syp4(2)+10) = vones(1:2)*[zeros(1,syp4(2)+10)]; % Box outline top,bot
      yp = [yp; yp3];
      if (0)  % 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);  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 = -ksharp*temp1(1+2*sr:lt);  temp3 = -ksharp*temp1(1:lt-2*sr);
         temp4 = (temp1(1+sr:lt-sr)+temp2+temp3)/(1-2*ksharp);  % Slimmed pulse.
         figure;   subplot('Position',[.2 .2 .6 .6]);
         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');  % xtemp,temp4(lt3:lt4),'r:',
         xlabel('Spacing in pixels'); ylabel('Density (D)');  axis tight;  zoom on;
      end  % End sharpening illustration.
      tsharp= (1-2*ksharp*cos(2*pi*fdel*fq2/dscan))./(1-2*ksharp);  % Sharpening MTF.
      tfrtot= tfrtot.*tsharp;  lfs = length(fstep);  fscan = fsharp;
   end  % End cosine boost for sharpening
   
   y6 = real(fscan(1,:));
   y7 = y4'*y6(1:ldisp) + (1-y4)'*.5*ones(size(y6(1:ldisp))); % Gradient; for combined plot.
   y7 = [y7; fscan(:,1:ldisp)];  % Combine tilted and gradient plots.
   y7 = min(y7,1);  y7 = max(y7,0);  % Limit the array to a range of {0,1}.
   sy7 = size(y7);  y8 = y7(1:5:161, 1:5:sy7(2));  sy8 = size(y8);  % Small box.
   
   y8(:,2:sy8(2)-1) = (y7(1:5:sy7(1),6:5:sy7(2)-5)+y7(1:5:sy7(1),5:5:sy7(2)-6) ...
      +y7(1:5:sy7(1),4:5:sy7(2)-7)+y7(1:5:sy7(1),7:5:sy7(2)-4) ...
      +y7(1:5:sy7(1),8:5:sy7(2)-3))/5;  % Anti aliasing!
   
   y9 = zeros([sy8(1)+2 sy8(2)+20]); y9(2:sy8(1)+1,11:sy8(2)+10) = y8;  sy9 = size(y9);
   yy1 = 16;  yx1 = round(.76*sy7(2));
   y7(yy1+1:yy1+sy9(1),yx1+1:yx1+sy9(2)) = y9;
   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.
   
   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);
   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);

   % 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]);
   image(255*y7);  colormap(cmap);  % get(gcf) gets properties.
   htxt = text(ldisp/2,9,['Scanner resolution = ' num2str(dscan,4) ' pixels/mm = ' ...
      num2str(dscan*25.4,4) '/in;  sinc^' num2str(sincpwr,1)]);
   set(htxt,'Color',[1 1 .8],'Fontsize',10,'HorizontalAlignment','center','FontWeight','bold')
   title(descript);  ylabel('Max contrast at bottom [0,1]; none at top [0.5]     ');
   subplot('Position',[.10 .09 .88 .31]);
   loglog(fq2(2:lf2)*fdel,100*tfrfn(2:lf2),'b-',fq2(2:lf2)*fdel,100*tfrtot(2:lf2),'k-');
   grid on;  zoom on;  hold on;
   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:');
   else loglog(fq2(2:lf2)*fdel,100*tfrscp(2:lf2),'b:');  end
   % Note f2 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.
   loglog(fq1,10.^((real(f2(1:ldisp))   -.5)/fmax+1),'m-');
   loglog(fq1,10.^((real(fscan(1,1:ldisp))-.5)/fmax+1),'r-');
   % if (ksharp>eps) loglog(fq2(2:lf2)*fdel,100*tsharp(2:lf2),'g-');  end % Test only.
   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
yp = yp./(max(max(yp))-min(min(yp)));  yp = yp-min(min(yp));  % Normalize to overshoot.
figure;  subplot('Position',[.01 .115 .98 .87]);
image(255*yp); axis off;  colormap(cmap);

if (nargin>=6) htxt = text(ldisp/30,4*lvon+lton+lv2,['Sharpen:  ' num2str(ksharp,3)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');  end
if (nargin>=5) htxt = text(ldisp/30,4*lvon+lv2,['Scan:  ' num2str(dscan,4) ' dpmm = ' ...
      num2str(dscan*25.4,4) ' dpi;  sinc^' num2str(sincpwr,1)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');  end
if (nargin>=3) htxt = text(ldisp/30,lvon+lv2,['Lens only:  flens = ' num2str(flens,3) ...
      ';  lord = ' num2str(lord,2)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
   htxt = text(ldisp/30,2*lvon+lv2,['Film only:  f50 = ' num2str(f50,3) ...
      ';  fboost = ' num2str(fboost,3)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
   htxt = text(ldisp/30,3*lvon+lv2,'Film + lens');
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
else  % Film only; no lens.
   htxt = text(ldisp/30,lvon+lv2,['Film:  f50 = ' num2str(f50,3) ...
      ';  fboost = ' num2str(fboost,3)]);
   set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold');
end

htxt = text(ldisp/30,lv2,['Original bands']);
set(htxt,'Color',[1 .2 0],'Fontsize',10,'FontWeight','bold')
subplot('Position',[.01 .11 .98 .005]);
loglog(fq1,ones(size(fq1)),'k-');
ax = axis;  ax(1) = ymin;  ax(2) = ymax;  axis(ax);
xlabel(xlblm);

% WRITE TIF FILE!
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));