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/exalg.cpp $
5 $LastChangedRevision: 164 $
6 $Date: 2015-05-14 00:49:54 -0400 (Thu, 14 May 2015) $
8 ***********************************************************
10 EXALG DETERMINES WHICH ALGORITHM IS TO BE EXECUTED, WHICH OF
11 THE STANDARD OPTIONS ARE IN EFFECT AND REPEATEDLY EXECUTES THE
12 ALGORITHM UNTIL THE TERMINATION CRITERION (SET BY STOP COMMAND)
70 // macro to compute the number of elements in an array
71 #define length(array) (sizeof(array)/sizeof(array[0]))
73 static int totalIters = 0;
78 INTEGER Length;INTEGER* Names;
80 INTEGER GetIndex(INTEGER word)
82 for (int i = 0; i < Length; i++)
92 INTEGER GetName(INTEGER index)
98 // Instantiate Algorithm classes
123 // Algorithms Classes
124 alg_class* AlgClasses[] =
126 &alb1, // (ALB1) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
127 &alb2, // (ALB2) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
128 &alb3, // (ALB3) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
129 &alb4, // (ALB4) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
130 &alb5, // (ALB5) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
131 &alp1, // (ALP1) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
132 &alp2, // (ALP2) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
133 &alp3, // (ALP3) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
134 &alp4, // (ALP4) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
135 &alp5, // (ALP5) DUMMY RECONSTRUCTION ALGORITHM - TO BE REPLACED BY USER
136 &emap, // (EMAP) maximum a-posteriori probability expectation maximization algorithm
137 &lino, // (LINO) LINOGRAM RECONSTRUCTION ALGORITHM
138 &art, // (ART) ALGEBRAIC RECONSTRUUCTION
139 &sart, // (SART) SIMULTANEOUS ALGEBRAIC RECONSTRUUCTION TECHNIQUE
140 &conv, // (CONV) GENERAL IMPLEMENTATION OF THE THREE BEST KNOWN CONVOLUTION METHODS
141 &sirt, // (SIRT) IMPLEMENTATION OF THE SIRT ALGORITHM
142 &dcon, // (DCON) A CONVOLUTION TYPE RECONSTRUCTION ALGORITHM FOR DIVERGENT
143 // BEAM PROJECTION DATA
144 &foru, // (FOUR) FOURIER RECONSTRUCTION METHOD
145 &rfl, // (RFL) RHO-FILTERED LAYERGRAM RECONSTRUCTION METHOD
146 &back, // (BACK) COMBINED BACK-PROJECTION RECONSTRUCTION ALGORITHM
147 &quad, // (QUAD) QUADRATIC ITERATIVE RECONSTRUCTION ALGORITHMS
148 &mart, // (MART) VARIOUS VERSIONS OF MULTIPLICATIVE ART ALGORITHM
149 &bdhk // (BDHK) ALGORITHM BASED ON BUTNARIU-DAVIDI-HERMAN-KAZANTSEV
150 // MINIMIZING TOTAL-VARIATION WITH PERTURBATIONS
153 #define NUMBER_OF_ALGORITHMS length(AlgClasses)
156 static const INTEGER algn[NUMBER_OF_ALGORITHMS] =
158 CHAR2INT('a', 'l', 'b', '1'),
159 CHAR2INT('a', 'l', 'b', '2'),
160 CHAR2INT('a', 'l', 'b', '3'),
161 CHAR2INT('a', 'l', 'b', '4'),
162 CHAR2INT('a', 'l', 'b', '5'),
163 CHAR2INT('a', 'l', 'p', '1'),
164 CHAR2INT('a', 'l', 'p', '2'),
165 CHAR2INT('a', 'l', 'p', '3'),
166 CHAR2INT('a', 'l', 'p', '4'),
167 CHAR2INT('a', 'l', 'p', '5'),
168 CHAR2INT('e', 'm', 'a', 'p'),
169 CHAR2INT('l', 'i', 'n', 'o'),
170 CHAR2INT('a', 'r', 't', ' '),
171 CHAR2INT('s', 'a', 'r', 't'),
172 CHAR2INT('c', 'o', 'n', 'v'),
173 CHAR2INT('s', 'i', 'r', 't'),
174 CHAR2INT('d', 'c', 'o', 'n'),
175 CHAR2INT('f', 'o', 'u', 'r'),
176 CHAR2INT('r', 'f', 'l', ' '),
177 CHAR2INT('b', 'a', 'c', 'k'),
178 CHAR2INT('q', 'u', 'a', 'd'),
179 CHAR2INT('m', 'a', 'r', 't'),
180 CHAR2INT('b', 'd', 'h', 'k')
183 static NameTable Algorithms = { NUMBER_OF_ALGORITHMS, (int*) algn };
185 // Blob Algorithms Names
186 static const INTEGER algn_blob[] =
188 CHAR2INT('a', 'l', 'b', '1'),
189 CHAR2INT('a', 'l', 'b', '2'),
190 CHAR2INT('a', 'l', 'b', '3'),
191 CHAR2INT('a', 'l', 'b', '4'),
192 CHAR2INT('a', 'l', 'b', '5'),
193 CHAR2INT('e', 'm', 'a', 'p'),
194 CHAR2INT('a', 'r', 't', ' '),
195 CHAR2INT('s', 'a', 'r', 't'),
196 CHAR2INT('s', 'i', 'r', 't'),
197 CHAR2INT('q', 'u', 'a', 'd'),
198 CHAR2INT('m', 'a', 'r', 't')
201 #define NUM_OF_ALG_BLOB length(algn_blob)
203 static NameTable BlobAlgorithms = { NUM_OF_ALG_BLOB, (int*) algn_blob };
207 static const INTEGER exec_codes[30] =
209 //do not modify the first 4 entries!
210 CHAR2INT('a', 'v', 'e', 'r'),
211 CHAR2INT('z', 'e', 'r', 'o'),
212 CHAR2INT('c', 'o', 'n', 't'),
213 CHAR2INT('p', 'h', 'a', 'n'),
214 //add any additional options to EXECUTE (first option) here
216 CHAR2INT('a', 'l', 'b', '1'),
217 CHAR2INT('a', 'l', 'b', '2'),
218 CHAR2INT('a', 'l', 'b', '3'),
219 CHAR2INT('a', 'l', 'b', '4'),
220 CHAR2INT('a', 'l', 'b', '5'),
221 CHAR2INT('a', 'l', 'p', '1'),
222 CHAR2INT('a', 'l', 'p', '2'),
223 CHAR2INT('a', 'l', 'p', '3'),
224 CHAR2INT('a', 'l', 'p', '4'),
225 CHAR2INT('a', 'l', 'p', '5'),
226 CHAR2INT('e', 'm', 'a', 'p'),
227 CHAR2INT('l', 'i', 'n', 'o'),
228 CHAR2INT('a', 'r', 't', ' '),
229 CHAR2INT('s', 'a', 'r', 't'),
230 CHAR2INT('c', 'o', 'n', 'v'),
231 CHAR2INT('s', 'i', 'r', 't'),
232 CHAR2INT('d', 'c', 'o', 'n'),
233 CHAR2INT('f', 'o', 'u', 'r'),
234 CHAR2INT('r', 'f', 'l', ' '),
235 CHAR2INT('b', 'a', 'c', 'k'),
236 CHAR2INT('q', 'u', 'a', 'd'),
237 CHAR2INT('m', 'a', 'r', 't'),
238 CHAR2INT('b', 'd', 'h', 'k'),
239 //add any additional reconstruction algorithms here
241 CHAR2INT('s', 'm', 'o', 'o'),
242 CHAR2INT('c', 'o', 'n', 't') };
244 //these are variables indicating indexes into the exec_codes array above
245 //the name ending in "St" indicates a starting index
246 //the name ending in "Cnt" indicates number of entries in given category
247 // !!! if you keep these indexes up to date after making changes to the
248 // !!! array above, there is not need for changing indexes in the rest of
249 // !!! the code below
250 //there are four EXECUTE options starting at index 0 (these are the
251 //first options - i.e. before the name of algorithm)
252 INTEGER executeOp1St = 0;
253 INTEGER executeOp1Cnt = 4;
254 //there are 24 algorithm names starting at index 4
257 //there are 2 EXECUTE options starting at index 28 (these are the
258 //second options - i.e. after the name of algorithm)
259 INTEGER executeOp2St = 27;
260 INTEGER executeOp2Cnt = 2;
264 BOOLEAN flagphan = FALSE;
288 static BOOLEAN pix_old_flag = FALSE; // basis of the previous run
289 static REAL delta_old; // delta of the previous blob run
290 static REAL support_old; // support of the previous blob run
291 static REAL shape_old; // shape of the previous blob run
293 static BOOLEAN called = FALSE;
294 static INTEGER* list;
299 INTEGER area; //size of picture
300 INTEGER wx; //size of picture in x direct
301 INTEGER wy; //size of picture in y direct
303 if (ProjFile.Open("prjfil") != 0)
305 fprintf(output, "\n **** unable to open prjfil");
306 fprintf(output, "\n **** program aborted\n");
310 ProjFile.ReadHeader(head, &NoisePar, &Spctrm, &Anglst);
312 // PERFORM INITIALIZATION ONLY ON THE FIRST CALL
316 Creacm.recon = new REAL[GeoPar.area];
319 Blob.recon = new REAL[Blob.area];
322 recon_copy = new REAL[GeoPar.area];
328 RaySel.numray = GeoPar.snrays;
332 RaySel.numray = GeoPar.usrays;
337 RaySel.initnr = GeoPar.fsnray;
338 RaySel.endray = GeoPar.lsnray;
342 RaySel.initnr = GeoPar.fusray;
343 RaySel.endray = GeoPar.lusray;
346 RaySel.begray = RaySel.initnr;
347 RaySel.curnr = RaySel.initnr;
349 if (RaySel.stype == 0)
351 RaySel.curnr = RaySel.initnr - RaySel.rayinc;
355 RaySel.curnp = RaySel.initnp;
357 if (RaySel.stype == 1)
359 RaySel.curnp = RaySel.initnp - RaySel.prjinc;
365 // INITIALIZE OPTION FLAGS
367 for (i = 0; i < 50; i++)
374 word = InFile.getwrd(FALSE, &eol, exec_codes, executeOp1Cnt + algCnt);
379 if (word == exec_codes[0])
382 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
384 else if (word == exec_codes[1])
386 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
388 else if (word == exec_codes[2])
391 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
393 else if (word == exec_codes[3])
396 if (Creacm.test == NULL)
398 fprintf(output, "\n **** error: in conjuction with the option phantom, picture command's option must be test");
399 fprintf(output, "\n **** algorithm execution aborted");
402 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
405 // LOOK UP ROUTINE NAME
406 // check algorithm name
407 if ((AlgIndex = Algorithms.GetIndex(word)) < 0)
409 fprintf(output, "\n **** unknown algorithm %s", int2str(word));
410 fprintf(output, "\n **** algorithm execution aborted");
417 BOOLEAN pix_basis = Blob.pix_basis;
419 //check for algorithm that goes with blobs
422 // check if selected algorithm is iterative
423 if (BlobAlgorithms.GetIndex(AlgName) < 0)
425 fprintf(output, "\n **** blob basis is not appropriate for algorithm %s, the algorithm is being executed using pixel basis instead.", int2str(AlgName));
428 Blob.pix_basis = true;
432 // ACCEPT AND DISPLAY COMMENT CARD
433 InFile.GetComment(head);
435 fprintf(output, "\n %s", head);
437 word = InFile.getwrd(FALSE, &eol, &(exec_codes[executeOp2St]),
440 // CHECK TO SEE IF OPTION HAS BEEN SPECIFIED
443 if (word == exec_codes[executeOp2St])
449 if (word == exec_codes[executeOp2St + 1])
455 thresh = InFile.getnum(TRUE, &eol);
456 w1 = InFile.getnum(FALSE, &eol);
457 w2 = InFile.getnum(FALSE, &eol);
458 w3 = InFile.getnum(FALSE, &eol);
460 // ECHO OPTION SELECTED
463 fprintf(output, "\n reconstruction is smoothed after");
468 fprintf(output, "\n reconstruction is contoured after");
471 InFile.listit(optflg);
477 fprintf(output, "\n threshold = %9.5f w1 = %9.5f w2 = %9.5f w3 = %9.5f", thresh, w1, w2, w3);
481 fprintf(output, "\n threshold = %9.5f w1 = %9.5f w2 = %9.5f w3 = %9.5f", thresh, w1, w2, w3);
488 fprintf(output, "\n threshold = %9.5f w1 = %9.5f w2 = %9.5f w3 = %9.5f", thresh, w1, w2, w3);
494 if (!pix_old_flag && !pix_basis)
495 { // previously blob, now blob
496 if (fabs(support_old - Blob.support) > Consts.zero || fabs(shape_old - Blob.shape) > Consts.zero || fabs(delta_old - Blob.delta) > Consts.zero)
498 fprintf(output, "\n **** error: current blob parameters must match those of the previous blob basis execution");
499 fprintf(output, "\n **** algorithm execution aborted");
505 if (!pix_old_flag && pix_basis)
506 { // previously blob, now pixel
507 Blob.blob2pix(Blob.recon, Creacm.recon);
511 if (pix_old_flag && !pix_basis)
512 { // previously pixel, now blob
513 Blob.recon = new REAL[Blob.area];
514 Blob.pix2blob(Creacm.recon, Blob.recon);
520 Blob.recon = new REAL[Blob.area];
522 // initiate recon with val or phantom image
523 if (flagphan == TRUE)
525 for (i = 0; i < GeoPar.area; i++)
527 Creacm.recon[i] = Creacm.test[i];
531 Blob.pix2blob(Creacm.recon, Blob.recon);
538 for (i = 0; i < GeoPar.area; i++)
540 Creacm.recon[i] = val;
545 for (i = 0; i < Blob.area; i++)
555 bldlst(&list, &weight);
556 recon = Creacm.recon;
563 Blob.bbldlst(&list, &weight);
575 // initialize selected algorithm
576 if (AlgClasses[AlgIndex]->Init() != 0)
578 fprintf(output, "\n **** algorithm %s failed to initialize", int2str(algn[AlgIndex]));
579 fprintf(output, "\n **** algorithm execution aborted");
586 // run until the termination criterion is met or an error occurs
591 iter = ((count - 1) % 50) + 1;
593 // execute single iteration of the selected algorithm
594 if (SuperSet.superiorizationEnabled)
596 tflag = executeSuperiorization(recon, list, weight, count, AlgClasses[AlgIndex]);
600 tflag = AlgClasses[AlgIndex]->Run(recon, list, weight, count);
604 fprintf(output, "\n algorithm executed in iteration %4i", count);
605 time = end_alg - begin;
606 fprintf(output, "\n %10.3f seconds for the execution of the algorithm", time);
611 fprintf(output, "\n **** algorithm indicated failure to continue");
613 // reset algorithm after failure
614 if (AlgClasses[AlgIndex]->Reset() != 0)
616 fprintf(output, "\n **** algorithm %s failed to reset", int2str(algn[AlgIndex]));
619 fprintf(output, "\n **** algorithm execution aborted");
623 if (Term.termInstance == NULL)
625 tflag = ((count - Term.iterations) >= 0);
629 tflag = Term.termInstance->Run(recon, list, weight, count);
634 // END THE ROUTINE IF THE TEST WAS SATISFIED
635 // CHECK IF OPTION SPECIFIED FOR LAST ITERATION
640 { // if currently blob, pixelize it
643 Blob.blob2pix(Blob.recon, Creacm.recon);
646 if ((optsw == 1) && (optflg[iter] != 0))
648 smooth(Creacm.recon, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
651 ////////// dump recon ////////////
652 fprintf(output, "\nexalg: algorithm %s result of iteration %i after smooth", int2str(algn[AlgIndex]), count);
653 for (i = 0; i < area; i++)
655 if ((i % 3) == 0) fprintf(output, "\n");
656 fprintf(output, " %36.30f", (double) recon[i]);
658 //////////////////////////////////
662 if ((optsw == 2) && (optflg[iter] != 0))
664 contur(Creacm.recon, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
667 ////////// dump recon ////////////
668 fprintf(output, "\nexalg: algorithm %s result of iteration %i after contur", int2str(algn[AlgIndex]), count);
669 for (i = 0; i < area; i++)
671 if ((i % 3) == 0) fprintf(output, "\n");
672 fprintf(output, " %36.30f", (double) recon[i]);
674 //////////////////////////////////
678 RecFile.WriteRec(head, int2str(AlgName), count, iter, Creacm.recon);
680 fprintf(output, "\n reconstruction completed after iteration %4i", count);
682 // set value for the "previous basis" in the next execution : TRUE if pixel
683 pix_old_flag = pix_basis;
684 // if current execution is blob basis, record the values of the blob parameters,
685 // so that IF the next iteration is also blob basis, they parameters should not be changed.
688 support_old = Blob.support;
689 shape_old = Blob.shape;
690 delta_old = Blob.delta;
692 // reset algorithm after termination
693 if (AlgClasses[AlgIndex]->Reset() != 0)
695 fprintf(output, "\n **** algorithm %s failed to reset", int2str(algn[AlgIndex]));
702 { // if currently blob, pixelize it
704 Blob.blob2pix(Blob.recon, Creacm.recon);
707 for (i = 0; i < GeoPar.area; i++) recon_copy[i] = Creacm.recon[i];
709 // EXECUTE SMOOTHING OR CONTOURING AS APPROPRIATE
710 if ((optsw == 1) && (optflg[iter] != 0))
712 smooth(recon_copy, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
715 ////////// dump recon ////////////
716 fprintf(output, "\nexalg: algorithm %s result of iteration %i after smooth", int2str(algn[AlgIndex]), count);
717 for (i = 0; i < area; i++)
719 if ((i % 3) == 0) fprintf(output, "\n"); fprintf(output, " %36.30f", (double) recon[i]);
721 //////////////////////////////////
725 if ((optsw == 2) && (optflg[iter] != 0))
727 contur(recon_copy, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
729 ////////// dump recon ////////////
731 fprintf(output, "\nexalg: algorithm %s result of iteration %i after contur", int2str(algn[AlgIndex]), count);
732 for (i = 0; i < area; i++)
734 if ((i % 3) == 0) fprintf(output, "\n");
735 fprintf(output, " %36.30f", (double) recon[i]);
738 //////////////////////////////////
741 RecFile.WriteRec(head, int2str(AlgName), count, iter, recon_copy);
744 fprintf(output, "\n iteration %4i completed", count);
746 fprintf(output, "\n %10.3f seconds for this iteration", time);
753 fprintf(output, "\n %10.3f seconds for this iteration", time);
758 fprintf(output, "\n %10.3f seconds for all iterations", time);
763 if (list != NULL) delete[] list;
764 if (weight != NULL) delete[] weight;
765 if (recon_copy != NULL) delete[] recon_copy;