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/dcon_dconbl.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 **********************************************************************
10 CONVOLUTION-TYPE METHOD FOR DIVERGENT DATA WITH BAND-LIMITING
13 CONVOLUTION RECONSTRUCTION TECHNIQUES FOR DIVERGENT BEAMS
14 G. T. HERMAN, A. V. LAKSHMINARAYANAN AND A. NAPARSTEK
15 COMPUTERS IN BIOLOGY AND MEDICINE
17 WRITTEN BY G. T. HERMAN, MARCH 1975
18 REVISED FOR THE CDC3500, MARCH 1976
19 REVISED FOR SNARK77 BY T. CHANG, MAY 1977
20 ***********************************************************
47 for(int i = 1; i < 10; i++) {
55 a *= (2 * i + 1) * 2 * i; // !
56 z *= x * x; // x ^ (2 * i)
64 #define Sin(X) sin((X))
68 void dcon_class::dconbl(REAL* recon, INTEGER interp, INTEGER missp, INTEGER nconr, REAL wconv, INTEGER method)
70 ////////////////////////////////////////////////////////
151 ////////////////////////////////////////////////////////
154 // ASSIGN SPACE FOR SMOOTHED PROJECTION DATA, ECHO SMOOTHING WEIGHTS
156 g = new REAL[GeoPar.nrays];
157 limup = GeoPar.nrays - 1; // index of last ray
159 if(!((wconv > 0.0) && (wconv <= 1.0))) wconv = (REAL) 1.0;
160 wconw = ((REAL) 1.0 - wconv) / (REAL) 2.0;
162 // ALLOCATE SPACE FOR CONVOLUTION INFORMATION
164 j = new REAL[GeoPar.nrays];
166 mdrym1 = GeoPar.midray - 1;
168 //mq = new REAL[mdrym1];
169 mq = new REAL[mdrym1 + 1]; //(JD 1/20/04)
174 alfa = GeoPar.pinc / GeoPar.stod;
175 alfasq = alfa * alfa;
176 a = alfa * GeoPar.radius;
177 j[GeoPar.midray] = a;
179 for(i = 1; i <= GeoPar.midray; i++)
181 if(GeoPar.tang) j[GeoPar.midray + i] = a / (REAL) sqrt((REAL) 1.0 + alfasq * SQR(((REAL)(i))));
182 if(GeoPar.arc) j[GeoPar.midray + i] = a * (REAL) cos( ((REAL)(i)) * alfa);
183 j[GeoPar.midray - i] = j[GeoPar.midray + i];
186 ////////////////////////////////////////////////////
188 fprintf(output,"\n***** test data 1 *****************************");
189 for(i = 0; i < GeoPar.nrays; i++)
191 if((i % 3) == 0) fprintf(output,"\n");
192 fprintf(output," %36.30f", j[i]);
194 fprintf(output,"\n");
196 ////////////////////////////////////////////////////
199 // CONVOLUTING FUNCTION Q(M)
201 for(m = 0; m < GeoPar.midray; m++)
203 n = 2 * m + 1; // 1, 3, 5, ...
204 if(GeoPar.tang) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * n * a));
206 if(GeoPar.arc) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * (REAL) sin((REAL) (n * alfa))));
213 if(GeoPar.arc) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * yyy));
215 // fprintf(output,"\n %36.30f %36.30f %36.30f", (double) mq[m], (double) xxx, (double) alfa);
216 fprintf(output,"\n %36.30f", (double) mq[m]);
220 ////////////////////////////////////////////////////
222 fprintf(output,"\n***** test data 2*****************************");
223 for(m = 0; m < GeoPar.midray; m++)
225 if((m % 3) == 0) fprintf(output,"\n");
226 fprintf(output," %36.30f", mq[m]);
228 fprintf(output,"\n");
230 ////////////////////////////////////////////////////
234 // CONSTANTS NEEDED FOR CONVOLUTION AND BACKPROJECTION
236 nraym1 = GeoPar.nrays - 1;
237 if(GeoPar.tang) concen = (REAL) 1.0 / ((REAL) 8.0 * SQR(a));
238 if(GeoPar.arc) concen = (REAL) 1.0 / ((REAL) 8.0 * SQR(alfa));
239 limmid = MIN0(nconr, mdrym1);
240 mod = new REAL[GeoPar.nrays];
242 ///midmod = mod + midray;
243 halfl = ((REAL)(GeoPar.midpix)) * GeoPar.pixsiz;
244 twodr = (REAL) 2.0 / GeoPar.radius;
245 fmidry = (REAL) GeoPar.midray + 1;
247 // INITIALIZE BETA0 TO BE USED IN CALCULATING WEIGHTS FOR
250 lastpr = ((GeoPar.prjnum - 1) / missp) * missp;
251 Anglst.getang(lastpr, &beta0, &sinth, &costh);
252 Anglst.getang(0, &theta, &sinth, &costh);
263 // THIS IS WHERE THE RECONSTRUCTION STARTS
264 // BUILD UP RECONSTRUCTED PICTURE ONE PROJECTION AT A TIME
266 for(np = 0; np < GeoPar.prjnum; np += missp)
269 ProjFile.ReadProj(np, g, GeoPar.nrays);
271 // debug ////////////////////////////////////////////////////
272 // input verification
274 fprintf(output,"\n**********************************");
275 for(nr = 0; nr < GeoPar.nrays; nr++)
277 if((nr % 3) == 0) fprintf(output,"\n");
278 fprintf(output," %36.30f", g[nr]);
281 //////////////////////////////////////////////////////
286 g[0] = (wconv * dtam + wconw * dtar) * j[0];
288 for(nr = 1; nr < nraym1; nr++)
293 g[nr] = (wconv * dtam + wconw * (dtar + dtal)) * j[nr];
298 g[nraym1] = (wconv * dtam + wconw * dtal) * j[nraym1];
301 stime += send - bend;
305 fprintf(output,"\n***** test data 3 *****************************");
306 for(nr = 0; nr < nraym1; nr++) {
307 if((nr % 3) == 0) fprintf(output,"\n");
308 fprintf(output," %36.30f", g[nr]);
310 fprintf(output,"\n");
314 // COMPUTE CONVOLUTED PROJECTION DATA
316 p = concen * g[GeoPar.midray];
319 for(nm = 1; nm <= limmid; nm += 2) {
320 p += (g[GeoPar.midray + nm] + g[GeoPar.midray - nm]) * mq[mqm];
324 mod[GeoPar.midray] = p;
326 // WE HAVE NOW TAKEN CARE OF THE CENTER RAY
328 for(m = 1; m <= GeoPar.midray; m++) {
329 inr = GeoPar.midray + m;
330 pr = concen * g[inr];
331 inl = GeoPar.midray - m;
332 pl = concen * g[inl];
338 for(nm = 1; nm <= nconr; nm += 2) {
341 for(mm = nm; mm <= nconr; mm += 2) {
344 pr += g[inm] * mq[mqm];
346 pl += g[inp] * mq[mqm];
354 pr += ( g[inp] + g[inm] ) * mq[mqm];
357 pl += ( g[inp] + g[inm] ) * mq[mqm];
365 mod[GeoPar.midray + m] = pr;
366 mod[GeoPar.midray - m] = pl;
372 fprintf(output,"\n***** test data mod *****************************");
373 for(m = 0; m < GeoPar.nrays; m++) {
374 if((m % 3) == 0) fprintf(output,"\n");
375 fprintf(output," %36.30f", mod[m]);
377 fprintf(output,"\n");
380 ////////////////////////////////////////////////////
384 ctime += cend - send;
386 // BACK-PROJECT CONVOLUTED PROJECTION DATA
391 // CALCULATE WEIGHTS FOR BETA INTEGRATION
394 if(np >= lastpr) Anglst.getang(0, &theta, &sinth, &costh);
395 if(np < lastpr) Anglst.getang(np+missp, &theta, &sinth, &costh);
397 w = (REAL) fmod(Consts.twopi + theta - beta0, Consts.twopi) / (REAL) 2.0;
401 // CALCULATE BACK-PROJECTION INFORMATION
403 rcosi = halfl * (sbeta - cbeta);
404 rcosy = - GeoPar.pixsiz * sbeta;
405 rcosx = GeoPar.pixsiz * cbeta;
406 rsini = -halfl * (sbeta + cbeta);
409 if(GeoPar.arc) w /= SQR(GeoPar.radius);
411 wjrfi = w - wtwodr * rsini;
412 wjrfx = -wtwodr * rsinx;
413 wjrfy = -wtwodr * rsiny;
414 anginy = (SQR(GeoPar.pixsiz / GeoPar.radius)) / alfa;
415 angiin = ((REAL)(GeoPar.midpix)) * anginy;
416 fkjrfx = rcosx / a - angiin;
417 subtr = fkjrfx * GeoPar.midpix;
419 cob = GeoPar.radius - halfl * cbeta;
420 wrsq = w * SQR(GeoPar.radius);
421 prini = GeoPar.radius + rsini;
426 for(i = 0; i < GeoPar.nelem; i++)
446 if(GeoPar.arc) fract = (REAL) atan(fract);
447 fkjrf = fract / alfa - subtr;
451 for(k = 0; k < GeoPar.nelem; k++)
461 if(GeoPar.arc) usq += SQR(rcos);
465 if(GeoPar.arc) fract = (REAL) atan(fract);
466 fkjrf = fract / alfa;
471 if(method < 0) fkjrf += fkjrfx;
473 frynum = fmidry - fkjrf;
477 fprintf(output,"\n frynum = %36.30f wjrf = %36.30f", frynum, wjrf);
481 recon[ind] += wjrf * qintp(frynum, mod, GeoPar.nrays, interp);
489 btime += bend - cend;
492 fprintf(output, "\n total time for obtaining smoothed projection data was %9.3f", stime); // changed precision to two digits - swr 1/21/06
493 fprintf(output, "\n total time for convoluting projection data was %16.3f", ctime); // ditto
494 fprintf(output, "\n total time for back-projecting was %28.3f", btime); // ditto
497 if(j != NULL) delete[] j;
498 if(mq != NULL) delete[] mq ; //wei 3/2005
500 delete[] mod; // bug 92 - Lajos - 03/02/2005
501 delete[] g; // bug 92 - Lajos - 03/02/2005