{$N+,E+,M 65500,0,655360}



PROGRAM OptiPro;

{***************************************************************************}

{*  PROGRAM:  OptiPro

{*  AUTHOR:  C1C David W. Croft

{*  ORGANIZATION:  CS-36, x4306

{*  PURPOSE:  This program demonstrates some optical processing techniques.

{*    It requires the files EGAVGA.BGI, STDHDR.TPU, and FFT.TPU to run.

{*    Of note is the procedure SampleTheorem2D.

{***************************************************************************}



uses

  crt,

  stdhdr,

  fft,

  graph;



const

  Forever: boolean = false;



VAR

  xr,yi, fxr,fyi: RecMat;

  lower,upper,nd,i,j,k: integer;

  alpha, bravo, charlie, delta, echo, foxtrot, golf, hotel: integer;

  option: char;



Procedure QuietPause;

{***************************************************************************}

{*  This procedure flushes the keyboard buffer and then waits for a keypress.

{***************************************************************************}

begin

  repeat

    if keypressed then

      option := readkey;

  until not(keypressed);

  readln;

end;



function ColorBar(x: integer): integer;

{***************************************************************************}

{*  This procedure reorders the color sequence to vary from dark to light.

{***************************************************************************}

var

  y: integer;

begin

  case x of

    0:  y := 0;

    1:  y := 8;

    2:  y := 1;

    3:  y := 9;

    4:  y := 2;

    5:  y := 10;

    6:  y := 3;

    7:  y := 11;

    8:  y := 5;

    9:  y := 13;

   10:  y := 6;

   11:  y := 4;

   12:  y := 12;

   13:  y := 14;

   14:  y := 7;

   15:  y := 15;

  end; {case}

  ColorBar := y;

end;



function sinc(x: extended):extended;

{***************************************************************************}

{*  Sin(x)/x.  Note:  this is not Sin(Pi*x)/(Pi*x).

{***************************************************************************}

var

  bob: extended;

begin

  if x = 0 then

    bob := 1

  else

    bob := sin(x)/x;

  sinc := bob;

end;



Procedure DrawFrames;

{***************************************************************************}

{*  This simply draws the window frames for procedures FillRect and

{*  FillRectSplit.  It called by procedure InitGraphics.

{***************************************************************************}

begin

  Rectangle(0, 0, 257, 257);

  Rectangle(382, 0, 639, 257);

end;



Procedure InitGraphics;

{***************************************************************************}

{*  Task:  This procedure switches from CRT mode text output to graphics.

{*  Note:  For this procedure to be generalized the last command,

{*    DrawFrames, would have to be deleted.

{*  Author:  C1C David W. Croft, CS-36, x4306

{***************************************************************************}

var

  GraphDriver: integer;

  GraphMode: integer;

  ErrorCode: integer;

begin

  DetectGraph(GraphDriver, GraphMode);

  InitGraph(GraphDriver, GraphMode, '');

  ErrorCode := GraphResult;

  if ErrorCode <> grOk then

    begin

      Writeln('Graphics error:  ',GraphErrorMsg(ErrorCode));

      Writeln('Program aborted...');

      Halt(1);

    end;

  DrawFrames;

end;



Procedure FillRect(xr, yi, fxr, fyi: RecMat; nd: integer);

{***************************************************************************}

{*  Task:  This procedure takes in two matrices of NDxND dimensions,

{*    with two sub-matrices each to represent the real and imaginary parts,

{*    and displays side by side a color-scaled graph of their magnitudes.

{*  Inputs:

{*    xr, yi:  the first matrix, its real and imaginary parts in matrix form

{*    fxr, fyi:  the second matrix

{*    nd:  the dimension of the matrices along a side

{***************************************************************************}

var

  blocksize: integer;

  alpha, bravo: integer;

  zm, fzm: RecMat;

  zbig, zwee, fzbig, fzwee: extended;

  szbig, szwee, sfzbig, sfzwee: string;

  scaleStr, fscaleStr: string;

  scale, fscale: extended;

  charlie: integer;

begin

  zbig := -1e4000; fzbig := -1e4000; zwee := 1e4000; fzwee := 1e4000;

  blocksize := 256 div nd;

  for alpha := 0 to nd-1 do

    for bravo := 0 to nd-1 do

      begin

        zm[alpha,bravo] :=sqrt(sqr(xr[alpha,bravo])+sqr(yi[alpha,bravo]) );

        fzm[alpha,bravo]:=sqrt(sqr(fxr[alpha,bravo])+sqr(fyi[alpha,bravo]) );

        if zm[alpha,bravo] > zbig then

          zbig := zm[alpha,bravo];

        if zm[alpha,bravo] < zwee then

          zwee := zm[alpha,bravo];

        if fzm[alpha,bravo] > fzbig then

          fzbig := fzm[alpha,bravo];

        if fzm[alpha,bravo] < fzwee then

          fzwee := fzm[alpha,bravo];

      end;

  str(zwee:7:4, szwee); str(zbig:7:4, szbig);

  str(fzwee:7:4, sfzwee); str(fzbig:7:4, sfzbig);

  scale := (zbig - zwee)/15;

  fscale := (fzbig - fzwee)/15;

  if scale < 0.01*zwee then

    scale := 0;

  if fscale < 0.01*fzwee then

    fscale := 0;

  SetTextStyle(0,1,1);

  SetFillStyle(SolidFill, 0);

  Bar(0, 260, 608, 349);

  for alpha := 0 to 15 do

    begin

      SetFillStyle(SolidFill, ColorBar(alpha));

      Bar(alpha*40, 300,  (alpha+1)*40 -1, 309);

      str((alpha*scale + zwee):5:2, scalestr);

      OutTextXY(alpha*40+8,260, scalestr);

      str((alpha*fscale + fzwee):5:2, fscalestr);

      OutTextXY(alpha*40+8,309, fscalestr);

    end;

  for alpha := 0 to nd-1 do

    for bravo := 0 to nd-1 do

      begin

        if scale = 0 then

          begin

            if zbig = 0 then

              charlie := 0

            else

              charlie := 15;

          end

        else

          charlie := trunc( (zm[alpha,bravo] - zwee)/ scale);

        SetFillStyle(SolidFill, ColorBar(charlie) );

        Bar(alpha*blocksize+1, bravo*blocksize+1, (alpha+1)*blocksize, (bravo+1)*blocksize );

        if fscale = 0 then

          begin

            if fzbig = 0 then

              charlie := 0

            else

              charlie := 15;

          end

        else

          charlie := trunc( (fzm[alpha,bravo] - fzwee)/ fscale);

        SetFillStyle(SolidFill, ColorBar(charlie));

        Bar(alpha*blocksize+383, bravo*blocksize+1,

          (alpha+1)*blocksize+382, (bravo+1)*blocksize );

      end;

end;



Procedure FillRectSplit(xr, yi: RecMat; nd: integer);

{***************************************************************************}

{*  Task:  This procedure takes in two matrices of NDxND dimensions

{*    and displays side by side a color-scaled graph of their values.

{*  Inputs:

{*    xr, yi:  the two matrices

{*    nd:  the dimension of the matrices along a side

{***************************************************************************}

var

  blocksize: integer;

  alpha, bravo: integer;

  zm, fzm: RecMat;

  zbig, zwee, fzbig, fzwee: extended;

  szbig, szwee, sfzbig, sfzwee: string;

  scaleStr, fscaleStr: string;

  scale, fscale: extended;

  charlie: integer;

begin

  zbig := -1e4000; fzbig := -1e4000; zwee := 1e4000; fzwee := 1e4000;

  blocksize := 256 div nd;

  for alpha := 0 to nd-1 do

    for bravo := 0 to nd-1 do

      begin

        zm[alpha,bravo] := xr[alpha,bravo];

        fzm[alpha,bravo] := yi[alpha,bravo];

        if zm[alpha,bravo] > zbig then

          zbig := zm[alpha,bravo];

        if zm[alpha,bravo] < zwee then

          zwee := zm[alpha,bravo];

        if fzm[alpha,bravo] > fzbig then

          fzbig := fzm[alpha,bravo];

        if fzm[alpha,bravo] < fzwee then

          fzwee := fzm[alpha,bravo];

      end;

  str(zwee:7:4, szwee); str(zbig:7:4, szbig);

  str(fzwee:7:4, sfzwee); str(fzbig:7:4, sfzbig);

  scale := (zbig - zwee)/15;

  fscale := (fzbig - fzwee)/15;

  if scale < 0.01*zwee then

    scale := 0;

  if fscale < 0.01*fzwee then

    fscale := 0;

  SetTextStyle(0,1,1);

  SetFillStyle(SolidFill, 0);

  Bar(0, 260, 608, 349);

  for alpha := 0 to 15 do

    begin

      SetFillStyle(SolidFill, ColorBar(alpha));

      Bar(alpha*40, 300,  (alpha+1)*40 -1, 309);

      str((alpha*scale + zwee):5:2, scalestr);

      OutTextXY(alpha*40+8,260, scalestr);

      str((alpha*fscale + fzwee):5:2, fscalestr);

      OutTextXY(alpha*40+8,309, fscalestr);

    end;

  for alpha := 0 to nd-1 do

    for bravo := 0 to nd-1 do

      begin

        if scale = 0 then

          begin

            if zbig = 0 then

              charlie := 0

            else

              charlie := 15;

          end

        else

          charlie := trunc( (zm[alpha,bravo] - zwee)/ scale);

        SetFillStyle(SolidFill, ColorBar(charlie) );

        Bar(alpha*blocksize+1, bravo*blocksize+1, (alpha+1)*blocksize, (bravo+1)*blocksize );

        if fscale = 0 then

          begin

            if fzbig = 0 then

              charlie := 0

            else

              charlie := 15;

          end

        else

          charlie := trunc( (fzm[alpha,bravo] - fzwee)/ fscale);



        SetFillStyle(SolidFill, ColorBar(charlie));

        Bar(alpha*blocksize+383, bravo*blocksize+1,

          (alpha+1)*blocksize+382, (bravo+1)*blocksize );

      end;

end;



procedure QuadSwitch(var fxr, fyi: RecMat; nd: integer);

{***************************************************************************}

{*  Task:  this procedure switches the quadrants of two NDxND matrices:

{*    Quad 1 to Quad 3 and Quad 2 to Quad 4.  It is particularly useful

{*    after perform a 2-dimensional Fast Fourier Transform or its inverse.

{*  Inputs:

{*    fxr, fyi:  the two matrices to be quad switched

{*    nd:  the dimensions of matrix along a side

{*  Author:  C1C David W. Croft, CS-36, x4306

{***************************************************************************}

var

  alpha, bravo: integer;

  charlie: integer;

  tempr, tempi: RecMat;

begin

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      begin

        tempr[alpha, bravo] := fxr[alpha, bravo];

        tempi[alpha, bravo] := fyi[alpha, bravo];

      end;

  charlie := nd div 2;

  for alpha := 0 to charlie - 1 do

    for bravo := 0 to charlie - 1 do

      begin

        fxr[alpha, bravo] := tempr[alpha + charlie, bravo + charlie];

        fyi[alpha, bravo] := tempi[alpha + charlie, bravo + charlie];

        fxr[alpha + charlie, bravo + charlie] := tempr[alpha, bravo];

        fyi[alpha + charlie, bravo + charlie] := tempi[alpha, bravo];

        fxr[alpha, bravo + charlie] := tempr[alpha + charlie, bravo];

        fyi[alpha, bravo + charlie] := tempi[alpha + charlie, bravo];

        fxr[alpha + charlie, bravo] := tempr[alpha, bravo + charlie];

        fyi[alpha + charlie, bravo] := tempi[alpha, bravo + charlie];

      end;

end;



Procedure LoadMatrix(var xr, yi: recMat; var nd: integer);

{***************************************************************************}

{*  Task:  this loads a 16x16 matrix designated by its real and imaginary

{*    parts with an arbitrary pattern.  It is called by procedure

{*    PhaseContrast.

{***************************************************************************}

begin

  nd := 16; {note:  must be a power of 2}

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      begin

        xr[alpha,bravo] := 0;

        yi[alpha,bravo] := 0;

      end;

  for alpha := 4 to 11 do

    for bravo := 4 to 11 do

      begin

        xr[alpha,bravo] := cos((bravo+1.0)*Pi) + 2;

        yi[alpha,bravo] := cos(bravo*Pi) + 2;

      end;

end;



Procedure PhaseContrast;

{***************************************************************************}

{*  Task:

{*    1.  Display the real and imaginary of an arbitrary optical signal.

{*    2.  Display its magnitude (visible to the eye -- no phase components).

{*    3.  Display the signal with an added uniform imaginary vector.

{*    4.  Display the magnitude of step 2 and the phase contrasted magnitude

{*        of step 3.

{***************************************************************************}

begin

  LoadMatrix(xr, yi, nd);

  {display original input:  real and imaginary}

  FillRectSplit(xr, yi, nd);

  readln;



  {display original input:  magnitude}

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      begin

        fxr[alpha,bravo] := 0;

        fyi[alpha,bravo] := 0;

      end;

  FillRect(xr,yi,fxr,fyi,nd);

  readln;



  {Display Corrupted:  real and imaginary}

  fxr := xr;

  fyi := yi;

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

     fyi[alpha,bravo] := fyi[alpha,bravo] + 2;

  FillRectSplit(fxr, fyi, nd);

  readln;



  {Display Original and Corrupted:  magnitude}

  FillRect(xr, yi, fxr, fyi, nd);

  readln;

end;



Procedure AmplitudeFilter;

{***************************************************************************}

{*  Task:

{*    1.  Display a signal and its fourier spatial transform.

{*    2.  Add noise and display the corrupted signal with its transform.

{*    3.  Display the filter and the fourier filtered corrupted signal.

{*    4.  Display the inverse fourier which is the recovered signal.

{*    5.  Repeat the above for a non-recoverable, non-filterable signal.

{***************************************************************************}

var

  delta, echo: integer;

begin

  nd := 16;

  for delta := 0 to 1 do

    begin

      echo := delta * 2;



      for alpha := 0 to nd - 1 do

        for bravo := 0 to nd - 1 do

          begin

            xr[alpha,bravo] := cos(alpha*Pi) + cos(bravo*Pi) + echo;

            yi[alpha,bravo] := 0;

          end;

      fxr := xr;

      fyi := yi;

      FFT2dCalc(fxr,fyi,nd,nd,0);

      QuadSwitch(fxr,fyi,nd);

      FillRect(xr, yi, fxr, fyi, nd);

      readln;



      for alpha := 0 to nd - 1 do

        for bravo := 0 to nd - 1 do

          xr[alpha,bravo] := xr[alpha,bravo] + 12;

      fxr := xr;

      fyi := yi;

      FFT2dCalc(fxr,fyi,nd,nd,0);

      QuadSwitch(fxr,fyi,nd);

      FillRect(xr, yi, fxr, fyi, nd);

      readln;



      {the filter}

      for alpha := 0 to nd - 1 do

        for bravo := 0 to nd - 1 do

          yi[alpha,bravo] := 1;

      for alpha := 7 to 8 do

        for bravo := 7 to 8 do

          yi[alpha,bravo] := 0;



      {the filtering}

      for alpha := 0 to nd - 1 do

        for bravo := 0 to nd - 1 do

          begin

            fxr[alpha,bravo] := fxr[alpha,bravo]*yi[alpha,bravo];

            fyi[alpha,bravo] := fyi[alpha,bravo]*yi[alpha,bravo];

          end;



      FillRect(yi, yi, fxr, fyi, nd);

      readln;



      xr := fxr;

      yi := fyi;

      FFT2dCalc(xr,yi,nd,nd,1);

      QuadSwitch(xr,yi,nd);

      FillRect(xr, yi, fxr, fyi, nd);

      readln;



    end;

end;



Procedure SanctifyIt(var xr: RecMat; nd: integer);

{***************************************************************************}

{*  Task:  expand a matrix of (nd div 2)x(nd div 2) dimensions into a matrix

{*    of NDxND dimensions but inserts zeroes (holes) where it was stretched.

{*    This matrix is used to demonstrate the blanks in a matrix which will

{*    be filled by the SampleTheorem2D procedure.

{***************************************************************************}

var

  temp: RecMat;

begin

  for alpha := 0 to (nd - 1) div 2 do

    for bravo := 0 to (nd - 1) div 2 do

      begin

        charlie := 2 * alpha;

        delta := 2 * bravo;

        xr[charlie+1, delta] := 0;

        xr[charlie, delta+1] := 0;

        xr[charlie+1, delta+1] := 0;

      end;

end;



Procedure StretchIt(var xr: RecMat; var nd: integer);

{***************************************************************************}

{*  Task:  this transforms a (ND div 2) by (ND div 2) matrix into a NDxND

{*    matrix by "stretching" it.  It will appear as the former matrix, just

{*    bigger.

{***************************************************************************}

var

  temp: RecMat;

begin

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      temp[alpha,bravo] := xr[alpha div 2, bravo div 2];

  xr := temp;

end;



Procedure SampleTheorem;

{***************************************************************************}

{*  Task:  this demonstrates the one-dimensional sampling theorem

{*    application on a sampled sine wave.  The user is asked how many samples

{*    he wants to take of the sine wave and then how many points that he

{*    wants the sampling theorem to evaluate between those sampled points.

{***************************************************************************}

type

  bigarray = array[0..199] of extended;

var

  bravo, charlie: extended;

  SampleData, OtherData, RecoveredData, DiffData: bigarray;

  sierra, samples: integer;

  hotel: extended;

begin

  clrscr;

  write('How many samples of the sine wave do you want to take?  ');

  readln(samples);

  write('How many intermediate values do you want estimate?  ');

  readln(delta);

  writeln;

  writeln('Sampled sine wave data:');

  sierra := delta div samples;

  for alpha := 0 to (delta - 1) div sierra do

    begin

      bravo := 2*Pi*alpha/(delta div sierra);

      charlie := sin(bravo);

      write(charlie:5:3,'  ');

      SampleData[alpha] := charlie;

    end;

  writeln;

  writeln;

  writeln('Actual intermediate values for sine wave:');

  for alpha := 0 to delta - 1 do

    begin

      bravo := Pi*(alpha+0.5)/(delta/2);

      charlie := sin(bravo);

      write(charlie:5:3,' ');

      OtherData[alpha] := charlie;

    end;

  writeln;

  writeln;

  writeln('Intermediate values retrieved using the Sampling Theorem:');

  for alpha := 0 to delta-1 do

    begin

      bravo := (alpha+0.5)/(delta * sierra);

      charlie := 0;

      for echo := 0 to (delta - 1) div sierra do

        begin

          hotel := pi*(delta*bravo - echo);

          charlie := charlie + SampleData[echo]

            * sinc(hotel);

        end;

      write(charlie:5:3,' ');

      RecoveredData[alpha] := charlie;

    end;

  writeln;

  readln;

end;



Procedure SampleTheorem2D(var xr: RecMat; nd: integer);

{***************************************************************************}

{*  TASK:  Given a formerly (ND div 2)x(ND div 2) matrix that was stretched

{*    into an NDxND matrix, use the sampling theorem to evaluate the correct

{*    values for the new intermediate points in the matrix.

{*  INPUTS:

{*    xr:  an NDxND matrix, formerly a (ND div 2)x(ND div 2) matrix that

{*      was expanded and with its holes filled with arbitrary values.

{*      Those arbitrary values will be discarded.

{*    nd:  the dimension of the matrix along a side.

{*  AUTHOR:  C1C David W. Croft, CS-36, x4306

{***************************************************************************}

var

  alphaX, alphaY, foxtrotX, foxtrotY, golfX, golfY: integer;

  bravoX, bravoY, hotelX, hotelY, indiaX, indiaY: extended;

  echoX, echoY: integer;

  bravo, charlie, julietX, julietY: extended;

  SampleData: RecMat;

  sierra, samples: integer;

  delta: integer;

begin

  SampleData := xr;

  sierra := 1;

  delta := nd;



  for alphaX := 0 to delta-1 do

    for alphaY := 0 to delta-1 do

    begin

      bravoX := (alphaX+0.5)/(delta * sierra);

      bravoY := (alphaY+0.5)/(delta * sierra);

      charlie := 0;

      for echoX := 0 to (delta - 1) div sierra do

        for echoY := 0 to (delta - 1) div sierra do

        begin

          hotelx := pi*(delta*bravoX - echoX);

          hotely := pi*(delta*bravoY - echoY);

          julietx := sinc(hotelx);

          juliety := sinc(hotely);

          charlie := charlie + SampleData[echoX,echoY]

            * julietx

            * juliety;

        end;

      xr[alphaX, alphaY] := charlie;

    end;

end;



Procedure Demo1SampleTheorem2D;

{***************************************************************************}

{*  Task:

{*    1.  Display an 8x8 sinusoidal matrix.

{*    2.  Display the matrix and its fourier transform.

{*    3.  Stretch the matrix out to a 16x16 matrix and display the empty

{*        holes along side a 16x16 matrix that was stretched and filled

{*        with the adjacent values.

{*    4.  Use the 2-D Sampling Theorem to fill the holes and display it.

{***************************************************************************}

begin

  InitGraphics;

  nd := 8;

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      begin

        xr[alpha,bravo] := sin(2*Pi*alpha/nd);

        yi[alpha,bravo] := 0;

      end;

  FillRectSplit(xr, yi, nd);

  QuietPause;



  fxr := xr;

  fyi := yi;

  FFT2DCalc(fxr,fyi,nd,nd,0);

  QuadSwitch(fxr, fyi, nd);

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;



  nd := 16;

  StretchIt(xr, nd);

  StretchIt(yi, nd);



  fxr := xr;

  fyi := yi;

  SanctifyIt(fxr, nd);

  SanctifyIt(fyi, nd);

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;



  fxr := xr;

  fyi := yi;

  SampleTheorem2D(fxr, nd);

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      fyi[alpha,bravo] := 0;

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;

end;



Procedure Demo2SampleTheorem2D;

{***************************************************************************}

{*  Task:

{*    1.  Display an 8x8 2-D sinc function matrix.

{*    2.  Display the matrix and its fourier transform.

{*    3.  Stretch the matrix out to a 16x16 matrix and display the empty

{*        holes along side a 16x16 matrix that was stretched and filled

{*        with the adjacent values.

{*    4.  Use the 2-D Sampling Theorem to fill the holes and display it.

{***************************************************************************}

begin

  InitGraphics;

  nd := 8;

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      begin

        xr[alpha,bravo] := sin(2*Pi*alpha/nd) + sin(2*Pi*bravo/nd);

        yi[alpha,bravo] := 0;

      end;

  FillRectSplit(xr, yi, nd);

  QuietPause;



  fxr := xr;

  fyi := yi;

  FFT2DCalc(fxr,fyi,nd,nd,0);

  QuadSwitch(fxr, fyi, nd);

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;



  nd := 16;

  StretchIt(xr, nd);

  StretchIt(yi, nd);



  fxr := xr;

  fyi := yi;

  SanctifyIt(fxr, nd);

  SanctifyIt(fyi, nd);

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;



  fxr := xr;

  fyi := yi;

  SampleTheorem2D(fxr, nd);

  for alpha := 0 to nd - 1 do

    for bravo := 0 to nd - 1 do

      fyi[alpha,bravo] := 0;

  FillRect(xr, yi, fxr, fyi, nd);

  QuietPause;

end;



{***************************************************************************}

{***************************************************************************}

BEGIN   {main program fft}

  repeat

    TextBackground(Black);

    clrscr;

    writeln('0.  Quit');

    writeln('1.  Sampling Theorem');

    writeln('2.  2-D Sampling Theorem Demo 1');

    writeln('3.  2-D Sampling Theorem Demo 2');

    writeln('4.  Phase Contrast');

    writeln('5.  Amplitude Filter');

    writeln;

    write('Option:  ');

    readln(option);

    if option = '0' then

      halt;

    if option = '1' then

      SampleTheorem

    else

      begin

        InitGraphics;

        nd := 16;

      end;

    case option of

      '2':  Demo1SampleTheorem2D;

      '3':  Demo2SampleTheorem2D;

      '4':  PhaseContrast;

      '5':  AmplitudeFilter;

      else

        write(#7);

    end; {case}

    RestoreCrtMode;

    TextMode(BW80);

  until Forever;

END.