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.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 **********************************************************************
11 THIS IS A CONVOLUTION TYPE RECONSTRUCTION ALGORITHM FOR
12 DIVERGENT BEAM PROJECTION DATA. VARIOUS KINDS OF REGULARIZING
13 FUNCTIONS ARE IMPLEMENTED.
15 TWO (OR MORE) DATA CARDS HAVE TO BE SUPPLIED BY THE USER
16 FOLLOWING EXECUTION COMMAND. THE INFORMATION ON AND AFTER
17 THE THIRD CARD ARE FILTER PARAMETERS, THEREFORE THEY ARE
21 1 INTERP,MISSP,NCONR,WCONV,METHOD
23 3 C,ALPHA -- IF(NF.NE.4HLINE
25 3 X1,Y1,X2,Y2,...,XL,YL
27 DESCRIPTION OF INPUT PARAMETERS :
28 INTERP - PARAMETER WHICH SPECIFIES INTERPOLATION METHOD,
29 THE VALUE OF INTERP MUST LIE BETWEEN -1 AND 6 AND
30 HAS THE FOLLOWING MEANING:
31 -1 -- MODIFIED CUBIC SPLINE INTERPOLATION,
32 0 -- BAND-LIMITING (SINC) INTERPOLATION,
33 POSITIVE -- INTERP-POINTS LAGRANGE INTERPOLATION.
34 MISSP - EVERY MISSP'TH PROJECTION IS TO BE USED IN
36 NCONR - NUMBER OF RAYS USED ON EACH SIDE OF CONVOLUTION,
37 FULL CONVOLUTION USED IF NCONR IS NEGATIVE.
38 WCONV - WEIGHT FOR CENTER RAY IN PRE-SMOOTHING.
39 METHOD - PARAMETER WHICH DETERMINES THE METHOD IN FINDING
40 THE RAY NUMBER AND WEIGHT FOR BACK PROJECTION.
41 0 -- ACCURATE BACK PROJECTION,
42 1 -- WEIGHT FOR BACKPROJECTION APPROXIMATED,
43 -1 -- BOTH WIEGHT AND RAY NUMBER APPROXIMATED.
44 NF - FILTER NAME, CAN BE ONE OF LINE, HAMMING, SINC,
45 EXPONENTIAL, PARABOLIC, COSINE, SHEPP, OR BAND.
46 C - RATIO OF FILTER BANDWIDTH TO NYQUIST CUTOFF FREQ.
47 ALPHA - PARAMETER USED IN HAMMING FILTER ONLY.
48 (XN,YN)- COORDINATES WHICH DEFINE THE FILTER FUNCTION
49 SEE SNARK05 MANUAL FOR DETAIL
50 EXAMPLES OF CONTROL CARD SEQUENCE :
61 IT IS ASSUMED THAT PROJECTION ANGLES ARE GIVEN IN INCREASING
62 ORDER AND THE LAST ANGLE IS LESS THAN THE FIRST ANGLE PLUS
64 IT IS THE RESPONSBILITY OF THE USER TO CHECK THAT USRAYS
65 IS LARGE ENOUGH FOR THE INTERPOLATION METHOD TO WORK.
67 1) FAST IMAGE RECONSTRUCTION BASED ON A RADON INVERSION
68 FORMULA APPROPRIATE FOR RAPIDLY COLLECTED DATA
69 G. T. HERMAN AND A. NAPARSTEK
70 SIAM JOURNAL ON APPLIED MATHEMATICS V.33 (1977)
71 2) FILTER SELECTION FOR FAN-BEAM CONVOLUTION RECONSTRUCTION
73 T. CHANG AND G. T. HERMAN
74 J. COMPUTER ASSISTED TOMOGRAPHY, TO APPEAR
75 WRITTEN BY T. CHANG, MAY 1977
76 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
95 BOOLEAN dcon_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
97 //fprintf(output,"\n xalg4");
99 static const INTEGER band = CHAR2INT('b','a','n','d');
101 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
102 // added missing options from "dcon_dconft.cpp". Lajos, Jan 13, 2005
103 static const INTEGER dcon_codes[8] =
105 CHAR2INT('b','a','n','d'),
106 CHAR2INT('l','i','n','e'),
107 CHAR2INT('h','a','m','m'),
108 CHAR2INT('p','a','r','a'),
109 CHAR2INT('e','x','p','o'),
110 CHAR2INT('s','i','n','c'),
111 CHAR2INT('c','o','s','i'),
112 CHAR2INT('s','h','e','p')
117 REAL theta,sinth,costh;
248 if(GeoPar.nrays <= 1) {
250 fprintf(output, "\n ***** nrays too small - job aborted *****");
255 // ALL BUT THE FIRST ITERATION ARE FOR SMOOTHING ONLY
258 if(iter > 1) return false;
260 // CHECK DATA COLLECTING GEOMETRY
263 fprintf(output, "\n *****geometry error*****");
264 fprintf(output, "\n\n this algorithm is written for fan beam only \n\n");
268 // READ CHECK AND ECHO USER SUPPLIED INFORMATION
270 interp = InFile.getint(TRUE, &eol);
272 if((interp < -1) || (interp > 6)) {
274 fprintf(output, "\n ***** illegal value for interp *****");
278 i = GeoPar.nrays - GeoPar.snrays;
279 if(((interp == -1) && (i < 2)) || ((interp > 2) && (i < (interp-2)))) {
281 fprintf(output, "\n ***** size of modified projection data too small");
282 fprintf(output, "\n for interpolation method specified *****");
287 missp = InFile.getint(FALSE, &eol);
289 if((missp <= 0) || (missp >= GeoPar.prjnum)) {
290 fprintf(output, "\n ***** illegal value for missp *****");
296 nconr = InFile.getint(FALSE, &eol);
297 wconv = InFile.getnum(FALSE, &eol);
298 if((wconv <= 0) || (wconv > 1)) {
299 fprintf(output, "\n ***** illegal value for wconv *****");
305 method = InFile.getint(FALSE, &eol);
308 fprintf(output, "\n ***** insufficient data item - end of line encountered *****");
313 if(ret_xalg4) return ret_xalg4;
316 fprintf(output, "\n %1i point lagrange interpolation", interp);
320 fprintf(output, "\n cubic spline interpolation");
324 fprintf(output, "\n band limiting (sinc) interpolation");
327 fprintf(output, "\n one out of every %5i projection(s) used for reconstruction", missp);
330 if((nconr < 0) || (nconr >= GeoPar.nrays)) nconr = GeoPar.nrays - 1;
332 fprintf(output, "\n %3i points used on each side of convolution", nconr);
334 wconw = (REAL) 0.5 * ((REAL) 1.0 - wconv);
335 fprintf(output, "\n pre-smoothing weight on center ray = %6.4f", wconv);
336 fprintf(output, "\n pre-smoothing weight on side rays = %6.4f", wconw);
340 fprintf(output, "\n weight for back-projection approximated");
344 fprintf(output, "\n accurate back-projection");
348 fprintf(output, "\n both weight and ray number for back-projection approximated");
351 nf = InFile.getwrd(TRUE, &eol, dcon_codes, 8);
353 fprintf(output, "\n ***** insufficient data item - end of line encountered *****");
358 // USE HERMAN, LAKSHMINARAYANAN AND NAPARSTEK'S ALGORITHM FOR
359 // BAND-LIMITING FILTER
362 fprintf(output, "\n band-limiting filter");
364 dconbl(recon, interp, missp, nconr, wconv, method);
368 // INITIALIZE CONSTANTS
370 nraym = GeoPar.nrays - 1;
371 fmdray = (REAL) GeoPar.midray + 1;
372 mdraym = GeoPar.midray;
373 mraym2 = GeoPar.midray - 1;
374 halfl = (REAL) GeoPar.midpix * GeoPar.pixsiz;
375 rcube = GeoPar.radius * GeoPar.radius * GeoPar.radius;
376 alfa = GeoPar.pinc / GeoPar.stod;
378 // INITIALIZE CONSTANTS FOR TANGENT DETECTOR CASE
379 // DETECN IS THE BASE FOR THE ARRAY WHICH STORES THE
380 // ARC DETECTOR LOCATIONS IN TANGENT DETECTOR UNIT
382 if(!(GeoPar.arc || (GeoPar.nrays <= 3))) {
384 alfa = (REAL) atan((REAL) mdraym * alfa) / (REAL) mdraym;
387 detecn = new REAL[mraym2];
390 for(i = 0; i < mraym2; i++) {
391 detecn[i] = (REAL)( (REAL) sin(sigma) / (REAL) cos(sigma)) / pcovst;
402 // ALLOCATE WORK AREAS FOR CONVOLUTING FUNCTION
403 // Q1 IS THE INDEX FOR THE FIRST ELEMENT IN Q1(M)
404 // Q2 IS THE INDEX FOR THE FIRST ELEMENT IN Q2(M)
406 q1 = new REAL[GeoPar.nrays+1]; //(JD 11/06)
407 q2 = new REAL[GeoPar.nrays+1]; //(JD 11/06)
409 // CONVOLUTING FUNCTION
411 if(dconft(q1, q2, nconr+1, alfa, nf)) {
414 if(detecn != NULL) delete[] detecn ;
415 if(q1 != NULL) delete[] q1 ;
416 if(q2 != NULL) delete[] q2 ;
417 if(ncos != NULL) delete[] ncos ;
418 if(g != NULL) delete[] g ;
419 if(g1 != NULL) delete[] g1;
420 if(modi != NULL) delete[] modi ; //wei 3/2005
427 // WEIGHT ARRAY FOR CONVOLUTION
428 // NCOS IS THE BASE FOR CONVOLUTION WEIGHT ARRAY
430 ncos = new REAL[mdraym];
434 for(i = 0; i < mdraym; i++) {
435 ncos[i] = (REAL) cos(sigma);
439 // ALLOCATE AREAS FOR PROJECTION DATA
440 // G IS THE BASE FOR PROJECTION DATA FOR ONE PROJECTION
441 // G1 IS THE BASE FOR WEIGHTED PROJECTION DATA
442 // MODI IS THE BASE FOR MODIFIED PROJECTION DATA
444 g = new REAL[GeoPar.nrays];
445 g1 = new REAL[GeoPar.nrays];
447 modi = new REAL[GeoPar.nrays];
449 limmid = (REAL) MIN0(nconr, mdraym);
451 //mdmod = modi + GeoPar.midray + 1;
452 //mdg = g + GeoPar.midray + 1;
453 //mdg1 = g1 + GeoPar.midray + 1;
454 mdmod = modi + GeoPar.midray; //(JD 11/06)
455 mdg = g + GeoPar.midray; //(JD 11/06)
456 mdg1 = g1 + GeoPar.midray; //(JD 11/06)
458 anginy = SQR(GeoPar.pixsiz / GeoPar.radius);
459 angiin = (REAL) GeoPar.midpix * anginy;
460 sbtrin = (REAL) GeoPar.midpix * angiin;
461 factor = (REAL) 0.25 * alfa * GeoPar.radius / Consts.pisq;
462 //lpj = ((GeoPar.prjnum-1) / missp) * missp + 1;
463 lpj = ((GeoPar.prjnum-1) / missp) * missp; //new indexing (from zero) (JD)
465 Anglst.getang(lpj, &beta0, &sinth, &costh);
466 //Anglst.getang(1, &theta, &sinth, &costh);
467 Anglst.getang(0, &theta, &sinth, &costh); //new indexing (from zero) (JD)
469 if((theta + Consts.twopi) <= beta0) {
471 fprintf(output, "\n ***** illegal arrangement of projection angles *****");
475 if(detecn != NULL) delete[] detecn ;
476 if(q1 != NULL) delete[] q1 ;
477 if(q2 != NULL) delete[] q2 ;
478 if(ncos != NULL) delete[] ncos ;
479 if(g != NULL) delete[] g ;
480 if(g1 != NULL) delete[] g1;
481 if(modi != NULL) delete[] modi ; //wei 3/2005
486 // THIS IS WHERE THE RECONSTRUCTION STARTS
487 // BUILD UP RECONSTRUCTED PICTURE ONE PROJECTION AT A TIME
489 for(np = 0; np < GeoPar.prjnum; np += missp) {
490 ProjFile.ReadProj(np, g, GeoPar.nrays);
492 // CONVERT THE PROJECTION DATA TO ARC GEOMETRY BY LINEAR
493 // INTERPOLATION IF NECESSARY
495 if(!(GeoPar.arc || GeoPar.nrays <= 3)) {
499 //for(i = 0; i < mraym2; i++){
500 for(i = 1; i <= mraym2; i++){ //(JD 11/06)
502 *(mdg1+i) = qintp(fmdray + dtc, g, GeoPar.nrays, 2);
503 *(mdg1-i) = qintp(fmdray - dtc, g, GeoPar.nrays, 2);
506 for(i = 1; i < nraym; i++){
512 // PRE-SMOOTHING OF PROJECTION DATA
517 g[0] = wconv * datm + wconw * dath;
519 //for(i = 1; i < nraym; i++) {
520 for(i = 1; i < nraym; i++) { //(JD 11/06)
524 g[i] = wconv * datm + wconw * (dath + datl);
527 g[GeoPar.nrays-1] = wconv * dath + wconw * datm;
530 // WEIGHTED PROJECTION DATA
534 //for(i = 0; i < mdraym; i++) {
535 for(i = 1; i <= mdraym; i++) { //(JD 11/06)
536 //*(mdg1+i) = *(mdg+i) * ncos[i];
537 //*(mdg1-i) = *(mdg-i) * ncos[i];
538 *(mdg1+i) = *(mdg+i) * ncos[i-1]; //(JD 11/06)
539 *(mdg1-i) = *(mdg-i) * ncos[i-1]; //(JD 11/06)
545 mdmod[0] = *q1 * (*mdg1) + (*q2) * (*mdg);
548 //for(i = 0; i < limmid; i++) {
549 for(i = 1; i <= limmid; i++) { //(JD 11/06)
550 mdmod[0] += *(q1+i) * ( *(mdg1+i) + *(mdg1-i) )+ *(q2+i) * ( *(mdg+i) + *(mdg-i));
554 //for(k = 0; k < mdraym; k++) {
555 for(k = 1; k <= mdraym; k++) { //(JD 11/06)
566 //for(i = 0; i < nconr; i++) {
567 for(i = 1; i <= nconr; i++) { //(JD 11/06)
570 //for(i1 = i; i1 < nconr; i1++) {
571 for(i1 = i; i1 <= nconr; i1++) { //(JD 11/06)
573 //if(!((mr1-i1) <= g1)) {
574 if(!((mr1-i1) <= (g1-1))) { //i.e. if (mr1-i1) <g1
577 pr1 += *mq1 * *(mr1-i1);
578 pr2 += *mq2 * *(mr -i1);
579 pl1 += *mq1 * *(ml1+i1);
580 pl2 += *mq2 * *(ml +i1);
581 //break; //(JD 11/06)
588 pr1 += *mq1 * (*(mr1+i) + *(mr1-i));
589 pr2 += *mq2 * (*(mr +i) + *(mr -i));
590 pl1 += *mq1 * (*(ml1+i) + *(ml1-i));
591 pl2 += *mq2 * (*(ml +i) + *(ml -i));
596 //*(mdmod + k) = pr1 + pr2 * ncos[k];
597 //*(mdmod - k) = pl1 + pl2 * ncos[k];
598 *(mdmod + k) = pr1 + pr2 * ncos[k-1]; // (JD 11/07/03)
599 *(mdmod - k) = pl1 + pl2 * ncos[k-1]; // (JD 11/07/03)
604 fprintf(output,"\n xalg4: modi");
605 for(i = 0; i < GeoPar.nrays; i++) {
606 if((i % 3) == 0) fprintf(output,"\n");
607 fprintf(output," %36.30f", (double) modi[i]);
614 // COMPUTE THE TRAPEZOIDAL WIDTH FOR BETA INTEGRATION
619 if(np < lpj) Anglst.getang(np + missp, &theta, &sinth, &costh);
620 if((np < lpj) && (theta <= beta1)) {
621 fprintf(output, "\n ***** illegal arrangement of projection angles *****");
624 if(detecn != NULL) delete[] detecn ;
625 if(q1 != NULL) delete[] q1 ;
626 if(q2 != NULL) delete[] q2 ;
627 if(ncos != NULL) delete[] ncos ;
628 if(g != NULL) delete[] g ;
629 if(g1 != NULL) delete[] g1;
630 if(modi != NULL) delete[] modi ; //wei 3/2005
635 //if(np >= lpj) Anglst.getang(1,&theta, &sinth, &costh);
636 if(np >= lpj) Anglst.getang(0,&theta, &sinth, &costh); // (JD 11/07/03)
637 w = (REAL) 0.5 * (REAL) fmod(Consts.twopi + theta - beta0, Consts.twopi);
644 cbhl = halfl * cbeta;
645 sbhl = halfl * sbeta;
646 cbpxs = GeoPar.pixsiz * cbeta;
647 sbpxs = GeoPar.pixsiz * sbeta;
649 ucosi = GeoPar.radius - sbhl - cbhl;
653 wbkpjy = (GeoPar.radius + (REAL) 2.0 * (cbhl + sbhl)) * bovrc;
654 wincy = (REAL) -2.0 * cbpxs * bovrc;
655 wincx = (REAL) -2.0 * sbpxs * bovrc;
659 angcos = GeoPar.radius - cbhl;
660 anginc = cbpxs / GeoPar.radius - angiin;
661 subtr = cbhl / GeoPar.radius - sbtrin;
666 for(iy = 0; iy < GeoPar.nelem; iy++) {
674 wbkpjy = wbkpjy + wincy;
681 ang = (REAL) atan2(angsin, angcos) - subtr;
684 for(ix = 0; ix < GeoPar.nelem; ix++) {
690 ang = (REAL) atan2(usin, ucos);
695 wbkpj = b / (usin * usin + ucos * ucos);
699 if(method < 0) ang += anginc;
701 kpr = fmdray - ang / alfa;
702 recon[ind] += qintp(kpr, modi, GeoPar.nrays, interp) * wbkpj;
709 btime += bend - cend;
711 fprintf(output, "\n total time for obtaining convoluting function was %10.3f seconds", ftime); // changed precision to three digits - swr 1/21/06
712 fprintf(output, "\n total time for convolution was %10.3f seconds", ctime); // ditto
713 fprintf(output, "\n total time for back-projection was %10.3f seconds", btime); // ditto
715 if(detecn != NULL) delete[] detecn ;
716 if(q1 != NULL) delete[] q1 ;
717 if(q2 != NULL) delete[] q2 ;
718 if(ncos != NULL) delete[] ncos ;
719 if(g != NULL) delete[] g ;
720 if(g1 != NULL) delete[] g1;
721 if(modi != NULL) delete[] modi ; //wei 3/2005