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/conv.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 **********************************************************************
10 THIS IS A GENERAL IMPLEMENTATION OF THE THREE BEST KNOWN
11 CONVOLUTION METHODS, NAMELY THOSE OF BRACEWELL + RIDDLE[1],
12 RAMACHANDRAN + LAKSHMINARAYANAN[2], AND SHEPP + LOGAN[3].
13 THE CONTROL CARD FOLLOWING THE EXECUTE CARD CONTAINS (IN THE
14 SNARK FREE FORMAT) THE FOLLOWING INFORMATION:
15 1) THE DESIRED FILTER ("BAND", "SINC", "COSINE" OR "HAMMING")
16 2) THE CUTOFF PARAMETER (OR THE VALUE OF ALPHA FOR HAMMING)
17 3) THE INTERPOLATION METHOD DESIRED
18 THE BAND LIMITING FILTER CORRESPONDS TO THE FILTER USED BY
19 BRACEWELL + RIDDLE(1) AND RAMACHANDRAN + LAKSHMINARAYANAN(2).
20 THE SINC FILTER IS GENERALIZATION OF THE FILTER USED BY SHEPP
22 THE COSINE FILTER IS A COMMON LOWPASS FILTER APPLIED TO THE
23 RECONSTRUCTION PROBLEM.
24 THE HAMMING FILTER IS THE GENERALIZED HAMMMING FILTER, A SPECI
25 CASE OF WHICH WAS USED BY CHESLER + REIDERER(4).
26 THE CUTOFF VALUE REPRESENTS THE FRACTION OF FOURIER SPACE
27 POWER TO INCLUDE IN THE RECONSTRUCTION. WHEN 'CUTOFF' IS ZERO,
28 IT IS TAKEN TO 1.0. WHEN 'CUTOFF' IS NEGATIVE, IT IS SET
29 TO '2.0*PRJNUM*PINC/(PI*NELEM*PIXSIZ)'. IF 'CUTOFF' IS GREATER
30 THAN 1.0, IT IS SET TO 1.0.
31 IF 'UNI' = .TRUE. (UNIFORM SPACING SPECIFIED),
32 THEN THE CUTOFF FORMS A CIRCULAR REGION IN FOURIER SPACE
33 WRERE THE DIAMETER OF THE CIRCLE IS 'CUTOFF/PIXSIZ'.
34 IF 'VRI' = .TRUE. THEN THE CUTOFF FORMS A SQUARE REGION WHERE
35 THE LENGTH OF A SIDE IS 'CUTOFF/PIXSIZ'.
36 EXAMPLE CONTROL CARD SEQUENCE
38 BAND LIMITING CUTOFF = -1.0 INTERPOLATION = 2
40 THIS VERSION OF THE CONVOLUTION METHOD WEIGHTS PROJECTION "NP"
41 BY "(THETA(NP+1)-THETA(NP-1))/2.0". IN ORDER TO ACCOMPLISH THIS
42 IT IS ASSUMED THAT THE PROJECTIONS ARE ARRANGED IN INCREASING
43 ORDER ON THETA AND THAT THETA(PRJNUM)-THETA(1) IS LESS THAN
44 PI IN THE PARALLEL CASE OR TWOPI IN THE DIVERGENT CASE.
47 [1] BRACEWELL + RIDDLE, INVERSION OF FAN-BEAM SCANS IN RADIO
48 ASTROMOMY, ASTROPHYS. J., V.150, PP.427-434 (1967).
49 [2] RAMACHANDRAN + LAKSHMINARAYANAN, THREE-DIMENSIONAL
50 RECONSTRUCTION FROM RADIOGRAPHS AND ELECTRON MICROGRAPHS:
51 APPLICATION OF CONVOLUTIONS INSTEAD OF FOURIER TRANSFORMS,
52 PROC. NAT. ACAD. SCI. U.S., V.68, PP.2236-2240 (1971).
53 [3] SHEPP + LOGAN, THE FOURIER RECONSTRUCTION OF A HEAD SECTION,
54 IEEE TRANS. NUCL. SCI., NS-21, PP.21-43 (1974).
55 [4] CHESLER + RIEDERER, RIPPLE SUPPRESSION DURING RECONSTRUCTION
56 IN TRANSVERSE TOMOGRAPHY, PHYS. MED. BIOL., V.20, PP.632-636
79 BOOLEAN conv_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
84 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
85 static const INTEGER fname[4] =
87 CHAR2INT('b', 'a', 'n', 'd'),
88 CHAR2INT('s', 'i', 'n', 'c'),
89 CHAR2INT('c', 'o', 's', 'i'),
90 CHAR2INT('h', 'a', 'm', 'm')
99 REAL* g = NULL; //wei 3.2005
101 REAL* q = NULL; //wei 3.2005
102 REAL* gp = NULL; //wei 3.2005
131 thmod = Consts.twopi;
135 word = InFile.getwrd(TRUE, &eol, fname, 4);
136 for (filter = 0; filter < 4; filter++)
138 if (word == fname[filter])
142 fprintf(output, "\n***** %s is not a valid filter name", int2str(word)); //(JD 12/31/03)
146 L30: Fourie.cutoff = InFile.getnum(FALSE, &eol);
147 Fourie.interp = InFile.getint(FALSE, &eol);
151 fprintf(output, "\n *** not enough arguments");
156 if (Fourie.cutoff < 0.0)
157 Fourie.cutoff = (REAL) 2.0 * GeoPar.prjnum * GeoPar.pinc
158 / ((REAL) Consts.pi * GeoPar.nelem * GeoPar.pixsiz);
159 if (Fourie.cutoff <= Consts.zero)
160 Fourie.cutoff = 1.0; // bug 151 - swr - 9/24/05
161 if (Fourie.cutoff > 1.0)
164 fprintf(output, "\n convolution reconstruction algorithm");
165 fprintf(output, "\n using filter %s", int2str(word));
166 fprintf(output, "\n cutoff (or alpha) = %10.4f", Fourie.cutoff);
168 if ((Fourie.interp < -1) || (Fourie.interp > 6))
170 fprintf(output, "\n *** invalid interpolation %4i", Fourie.interp);
175 if (Fourie.interp > 0)
177 fprintf(output, "\n %2i point lagrange interpolation",
180 if (Fourie.interp < 0)
182 fprintf(output, "\n cubic spline interpolation");
184 if (Fourie.interp == 0)
186 fprintf(output, "\n band limiting (sinc) interpolation");
189 ndelta = 2 * Fourie.interp;
190 if (Fourie.interp < 0)
192 gpnum = GeoPar.snrays + ndelta;
194 // G IS THE BASE OF PROJECTION DATA
196 g = new REAL[GeoPar.usrays];
198 // Q IS THE BASE OF THE CONVOLUTING FUNCTION VALUES
200 qnum = GeoPar.nrays + ndelta;
204 qinit(q, qnum, filter, Fourie.cutoff);
206 // GP IS THE BASE OF THE CONVOLUTED PROJECTIONS
208 gp = new REAL[gpnum];
210 // BEGIN THE RECONSTRUCTION
212 Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
214 theta0 = theta - thmod;
215 for (np = 0; np < GeoPar.prjnum; np++)
218 ProjFile.ReadProj(np, g, GeoPar.usrays);
220 Anglst.getang(np, &theta1, &sinth, &costh);
222 // COMPUTE THE WEIGHT FOR THE PROJECTION
226 if (np == GeoPar.prjnum - 1)
229 Anglst.getang(np2, &theta2, &sinth2, &costh2);
231 if (np == (GeoPar.prjnum - 1))
233 w = (theta2 - theta0) / thfac;
236 fprintf(output, "\n *** angles not in increasing order");
243 delete[] gp; //wei 3/2005
248 Fourie.rinc = GeoPar.pinc;
251 Fourie.rinc = GeoPar.pinc * (REAL) MAX0(fabs(sinth), fabs(costh));
253 Fourie.rinc = GeoPar.pinc * GeoPar.radius / GeoPar.stod;
262 // COMPUTE CONVOLUTED PROJECTIONS
264 izro = (GeoPar.usrays - gpnum) / 2;
268 for (nr = 0; nr < gpnum; nr++)
274 iz = MIN0(izro, GeoPar.usrays);
277 for (i = 0; i < iz; i++)
287 if (iz <= GeoPar.usrays)
289 for (i = iz; i <= GeoPar.usrays; i++)
300 // BACK PROJECT THE DATA
301 bckprj(recon, GeoPar.nelem, gp, gpnum, sinth, costh,
302 GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
311 delete[] gp; //wei 3/2005