src/snark/qfilt.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/qfilt.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  COMPUTE THE VALUE OF THE FILTER FUNCTION FOR THE RHO-FILTERED
  11  LAYERGRAM ALGORITHM
  12  */
  13
  14 /*
  15  THE IMPLEMENTED FILTERS ARE-
  16  FILTER=0     BRACEWEL  ABS(R)*BOX(R/(2.0*CUTOFF))
  17  FILTER=1     SHEPP     ABS(R)*SINC(R/(2.0*CUTOFF))
  18  *BOX(R/(2.0*CUTOFF))
  19  FILTER=2     COSINE    ABS(R)*COS(PI*R/(2.0*CUTOFF))
  20  *BOX(R/(2.0*CUTOFF))
  21  FILTER=3     HAMMING   ABS(R)*(ALPHA+(1.0-ALPHA)*COS(PI*R))
  22  *BOX(R/2.0)
  23  WHERE ALPHA = CUTOFF
  24  */
  25
  26 #include <cstdlib>
  27 #include <cstdio>
  28 #include <cmath>
  29
  30 #include "blkdta.h"
  31 #include "qfilt.h"
  32 #include "consts.h"
  33 #include "uiod.h"
  34
  35 REAL qfilt(REAL r, REAL cutoff, INTEGER filter)
  36 {
  37         if (!((filter < 0) || (filter > 3)))
  38         {
  39                 switch (filter)
  40                 {
  41
  42                 case 0:
  43                         if (fabs(r) > cutoff)
  44                                 return 0.0;
  45                         return (REAL) fabs(r);
  46
  47                 case 1:
  48                         if (fabs(r) > cutoff)
  49                                 return 0.0;
  50                         return (REAL) (cutoff * sin(Consts.pid2 * fabs(r) / cutoff)
  51                                         / Consts.pid2);
  52
  53                 case 2:
  54                         if (fabs(r) > cutoff)
  55                                 return 0.0;
  56                         return (REAL) (fabs(r) * cos(Consts.pid2 * r / cutoff));
  57
  58                 case 3:
  59                         if (fabs(r) > 1.0)
  60                                 return 0.0;
  61                         return (REAL) (fabs(r)
  62                                         * (cutoff + (1.0 - cutoff) * cos(Consts.pi * r)));
  63                 }
  64         }
  65
  66         // INVALID FILTER, OUTPUT ERROR MESSAGE AND STOP
  67
  68         fprintf(output, "\n **** filter = %d is not a valid filter type for qfilt routine", filter);
  69         fprintf(output, "\n **** program aborted\n");
  70         exit(500);
  71 }