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/sirt.cpp $
5 $LastChangedRevision: 97 $
6 $Date: 2014-07-02 20:07:41 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 IMPLEMENTATION OF THE SIRT ALGORITHM OF PETER GILBERT
11 TOGETHER WITH EXTENSIONS AND MODIFICATIONS BY
12 LAKSHMINARAYANAN AND LENT OF SUNY BUFFALO
13 MAY, 1977; REVISED JANUARY, 1978.
14 SEE TECHNICAL REPORT 95 OF DEPT OF COMPUTER SCIENCE
15 PROGRAMMED IN ANSI FORTRAN FOR PORTABILITY
16 USES AUXILIARY FILES USER1 AND USER2 TO REDUCE CORE
18 METHOD...CONTROLS ALGORITHM SELECTION
19 IF 'GSIRT' USE LAK-LENT GENERALIZATION OF GILBERT SIRT
20 SEE EQUATION (20) OF REPORT 95
21 RHO(Q+1) = RHO(Q) + 1.0/SIGMA *(B-(DN)**-1 * PTRANS * P *
24 B = (DL)**-1 * P TRANS * R
25 DN = P TRANS * P * U, U A VECTOR OF ONES
26 DL IS A WEIGHTED SUM OF RAY LENGTHS
27 R IS VECTOR OF PROJECTION DATA
29 FOR 'LSIRT', MODIFIED SIRT ALGORITHMS OF LAK AND LENT
31 GENERAL FORM OF MODIFIED SIRT IS (SEE EQUATION (23) )
32 RHO(Q+1) = RHO(Q) + 1.0/SIGMA * (B - DV**-1 * PTRANS *
33 Z * V**2 * P * RHO(Q))
35 N(J) = NUMBER OF PIXELS IN RAY J
36 L(J) = LENGTH OF RAY J
38 FOR THE LSIRT ALGORITHMS, (Z * V**2)(J) TAKES THE FOLLOWING FORM
40 DV(I) TAKES THE FOLLOWING FORM
41 SUM OVER ALL RAYS J PASSING THROUGH PIXEL I OF N(J)**POWER1
42 *** DV IS DESIGNATED THROUGHOUT THIS PROGRAM BY THE SYMBOL DN ***
43 B (IN THE TECHNICAL REPORT) HAS THE FOLLOWING FORM
44 DV**-1 * PTRANS * Z * V * R
46 (Z * V * R)(J) = R(J)/(L(J)*W(J)) * N(J)**POWER1
47 VALUES OF LSIRT MODIFIER
48 LSIRT = 1 POWER1=1 POWER2=0
49 LSIRT = 2 POWER1 = 0 POWER2 = -1
50 LSIRT = 3 POWER1 = -1 POWER2 = -2
51 TO GET CONSTRAINED SIRT, USE THE SNARK MODE CARD
55 THE FOLLOWING SEQUENCE OF KEYWORDS IS AVAILABLE TO THE USER
58 THIS KEYWORD MUST BE PRESENT
59 IT MAY BE FOLLOWED BY EITHER 'GSIRT' OR 'LSIRT TYPE'.
60 TYPE IS AN INTEGER MODIFIER WHICH MUST TAKE ON VALUES
64 IF THIS IS PRESENT, IT MUST BE FOLLOWED BY THE FLOATING-
65 POINT VALUE OF THE RELAXATION PARAMETER;
67 IF THIS IS PRESENT, IT MUST BE FOLLOWED BY THE FLOATING-
69 RELAXATION FACTOR = 1.0/SIGMA
71 ***** RELAX AND SIGMA MUST NOT BOTH BE USED *****
73 IF NEITHER RELAX OR SIGMA IS USED, RELAXATION FACTOR =
77 MAY BE FOLLOWED BY 'BACKPROJECTED'
78 THIS CAUSES A BACKPROJECTED IMAGE TO BE USED AS THE
80 IF START IS NOT PRESENT, THE EXEC COMMAND OPTION IS USED
83 THIS CAUSES ALL RECONSTRUCTIONS (OTHER THAN THE INITIAL GUESS)
84 TO BE NORMALIZED TO AVEDEN. THIS STEP IS OMITTED IF
85 AVEDEN IS TOO SMALL OR IF THE TOTAL DENSITY OF THE
86 ESTIMATED IMAGE IS LOW
88 EXAMPLE INPUT SEQUENCES FOLLOW
90 METHOD GSIRT SIGMA = .5 START NORMAL
91 METHOD LSIRT 1 RELAX = 3.
93 **** THIS SIRT IS FOR STRIP PROJECTION DATA ONLY ****
95 ADDED FORMAT 1090 FOR WRITING INTO USER1 BAND USER2..7/21/88.SADA
113 #include "projfile.h"
119 #define nrdis(rinc) MIN0(((INTEGER)(size * (fabs(costh) + fabs(sinth))/(rinc) - 0.001)), limit)
121 INTEGER sirt_class::Init()
124 hmeth = CHAR2INT('m', 'e', 't', 'h');
125 hgsir = CHAR2INT('g', 's', 'i', 'r');
126 hlsir = CHAR2INT('l', 's', 'i', 'r');
127 hstar = CHAR2INT('s', 't', 'a', 'r');
128 hnorm = CHAR2INT('n', 'o', 'r', 'm');
129 hrela = CHAR2INT('r', 'e', 'l', 'a');
130 hsigm = CHAR2INT('s', 'i', 'g', 'm');
131 hback = CHAR2INT('b', 'a', 'c', 'k');
139 sbase = NULL; //wei, 3/2005
144 INTEGER sirt_class::Reset()
152 if (gen_inv_dn != NULL)
154 if (vector_b != NULL)
159 delete[] sbase; //wei, 3/2005
164 BOOLEAN sirt_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
193 static BOOLEAN normal;
194 BOOLEAN bkstrt, gsirt, sirt, eol;
196 INTEGER limit, first, last;
198 static INTEGER lsirt; // set static hstau
213 // STATEMENT FUNCTION ENABLING TO DISCARD RAYS THAT DONT
214 // INTERSECT THE RECONSTRUCTION REGION IN THE SUMMATIONS
215 // DESCRIBED IN THE PREFACE
220 //major changes in the input reading sequence, using new getwrd funcion. hstau 1/2004
222 //for stripe mode only!
225 fprintf(output, "\n *** this sirt algorithm is for strip data only, try quad algorithm for line data**");
229 // SELECT THE SNARK OR USER SET OF RAYS
230 limit = GeoPar.lusray - GeoPar.midray;
232 if (GeoPar.snrays != GeoPar.usrays)
235 // CHECK FOR SNARK OR USER
239 // SNARK OPTION IN EFFECT
241 fprintf(output, "\n snark rays selected");
243 limit = GeoPar.lsnray - GeoPar.midray;
247 // USER OPTION IN EFFECT
248 fprintf(output, "\n user rays selected");
255 word = InFile.getwrd(TRUE, &eol, exp0, 1);
259 fprintf(output, "\n **** error - keyword method is missing***");
265 word = InFile.getwrd(FALSE, &eol, exp1, 2);
269 fprintf(output, "\n **** error - one and only one of the following methods must be specified: gsirt or lsirt***");
273 gsirt = word == hgsir;
274 sirt = word == hlsir;
278 fprintf(output, "\n gsirt method");
282 lsirt = InFile.getint(FALSE, &eol);
285 fprintf(output, "\n **** error - must specify lsirt type ***");
288 if (!(lsirt > 0 && lsirt < 4))
292 fprintf(output, "\n lsirt method");
303 { hrela, hsigm, hstar, hnorm };
304 word = InFile.getwrd(FALSE, &eol, exp2, 4);
308 // finished reading input
312 relax = InFile.getnum(FALSE, &eol);
315 fprintf(output, "\n **** error - must specify relax ***");
318 if (relax < Consts.zero)
320 fprintf(output, "\n **** error - relax too small or negative ***");
323 sigma = (REAL) 1.0 / relax;
327 word = InFile.getwrd(FALSE, &eol, exp2a, 2);
330 else if (word == hnorm)
336 else if (word == hsigm)
338 sigma = InFile.getnum(FALSE, &eol);
341 fprintf(output, "\n **** error - must specify sigma ***");
344 if (sigma < Consts.zero)
346 fprintf(output, "\n **** error - sigma too small or negative ***");
351 word = InFile.getwrd(FALSE, &eol, exp2b, 2);
354 else if (word == hnorm)
361 else if (word == hstar)
366 word = InFile.getwrd(FALSE, &eol, exp3, 2);
367 if (eol || word == hback)
381 fprintf(output, "\n sigma = inverse of relaxation = %10.4f\n",
383 // removed from the echoing hstau.2003
385 // NONCONVERGENCE WARNING ON SIGMA
388 fprintf(output, "\n ***small sigma value. probable divergence.***"); //hstau
391 if (Modefl.lofl || Modefl.upfl)
393 fprintf(output, "\n this run is constrained sirt");
397 fprintf(output, "\n this run is unconstrained sirt");
402 fprintf(output, "\n this run is normalized sirt");
406 fprintf(output, "\n this run is unnormalized sirt");
409 // CHECK ON SMALL TOTAL DENSITY AND ISSUE WARNING
410 totden = ((REAL) (GeoPar.area)) * GeoPar.aveden;
412 if (fabs(totden) <= Consts.zero)
415 fprintf(output, "\n ***average density too small; no normalization***"); //hstau
420 fprintf(output, "\n starting picture is back-projected density \n");
424 fprintf(output, "\n starting picture is as in execute command");
427 // END INPUT CHECK AND ECHO
435 gen_inv_dn = new REAL[area];
436 vector_b = new REAL[area];
438 pbase = new REAL[GeoPar.nrays];
440 sbase = new INTEGER[GeoPar.prjnum];
442 size = GeoPar.nelem * GeoPar.pixsiz / (REAL) 2.0;
443 Fourie.rinc = GeoPar.pinc;
445 // GET DN AND B ON USER1
446 // GENERATE DN AND DL
447 // IF 'LSIRT' DONT GENERATE DL,
448 // GET P TRANSPOSE N*R/L DIRECTLY
449 // WRITTEN FOR VARIABLE RAY WIDTH
451 for (i = 0; i < area; i++)
457 for (np = 0; np < GeoPar.prjnum; np++)
459 Anglst.getang(np, &theta, &sinth, &costh);
462 ProjFile.ReadProj(np, pbase, GeoPar.nrays);
466 Fourie.rinc = GeoPar.pinc
467 * (REAL) MAX0(fabs(costh), fabs(sinth));
469 first = GeoPar.midray - nrdis(Fourie.rinc);
472 //last = GeoPar.nrays + 1 - first;
473 last = GeoPar.nrays - 1 - first; //should be 2 less. hstau
475 for (nr = first; nr <= last; nr++)
479 ray(np, nr, list, weight, &numb, &snorm);
483 Blob.bray(np, nr, list, weight, &numb, &snorm);
487 posit(np, nr, &ax, &ay, &cth, &sth);
489 * raylen(SOT_rect, ax, ay, cth, sth, 0., 0., size,
493 raynar = (REAL) (numb);
502 raynar = float(numb);
503 raylar = dl * raynar;
512 raynar = ((REAL) 1.0) / numb;
513 raylar = dl * raynar;
518 for (nb = 0; nb < numb; nb++)
533 // TRYING FOR COMPATIBILITY IN CASE SOME PIXELS DO NOT
534 // CONTRIBUTE TO THE PROJECTION DATA
536 for (i = 0; i < area; i++)
540 if (w2[i] > Consts.zero)
541 w2[i] = ((REAL) 1.0) / w2[i];
546 // W2 NOW HOLDS GENERALIZED INVERSE OF DN
549 for (i = 0; i < area; i++)
551 gen_inv_dn[i] = w2[i];
556 for (i = 0; i < area; i++)
561 for (np = 0; np < GeoPar.prjnum; np++)
563 ProjFile.ReadProj(np, pbase, GeoPar.nrays);
566 last = GeoPar.nrays - 1 - first; //should be 2 less. hstau
569 for (nr = first; nr <= last; nr++)
574 ray(np, nr, list, weight, &numb, &snorm);
578 Blob.bray(np, nr, list, weight, &numb, &snorm);
582 for (nb = 0; nb < numb; nb++)
592 for (i = 0; i < area; i++)
606 // WORK HOLDS PART OF B; GET GENERALIZED INVERSE IN w2
608 for (i = 0; i < area; i++)
610 if (w2[i] > Consts.zero)
611 w2[i] = ((REAL) 1.0) / w2[i];
614 for (i = 0; i < area; i++)
616 gen_inv_dn[i] = w2[i];
619 // SET UP B VECTOR FOR 'LSIRT'
621 for (i = 0; i < area; i++)
623 w2[i] = w1[i] * w2[i];
627 for (i = 0; i < area; i++)
634 for (i = 0; i < area; i++)
638 // ITERATIONS BEGIN HERE
641 for (i = 0; i < area; i++)
648 for (i = 0; i < area; i++)
650 w1[i] = recon[i] + w1[i] / sigma;
653 for (i = 0; i < area; i++)
658 // GOING TO GENERATE P TRANSPOSE P RHO IN WORK
660 for (i = 0; i < area; i++)
665 for (np = 0; np < GeoPar.prjnum; np++)
668 last = GeoPar.nrays - 1 - first; // should be 2 less. hstau
671 for (nr = first; nr <= last; nr++)
676 ray(np, nr, list, weight, &numb, &snorm);
680 Blob.bray(np, nr, list, weight, &numb, &snorm);
686 for (nb = 0; nb < numb; nb++)
695 psum /= float(numb * numb);
697 for (nb = 0; nb < numb; nb++)
699 w1[list[nb]] += psum;
705 for (i = 0; i < area; i++)
707 recon[i] = gen_inv_dn[i];
711 // RECON HOLDS GENERALIZED INVERSE OF DN
712 // WORK HOLDS PTRANSPOSE P RHO
714 for (i = 0; i < area; i++)
719 for (i = 0; i < area; i++)
725 for (i = 0; i < area; i++)
727 ronew = recon[i] - w1[i] / sigma;
729 ronew = MAX0(ronew, Modefl.lower);
731 ronew = MIN0(ronew, Modefl.upper);
736 if (i % 2 == Blob.pr)
753 if (fabs(psum) > Consts.zero)
755 ratio = totden / psum;
757 for (i = 0; i < area; i++)
766 "\n ***current total density too small; no normalization for this iteration***"); //hstau