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/rfl.cpp $
5 $LastChangedRevision: 94 $
6 $Date: 2014-07-02 18:53:30 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 RHO-FILTERED LAYERGRAM RECONSTRUCTION METHOD
11 THE CONTROL CARD FOLLOWING THE EXECUTE CARD CONTAINS (IN THE
12 SNARK FREE FORMAT) THE FOLLOWING INFORMATION:
13 1) THE DESIRED FILTER ("BAND", "SINC", "COSINE" OR "HAMMING")
14 2) THE CUTOFF PARAMETER (OR THE VALUE OF ALPHA FOR HAMMING)
15 3) THE INTERPOLATION METHOD DESIRED
16 4) BACK PROJECTION REGION SIZE (OPTIONAL)
17 THE FILTERS ARE CIRCULARLY SYMMETRIC VERSIONS OF THE
18 CORRESPONDING FILTERS USED IN THE CONVOLUTION ALGORITHM.
19 THE BAND LIMITING FILTER CORRESPONDS TO THE FILTER USED BY
20 BRACEWELL + RIDDLE[1] AND RAMACHANDRAN + LAKSHMINARAYANAN[2].
21 THE SINC FILTER IS A GENERALIZATION OF THE FILTER USED BY SHEPP
23 THE COSINE FILTER IS A COMMON LOWPASS FILTER APPLIED TO THE
24 RECONSTRUCTION PROBLEM.
25 THE HAMMING FILTER IS THE GENERALIZED HAMMING FILTER, A SPECIAL
26 CASE OF WHICH WAS USED BY CHESLER + RIEDERER[4].
27 THE CUTOFF VALUE REPRESENTS THE FRACTION OF THE FOURIER SPACE
28 POWER RELATIVE TO THE RAY SPACING TO INCLUDE IN THE RECONSTRUCTION
29 A CUTOFF AT THE FOLDING FREQUENCY OF THE 2D FOURIER TRANSFORM
30 CORRESPONDS TO CUTOFF=PINC/PIXSIZ. WHEN "CUTOFF" IS GIVEN AS 0.0,
31 IT IS SET TO 1.0. WHEN "CUTOFF" IS NEGATIVE, IT IS SET TO
32 "2.0*PRJNUM*PINC/(PI*NELEM*PIXSIZ)".
33 THE BACK PROJECTION REGION SIZE MAY BE USER SPECIFIED AS
34 ANY EVEN NUMBER GREATER THAN NELEM (WITHIN REASON), OR IF NOT
35 GIVEN DEFAULTS TO NELEM+1.
36 EXAMPLE CONTROL CARD SEQUENCE
38 BAND LIMITING CUTOFF = -1.0 INTERPOLATION = 2
40 THIS VERSION OF THE RFL 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.
46 [1] BRACEWELL + RIDDLE, INVERSION OF FAN-BEAM SCANS IN RADIO
47 ASTROMOMY, ASTROPHYS. J., V.150, PP.427-434 (1967).
48 [2] RAMACHANDRAN + LAKSHMINARAYANAN, THREE-DIMENSIONAL
49 RECONSTRUCTION FROM RADIOGRAPHS AND ELECTRON MICROGRAPHS:
50 APPLICATION OF CONVOLUTIONS INSTEAD OF FOURIER TRANSFORMS,
51 PROC. NAT. ACAD. SCI. U.S., V.68, PP.2236-2240 (1971).
52 [3] SHEPP + LOGAN, THE FOURIER RECONSTRUCTION OF A HEAD SECTION,
53 IEEE TRANS. NUCL. SCI., NS-21, PP.21-43 (1974).
54 [4] CHESLER + RIEDERER, RIPPLE SUPPRESSION DURING RECONSTRUCTION
55 IN TRANSVERSE TOMOGRAPHY, PHYS. MED. BIOL., V.20, PP.632-636
80 BOOLEAN rfl_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
82 REAL* reco = NULL; //wei, 3/2005
123 INTEGER n[3], filter;
125 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 16, 2004
126 static const INTEGER fname[4] =
128 CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'),
129 CHAR2INT('c', 'o', 's', 'i'), CHAR2INT('h', 'a', 'm', 'm') };
140 thmod = Consts.twopi;
143 word = InFile.getwrd(TRUE, &eol, fname, 4);
145 for (filter = 0; filter < 4; filter++)
147 if (word == fname[filter])
151 fprintf(output, "\n***** %4s is not a valid filter name", int2str(word));
154 L30: Fourie.cutoff = InFile.getnum(FALSE, &eol);
155 Fourie.interp = InFile.getint(FALSE, &eol);
159 fprintf(output, "\n *** not enough parameters");
163 melen = InFile.getint(FALSE, &eol);
164 if (Fourie.cutoff < 0.0)
165 Fourie.cutoff = (REAL) 2.0 * GeoPar.prjnum * GeoPar.pinc / ((REAL) Consts.pi * GeoPar.nelem * GeoPar.pixsiz);
166 if (Fourie.cutoff <= Consts.zero)
168 if (Fourie.cutoff > 1.0)
170 fprintf(output, "\n rho-filtered layergram reconstruction algorithm");
171 fprintf(output, "\n using filter %4s", int2str(word));
172 fprintf(output, "\n cutoff (or alpha) = %9.4f", Fourie.cutoff);
174 if ((Fourie.interp < -1) || (Fourie.interp > 6))
176 fprintf(output, "\n invalid interpolation %4i", Fourie.interp);
180 if (Fourie.interp > 0)
182 fprintf(output, "\n %2i point lagrange interpolation", Fourie.interp);
184 if (Fourie.interp < 0)
186 fprintf(output, "\n cubic spline interpolation");
188 if (Fourie.interp == 0)
190 fprintf(output, "\n band limiting (sinc) interpolation");
194 melen = GeoPar.nelem + 1;
197 if (melen < GeoPar.nelem)
199 fprintf(output, "\n***** melen is less than nelem, melen= %4i", melen);
203 fprintf(output, "\n backprojection is performed on a %5i x %5i region", melen, melen);
209 if ((2 * n[1]) != melen)
211 fprintf(output, "\n***** melen must be even, melen= %4i", melen);
215 narea = (melen + 2) * melen;
217 reco = new REAL[narea];
219 for (nn = 0; nn < narea; nn++)
224 // G IS THE BASE OF PROJECTION DATA
226 ndelta = 2 * Fourie.interp;
227 if (Fourie.interp < 0)
229 gnum = GeoPar.numray(melen, GeoPar.pixsiz) + ndelta;
233 // BEGIN THE RECONSTRUCTION
235 Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
237 theta0 = theta - thmod;
239 for (np = 0; np < GeoPar.prjnum; np++)
242 ProjFile.ReadProj(np, g, gnum);
244 Anglst.getang(np, &theta1, &sinth, &costh);
246 // COMPUTE THE WEIGHT FOR THE PROJECTION
250 if (np == (GeoPar.prjnum - 1))
253 Anglst.getang(np2, &theta2, &sinth2, &costh2);
255 if (np == (GeoPar.prjnum - 1))
258 w = (theta2 - theta0) / thfac;
262 fprintf(output, "\n *** angles are not in increasing order");
264 delete[] reco; //wei 3/2005
268 Fourie.rinc = GeoPar.pinc;
270 Fourie.rinc *= GeoPar.radius / GeoPar.stod;
272 Fourie.rinc = (REAL) (Fourie.rinc * MAX0(fabs(sinth), fabs(costh)));
277 for (nr = 0; nr < gnum; nr++)
281 // BACK PROJECT THE DATA
282 bckprj(reco, melen, g, gnum, sinth, costh, GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
286 // COMPUTE 2-D FOURIER TRANSFORM
287 delete[] g; // bug 92 - Lajos - 03/02/2005
291 // MULTIPLY BY FILTER
292 topinc = (REAL) (2.0 * GeoPar.pinc);
294 scale = topinc / (melen * GeoPar.pixsiz);
296 nrow = melen / 2 + 1;
300 for (i = 0; i < ncol; i++)
302 ip = MIN0(i, melen - i);
303 for (j = 0; j < nrow; j++)
307 r = (REAL) (sqrt(((REAL) (ip * ip + jp * jp))) * scale);
308 factr = qfilt(r, Fourie.cutoff, filter) / topinc;
311 reco[ind + 1] *= factr;
315 // INVERSE FOURIER TRANSFORM
318 // SCALE TO ESTIMATED AVERAGE DENSITY
320 offset = (melen - 1 - GeoPar.nelem) / 2;
323 for (i = 0; i < GeoPar.nelem; i++)
325 index = (i + offset) * melen + offset;
326 for (j = 0; j < GeoPar.nelem; j++)
328 recon[ind] = reco[index];
335 diff = GeoPar.aveden - sum / GeoPar.area;
337 for (i = 0; i < GeoPar.area; i++)
343 delete[] reco; //wei,3/2005