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/art.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 THIS IS A GENERAL INPLEMENTATION OF ALGEBRAIC RECONSTRUUCTION
11 THIS SUBROUTINE IS CAPABLE OF EXECUTION THE FOLLOWING METHODS:
14 A RELAXATION METHOD FOR RECONSTRUCTING OBJECTS FROM NOISY
15 X-RAY, G.T. HERMAN, MATH. PROG.,V8(1975),1-19.
17 ITERATIVE QUADRATIC OPTIMIZATION WITH APPLICATION IN DI-
18 AGNOSTIC RADIOLOGY, G.T. HERMAN AND A. LENT, SUNYAB
19 DEPARTMENT OF COMPUTER SCIENCE, INTERNAL RE-
22 A COMPUTER IMPLEMENTATION OF A BAYESIAN ANALYSIS OF IMAGE
23 RECONSTRUCTION, G.T. HERMAN , A. LENT, H. HURWITZ, AND
24 H.P. LUNG, SUNYAB DEPARTMENT OF COMPUTER SCIENCE INTERNAL
26 WITH THE FOLLOWING OPTIONAL MODIFICATIONS ON THE RECONSTRUCTION
28 REFERENCE IS THE SAME AS ART4.
29 BINARY PATTERNS (BART)
30 RECONSTRUCTION OF BINARY PATTERNS FROM A FEW PROJECTIONS,
31 G.T. HERMAN, INTERNATIONAL COMPUTING SYMPOSIUM 1973,
35 KOUNT - NUMBER OF RAYS USED IN EACH ITERATION
36 HALFWY- EQUALS TO (UPPER+LOWER)/2.0
37 THTOL - EQUALS TO (UPPER-LOWER)/1000.0
38 THRLO - (1) STARTING FROM LOWER, AND IS INCREMENTED BY
39 (UPPER-LOWER)/(KOUNT*2**(ITER+1)) AFTER EACH RECON-
40 STRUCTION LOOP IN THIS ROUTINE WHEN BART IS
41 SPECIFIED (BART WILL BE EXPLAINED LATER)
42 (2) EQUALS TO LOWER OTHERWISE
43 THRUP - (1) STARTING FROM UPPER, AND IS DECREMENTED BY THE SAME
44 AMOUNT AS INDICATED IN THRLO (1) AFTER EACH RECON-
45 STRUCTION LOOP WHEN BART IS SPECIFIED
46 (2) EQUALS TO UPPER OTHERWISE
48 THE FOLLOWING KEYWORDS ARE AVAILABLE FOR THE USER TO IMPLEMENT
52 RECONSTRUCTION WILL BE BASED ON ART3 METHOD
54 RECONSTRUCTION WILL BE BASED ON ART4 METHOD
56 RECONSTRUCTION WILL BE BASED ON BAYESIAN METHOD
58 EITHER "ART2" OR "BOUND" OR "BART" MUST BE SPECIFIED AFTER
59 THIS KEYWORD. WHEN "ART2" IS SPECIFIED, THE ART2 WAY
60 CONSTRAINT IS USED IN THE RECONSTRUCTION. IF "BOUND" IS
61 SPECIFIED, THE RECONSTRUCTION IS BOUNDED BY THE BOUNDARIES
62 "UPPER" AND "LOWER". WHEN "BART" IS SPECIFIED, THE RECON-
63 STRUCTED PICTURE IS ASSUMED TO BE A BINARY PATTERN(WITH
64 GRAY LEVELS "UPPER" AND "LOWER"). AND THE PICTURE ELEMENTS AR
65 COMPARED WITH THRLO, THRUP, AND HALFWY WITH ASSOCIATED CON-
66 TRAINT. WHEN LOFL IS ON, IF THE ORIGINAL PICTURE ELEMENT IS
67 LESS THAN THRLO AND ITS RECONSTRUCTION IS GREATER THAN HALFWY,
68 THEN THIS ELEMENT IS CHANGED TO HALFWY-THTOL. WHEN UPFL IS
69 ON, IF THE ORIGINAL PICTURE ELEMENT IS GREATER THAN THRUP
70 AND ITS RECONSTRUCTION IS LESS THAN HALFWY THEN IT IS
71 CHANGED TO HALFWY+THTOL.
73 THE USER CAN SPECIFY THE VALUE OF THE RELAXATION PARAMETER TO
74 BE USED IN THE RECONSTRUCTION BY THIS OPTION. THE RELAX-
75 ATION PARAMETER CAN BE EITHER CONSTANT OR VARIABLE. IF A
76 CONSTANT IS DESIRED THE USER MUST SPECIFY "RELAX CONSTANT"
77 FOLLOWED BY A FLOATING POINT NUMBER. IF THE PARAMETER IS
78 VARIABLE THEN "RELAX VARIABLE" MUST BE SPECIFIED, AND A
79 FUNCTION RSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT)
80 FOR CALCULATING THE VARIABLE PARAMETER MUST
81 BE SUPPLIED BY THE USER, WHERE THE CALCULATED
82 RELAXATION PARAMETER IS RETURNED BY ASSIGNING THE VALUE TO
83 RSET. THE DEFAULT IS CONSTANT 1.0.
85 THE MODIFIERS FOR CONRELAX IS THE SAME AS OPTION "RELAX".
86 IN THE VARIABLE CASE THE NAME OF THE USER SUPPLIED
87 FUNCTION IS CRSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT).
89 AN INTEGER MUST BE GIVEN AFTER "NORM". THE USER HAS HIS
90 FREEDOM OF CHOOSING THE TYPE OF NORM (>0) USED IN THE ITER-
91 ATIVE STEP BY SPECIFYING THE INTEGER NUMBER OF THE NORM.
92 THE DEFAULT IS 2-NORM.
94 AN INTEGER MUST BE SPECIFIED AFTER THIS KEYWORD, AND THE
95 VALUE OF THE INTEGER IS THE NUMBER OF RAYS USED FOR THE RE-
96 CONSTRUCTION FOR EACH SNARK05 ITERATION. THE RAYS CHOSEN
97 ONE AFTER ANOTHER ARE DETERMINED BY THE SNARK05 "SELECT"
98 COMMAND. THE DEFAULT IS A NUMBER EQUIVALENT TO
101 THE FINAL RECONSTRUCTION OF EACH ITERATION
102 WILL BE NORMALIZED SO THAT THE AVERAGE DENSITY OF THE PICT-
103 URE EQUALS TO AVEDEN; OTHERWISE THE RECONSTRUCTION IS
106 THERE ARE THREE POSSIBILITIES TO DEFINE THE TOLERANCE FOR
108 (1) "TOLERANCE FIXED T", WHERE T IS A USER SUPPLIED FL-
109 OATING POINT NUMBER, AND IT IS USED AS THE TOLERANCE DUR-
110 ING THE RECONSTRUCTION.
111 (2) "TOLERANCE VARIABLE", THE TOLERANCE IS CHANGED FOR EACH
112 SNARK05 ITERATIVE STEP. THE USER HAS TO SUPPLY A FUNC-
113 TION TSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT)
114 SUCH THAT THE CALCULATED TOLERANCE IS ASSIGNED TO
116 (3) "TOLERANCE NOISE T", THE TOLERANCE IS DEPENDENT ON THE
117 NOISE IN THE PROJECTION DATA. IT IS CALCULATED BY
118 T*SQRT(VARIAN) WHERE VARIAN IS OBTAINED BY
119 CALLING FUNCTION ERRFAC.
141 INTEGER art_class::Init()
143 //Xalg1User3.Open(); //Ran, bug 266
144 //ArtIn.delta_1 = NULL; //wei, 3/2005 //Ran, bug 266
148 INTEGER art_class::Reset()
150 // Xalg1User3.Close(); //Ran, bug 266
151 // if (ArtIn.delta_1 != NULL) delete[] ArtIn.delta_1; //wei,3/2005 //Ran, bug 266
153 delete[] ArtIn.totalRays; //wei, bug 266
157 BOOLEAN art_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
160 static INTEGER iflagr;
161 static INTEGER iflagt;
162 static INTEGER iflagc;
189 static BOOLEAN flag; //may not need to be static hstau
213 readin(recon, list, weight, &erflag);
220 ArtIn.thrlo = Modefl.lower;
221 ArtIn.thrup = Modefl.upper;
223 ArtIn.halfwy = (Modefl.upper + Modefl.lower) / (REAL) 2.0;
224 ArtIn.thtol = (Modefl.upper - Modefl.lower) / (REAL) 1000.0;
225 dist = Modefl.upper - Modefl.lower;
230 Fourie.rinc = GeoPar.pinc;
231 rinc2 = GeoPar.pinc * GeoPar.pinc;
240 for (i = 0; i < iter + 1; i++)
243 thrinc = dist / (ArtIn.kount * pof2);
248 for (i = 0; i < ArtIn.kount; i++)
251 raysum = Anglst.prdta(np, nr);
256 // GET VARIABLE RELAXATION PARAMETER
258 ArtIn.relax = rset(recon, iter, np, nr, iflagr, raysum, list,
268 // GET VARIABLE TOLERANCE BASED ON NOISE IN PROJECTION
270 if (!(GeoPar.line || GeoPar.uni))
272 Anglst.getang(np, &theta, &sinth, &costh);
273 Fourie.rinc = GeoPar.pinc * MAX0(sinth, costh);
274 rinc2 = Fourie.rinc * Fourie.rinc;
276 v = errfac(raysum, Fourie.rinc, rinc2);
277 ArtIn.toler = ArtIn.ct * (REAL) sqrt(v);
284 // GET VARIABLE TOLERANCE
286 ArtIn.toler = tset(recon, iter, np, nr, iflagt, raysum, list,
292 if (ArtIn.toler < 0.0)
295 "\n ***tolerance calculated at iter= %5i np= %5i nr= %5i tolerance= %9.5f***",
296 iter, np, nr, ArtIn.toler); //hstau
299 "\n ***error - tolerance can not be less than 0.0***"); //hstau
310 ArtIn.cr = crset(recon, iter, np, nr, iflagc, raysum, list, weight,
320 psrsum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
321 GeoPar.line, ArtIn.art2f);
325 psrsum = Blob.bpseudo(recon, np, nr, list, weight, &numb, &snorm,
326 GeoPar.line, ArtIn.art2f);
332 // CALCULATE P-NORM WHEN P IS OTHER THAN 2
337 for (j = 0; j < numb; j++)
339 snorm += (REAL) pow(weight[j], ArtIn.p);
342 if (snorm > Consts.zero)
344 diff = raysum - psrsum;
345 abdiff = (REAL) fabs(diff);
348 bkproj(recon, diff, abdiff, ArtIn.totalRays, np, nr, list,
349 weight, numb, snorm);
353 // bug 222 - swr - 04/08/07
354 ArtIn.thrlo += thrinc;
355 ArtIn.thrup -= thrinc;
364 for (i = 0; i < area; i++)
367 temp = clip(temp, ArtIn.thrlo, ArtIn.thrup, 1.0);
370 if (i % 2 == Blob.pr)
381 averec = sum / float(area);
386 avedif = GeoPar.aveden - averec; //aveden in both basis coincide
388 // IN THE FOLLOWING LOOP OF 820 THE COMMENTED LINE WAS
389 // REPLACED WITH THE CODE INSIDE LOOP 820....10/13/87
390 // Dr. GABOR T HERMAN
392 for (i = 0; i < area; i++)
395 pictk = recon[i] + avedif;
398 pictk = clip(pictk, ArtIn.thrlo, ArtIn.thrup, 1.0);
401 if (Modefl.lofl && (recon[i] <= ArtIn.thrlo)
402 && (pictk >= ArtIn.halfwy))
403 pictk = ArtIn.halfwy - ArtIn.thtol;
405 if (Modefl.upfl && (recon[i] >= ArtIn.thrup)
406 && (pictk <= ArtIn.halfwy))
407 pictk = ArtIn.halfwy + ArtIn.thtol;