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) $
8 ***********************************************************
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)
28 void qinit(REAL* table, INTEGER isize, INTEGER filter, REAL c)
44 // BRACEWELL BAND LIMITING FILTER
47 factr = c * c / (REAL) 4.0;
50 for (i = 1; i < isize; i++)
54 * ((REAL) 2.0 * sinc(mc) - SQR(sinc(mc / ((REAL) 2.0))));
58 // SHEPP LOW PASS FILTER
61 factr = SQR(c / Consts.pi);
62 table[0] = (REAL) 2.0 * factr;
65 for (i = 1; i < isize; i++)
68 amc = (REAL) 0.5 + mc;
69 smc = (REAL) 0.5 - mc;
71 * (sinc((REAL) 0.5 * amc) * (REAL) cos(Consts.pid2 * smc)
72 + sinc((REAL) 0.5 * smc)
73 * (REAL) cos(Consts.pid2 * amc));
77 // LOW PASS FILTER (COS)
80 factr = SQR(c / Consts.pi);
81 table[0] = factr * (Consts.pi - (REAL) 2.0);
82 factr *= Consts.pid2 * Consts.pid2;
84 for (i = 1; i < isize; i++)
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);
93 * (sinc(amc) + sinc(smc)
94 - (REAL) 0.5 * (SQR(h1) + SQR(h2)));
98 // GENERALIZED HAMMING
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;
107 for (k = 2; k < nsize; k += 2)
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));
114 table[isize - 1] = cm1 * factr / (REAL) 2.0
115 * ((REAL) pow((isize - 2), (-2)) + (REAL) pow((isize), (-2)));
119 fprintf(output, "\n **** filter = %d is not a valid filter type for qinit routine", filter);
120 fprintf(output, "\n **** program aborted\n");