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/back.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 COMBINED BACK-PROJECTION RECONSTRUCTION ALGORITHM.
12 THIS SUBROUTINE COMPUTES A BACK-PROJECTION RECONSTRUCTION
13 USING EITHER THE EQUATIONS OF THE LINES AS GIVEN BY RAY OR WRAY
14 OR BY INTERPOLATING BETWEEN THE GIVEN PROJECTION DATA TO
15 APPROXIMATE THE RADON TRANSFORM AT A PARTICULAR POINT. THE
16 FORMER IS SPECIFIED BY THE 'DISCRETE' PARAMETER, THE LATTER BY
19 IN THE DISCRETE CASE, THE USER HAS THE ADDITIONAL OPTIONS OF
20 SPECIFYING IF THE RAYSUM IS TO BE DIVIDED BY 'SNORM' BEFORE
21 BEING BACK-PROJECTED ('DISTRIBUTED') AND IF PIXELS WHICH ARE
22 INTERCEPTED BY A ZERO OR NEGATIVE RAY SHOULD BE REMOVED FROM
23 THE COMPUTATIONS ('ZEROIZING').
25 IN THE CONTINUOUS CASE, THE DESIRED INTERPOLATION METHOD MUST
26 BE SPECIFIED. WARNING- THE CONTINUOUS CASE WEIGHTS PROJECTION
27 "NP" BY "(THETA(NP+1)-THETA(NP-1))/2.0". IN ORDER TO ACCOMPLISH
28 THIS IT IS ASSUMED THAT THE PROJECTIONS ARE ARRANGED IN INCREASING
29 ORDER ON THETA AND THAT THETA(PRJNUM)-THETA(1) IS LESS THAN
30 PI IN THE PARALLEL CASE OR TWOPI IN THE DIVERGENT CASE.
32 THE RESULT OF THE BACK-PROJECTION WILL BE NORMALIZIED TO
33 THE CORRECT AVERAGE DENSITY IF EITHER 'ADDITIVE' OR 'MULTIPLICATIV
34 HAS BEEN SPECIFIED. THE FORMER PERFORMS THE NORMALIZATION BY
35 THE ADDITION OF A CONSTANT TO THE ENTIRE PICTURE, THE LATTER BY
59 BOOLEAN back_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
63 REAL* g = NULL; //wei 3/2005
64 BOOLEAN disc, zeros, dist, add, mult;
66 static const INTEGER hcont = CHAR2INT('c', 'o', 'n', 't');
67 static const INTEGER hdisc = CHAR2INT('d', 'i', 's', 'c');
68 static const INTEGER hdist = CHAR2INT('d', 'i', 's', 't');
69 static const INTEGER hzero = CHAR2INT('z', 'e', 'r', 'o');
70 static const INTEGER haddi = CHAR2INT('a', 'd', 'd', 'i');
71 static const INTEGER hmult = CHAR2INT('m', 'u', 'l', 't');
73 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
74 static const INTEGER back_codes[6] =
76 CHAR2INT('c', 'o', 'n', 't'),
77 CHAR2INT('d', 'i', 's', 'c'),
78 CHAR2INT('z', 'e', 'r', 'o'),
79 CHAR2INT('d', 'i', 's', 't'),
80 CHAR2INT('a', 'd', 'd', 'i'),
81 CHAR2INT('m', 'u', 'l', 't') };
109 neginf = -Consts.infin;
110 zflag = neginf + neginf;
116 // PROCESS CONTROL CARD
118 word = InFile.getwrd(TRUE, &eol, back_codes, 2);
123 fprintf(output, "\n **** neither continuous or discrete specified");
129 word = InFile.getwrd(FALSE, &eol, &(back_codes[2]), 4);
134 word = InFile.getwrd(FALSE, &eol, &(back_codes[3]), 3);
140 word = InFile.getwrd(FALSE, &eol, &(back_codes[4]), 2);
147 Fourie.interp = InFile.getint(FALSE, &eol);
148 if ((eol) || (Fourie.interp < -1) || (Fourie.interp > 6))
152 fprintf(output, "\n **** invalid interpolation method");
157 word = InFile.getwrd(FALSE, &eol, &(back_codes[4]), 2);
170 fprintf(output, "\n discrete back-projection");
175 fprintf(output, "\n continuous back-projection");
177 if (Fourie.interp == -1)
179 fprintf(output, "\n modified cubic spline interpolation");
182 if (Fourie.interp == 0)
184 fprintf(output, "\n band limiting sinc interpolation");
187 if (Fourie.interp > 0)
189 fprintf(output, "\n %2i order lagrange interpolation",
196 fprintf(output, "\n distributed back-projection");
201 fprintf(output, "\n zeroizing is in effect");
206 fprintf(output, "\n additive normalization");
211 fprintf(output, "\n multiplicative normalization");
219 nsize = GeoPar.snrays + 2 * Fourie.interp;
220 if (Fourie.interp == -1)
221 nsize = GeoPar.snrays + 6;
229 thmod = Consts.twopi;
232 Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
233 theta0 = theta - thmod;
235 for (np = 0; np < GeoPar.prjnum; np++)
238 ProjFile.ReadProj(np, g, nsize);
240 Anglst.getang(np, &theta1, &sinth, &costh);
244 if (np == (GeoPar.prjnum - 1))
247 Anglst.getang(np2, &theta2, &sinth2, &costh2);
249 if (np == (GeoPar.prjnum - 1))
252 w = (theta2 - theta0) / thfac;
257 "\n **** projection angles not in increasing order");
260 delete[] g; //wei 3/2005
264 Fourie.rinc = GeoPar.pinc;
267 Fourie.rinc *= GeoPar.radius / GeoPar.stod;
269 Fourie.rinc *= (REAL) MAX0(fabs(sinth), fabs(costh));
275 for (nr = 0; nr < nsize; nr++)
280 bckprj(recon, GeoPar.nelem, g, nsize, sinth, costh,
281 GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
290 g = new REAL[GeoPar.nrays];
295 // FLAG PIXELS INTERCEPTED WITH ZERO OR NEGATIVE RAYS
297 for (np = 0; np < GeoPar.prjnum; np++)
300 ProjFile.ReadProj(np, g, GeoPar.nrays);
302 for (nr = GeoPar.fsnray; nr <= GeoPar.lsnray; nr++)
303 { // changed < to <=. lajos Nov 17, 2004
304 if (g[nr] <= Consts.zero)
308 ray(np, nr, list, weight, &numb, &snorm);
310 wray(np, nr, list, weight, &numb, &snorm);
315 for (nb = 0; nb < numb; nb++)
325 // PERFORM DISCRETE BACK-PROJECTION
328 for (np = 0; np < GeoPar.prjnum; np++)
331 ProjFile.ReadProj(np, g, GeoPar.nrays);
332 for (nr = GeoPar.fsnray; nr <= GeoPar.lsnray; nr++)
333 { // changed < to <=. lajos Nov 17, 2004
336 ray(np, nr, list, weight, &numb, &snorm);
338 wray(np, nr, list, weight, &numb, &snorm);
348 // CORRECT SNORM BY REMOVING FLAGGED PIXELS
350 for (nb = 0; nb < numb; nb++)
353 if (recon[k] > neginf)
354 snorm += weight[nb] * weight[nb];
359 for (nb = 0; nb < numb; nb++)
362 recon[k] += amt * weight[nb];
370 // PERFORM NORMALIZATION
377 delete[] g; //wei 3/2005
380 // NO NORMALIZATION REQUESTED - JUST CLEAR FLAGGED PIXELS
382 for (k = 0; k < GeoPar.area; k++)
384 if (recon[k] <= neginf)
389 delete[] g; //wei 3/2005
392 // COMPUTE TOTAL DENSITY AND NUMBER OF UNFLAGGED PIXELS
398 for (k = 0; k < GeoPar.area; k++)
401 if (recon[k] > neginf)
412 // ADDITIVE CORRECTION TO AVERAGE DENSITY
414 corr = GeoPar.aveden - sum / ((REAL) (knt));
416 for (k = 0; k < GeoPar.area; k++)
419 if (recon[k] <= neginf)
423 delete[] g; //wei 3/2005
426 // MULTIPLICATIVE CORRECTION TO AVERAGE DENSITY
429 corr = ((REAL) (knt)) * GeoPar.aveden / sum;
431 for (k = 0; k < GeoPar.area; k++)
433 if (recon[k] <= neginf)
438 delete[] g; //wei 3/2005