src/snark/qinit.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/qinit.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  INITIALIZE FILTER FUNCTION TABLE
  11  FILTER = 1     FOR BAND LIMITING FILTER
  12  = 2     FOR LOW PASS SINC FILTER
  13  = 3     FOR LOW PASS COS  FILTER
  14  = 4     FOR GENERALIZED HAMMING FILTER
  15  (CUTOFF BECOMES PARAMETER ALPHA)
  16  */
  17
  18 #include <cstdlib>
  19 #include <cstdio>
  20 #include <cmath>
  21
  22 #include "blkdta.h"
  23 #include "sinc.h"
  24 #include "qinit.h"
  25 #include "consts.h"
  26 #include "uiod.h"
  27
  28 void qinit(REAL* table, INTEGER isize, INTEGER filter, REAL c)
  29 {
  30         REAL mc;
  31
  32         REAL factr;
  33         INTEGER i;
  34         REAL amc;
  35         REAL smc;
  36         REAL h1;
  37         REAL h2;
  38         REAL cm1;
  39         INTEGER nsize;
  40         INTEGER k;
  41
  42         switch (filter)
  43         {
  44         // BRACEWELL BAND LIMITING FILTER
  45
  46         case 0:
  47                 factr = c * c / (REAL) 4.0;
  48                 table[0] = factr;
  49
  50                 for (i = 1; i < isize; i++)
  51                 {
  52                         mc = i * c;
  53                         table[i] = factr
  54                                         * ((REAL) 2.0 * sinc(mc) - SQR(sinc(mc / ((REAL) 2.0))));
  55                 }
  56                 return;
  57
  58                 // SHEPP LOW PASS FILTER
  59
  60         case 1:
  61                 factr = SQR(c / Consts.pi);
  62                 table[0] = (REAL) 2.0 * factr;
  63                 factr *= Consts.pid2;
  64
  65                 for (i = 1; i < isize; i++)
  66                 {
  67                         mc = i * c;
  68                         amc = (REAL) 0.5 + mc;
  69                         smc = (REAL) 0.5 - mc;
  70                         table[i] = factr
  71                                         * (sinc((REAL) 0.5 * amc) * (REAL) cos(Consts.pid2 * smc)
  72                                                         + sinc((REAL) 0.5 * smc)
  73                                                                         * (REAL) cos(Consts.pid2 * amc));
  74                 }
  75                 return;
  76
  77                 // LOW PASS FILTER (COS)
  78
  79         case 2:
  80                 factr = SQR(c / Consts.pi);
  81                 table[0] = factr * (Consts.pi - (REAL) 2.0);
  82                 factr *= Consts.pid2 * Consts.pid2;
  83
  84                 for (i = 1; i < isize; i++)
  85                 {
  86                         mc = i * c;
  87                         amc = (REAL) 0.5 + mc;
  88                         smc = (REAL) 0.5 - mc;
  89                         h1 = sinc((REAL) 0.5 * amc);
  90                         h2 = sinc((REAL) 0.5 * smc);
  91                         table[i] =
  92                                         factr
  93                                                         * (sinc(amc) + sinc(smc)
  94                                                                         - (REAL) 0.5 * (SQR(h1) + SQR(h2)));
  95                 }
  96                 return;
  97
  98                 // GENERALIZED HAMMING
  99
 100         case 3:
 101                 cm1 = c - (REAL) 1.0;
 102                 factr = (REAL) 1.0 / (Consts.pi * Consts.pi);
 103                 table[0] = c / (REAL) 4.0 + cm1 * factr;
 104                 table[1] = -c * factr - cm1 / (REAL) 8.0;
 105                 nsize = isize - 1;
 106
 107                 for (k = 2; k < nsize; k += 2)
 108                 {
 109                         table[k] = cm1 * factr / (REAL) 2.0
 110                                         * ((REAL) pow((k - 1), (-2)) + (REAL) pow((k + 1), (-2)));
 111                         table[k + 1] = -c * factr * (REAL) pow((k + 1), (-2));
 112                 }
 113
 114                 table[isize - 1] = cm1 * factr / (REAL) 2.0
 115                                 * ((REAL) pow((isize - 2), (-2)) + (REAL) pow((isize), (-2)));
 116                 return;
 117
 118         default:
 119                 fprintf(output, "\n **** filter = %d is not a valid filter type for qinit routine", filter);
 120                 fprintf(output, "\n **** program aborted\n");
 121                 exit(500);
 122                 //}
 123         }
 124 }