src/snark/exalg.cpp

   1 /*
   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) $
   7  $Author: zye $
   8  ***********************************************************
   9
  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)
  13  IS SATISFIED.
  14  */
  15
  16 #include <cstdio>
  17 #include <cstdlib>
  18 #include <cmath>
  19 #include "blkdta.h"
  20 #include "creacm.h"
  21 #include "geom.h"
  22 #include "term.h"
  23 #include "raysel.h"
  24 #include "uiod.h"
  25 #include "int2str.h"
  26 #include "second.h"
  27 #include "bldlst.h"
  28 #include "projfile.h"
  29 #include "alb1.h"
  30 #include "alb2.h"
  31 #include "alb3.h"
  32 #include "alb4.h"
  33 #include "alb5.h"
  34 #include "alp1.h"
  35 #include "alp2.h"
  36 #include "alp3.h"
  37 #include "alp4.h"
  38 #include "alp5.h"
  39 #include "emap.h"
  40 #include "lino.h"
  41 #include "art.h"
  42 #include "sart.h"
  43 #include "conv.h"
  44 #include "sirt.h"
  45 #include "dcon.h"
  46 #include "foru.h"
  47 #include "rfl.h"
  48 #include "back.h"
  49 #include "quad.h"
  50 #include "mart.h"
  51 #include "trm1.h"
  52 #include "trm2.h"
  53 #include "trm3.h"
  54 #include "trm4.h"
  55 #include "trm5.h"
  56 #include "smooth.h"
  57 #include "contur.h"
  58 #include "recfile.h"
  59 #include "infile.h"
  60 #include "blob.h"
  61 #include "uiod.h"
  62 #include "exalg.h"
  63 #include "consts.h"
  64 #include "bdhk.h"
  65 #include "pseudo.h"
  66 #include "superior.h"
  67
  68 using namespace std;
  69
  70 // macro to compute the number of elements in an array
  71 #define length(array) (sizeof(array)/sizeof(array[0]))
  72
  73 static int totalIters = 0;
  74
  75 class NameTable
  76 {
  77 public:
  78         INTEGER Length;INTEGER* Names;
  79
  80         INTEGER GetIndex(INTEGER word)
  81         {
  82                 for (int i = 0; i < Length; i++)
  83                 {
  84                         if (word == Names[i])
  85                         {
  86                                 return i;
  87                         }
  88                 }
  89                 return -1;
  90         }
  91
  92         INTEGER GetName(INTEGER index)
  93         {
  94                 return Names[index];
  95         }
  96 };
  97
  98 // Instantiate Algorithm classes
  99 alb1_class alb1;
 100 alb2_class alb2;
 101 alb3_class alb3;
 102 alb4_class alb4;
 103 alb5_class alb5;
 104 alp1_class alp1;
 105 alp2_class alp2;
 106 alp3_class alp3;
 107 alp4_class alp4;
 108 alp5_class alp5;
 109 emap_class emap;
 110 lino_class lino;
 111 art_class art;
 112 sart sart;
 113 conv_class conv;
 114 sirt_class sirt;
 115 dcon_class dcon;
 116 foru_class foru;
 117 rfl_class rfl;
 118 back_class back;
 119 quad_class quad;
 120 mart_class mart;
 121 bdhk_class bdhk;
 122
 123 // Algorithms Classes
 124 alg_class* AlgClasses[] =
 125 {
 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
 151 };
 152
 153 #define NUMBER_OF_ALGORITHMS length(AlgClasses)
 154
 155 // Algorithms Names
 156 static const INTEGER algn[NUMBER_OF_ALGORITHMS] =
 157 {
 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')
 181 };
 182
 183 static NameTable Algorithms = { NUMBER_OF_ALGORITHMS, (int*) algn };
 184
 185 // Blob Algorithms Names
 186 static const INTEGER algn_blob[] =
 187 {
 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')
 199 };
 200
 201 #define NUM_OF_ALG_BLOB length(algn_blob)
 202
 203 static NameTable BlobAlgorithms = { NUM_OF_ALG_BLOB, (int*) algn_blob };
 204
 205 void exalg()
 206 {
 207         static const INTEGER exec_codes[30] =
 208         {
 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
 215
 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
 240
 241                 CHAR2INT('s', 'm', 'o', 'o'),
 242                 CHAR2INT('c', 'o', 'n', 't') };
 243
 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
 255         INTEGER algSt = 4;
 256         INTEGER algCnt = 22;
 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;
 261
 262         INTEGER word;
 263         BOOLEAN eol;
 264         BOOLEAN flagphan = FALSE;
 265         CHAR head[81];
 266         INTEGER count;
 267
 268         INTEGER optsw;
 269         INTEGER optflg[51];
 270
 271         INTEGER AlgIndex;
 272
 273         BOOLEAN contin;
 274
 275         INTEGER i;
 276         REAL val;
 277         INTEGER AlgName;
 278         REAL thresh;
 279         REAL w1;
 280         REAL w2;
 281         REAL w3;
 282         INTEGER iter;
 283         REAL begin, start;
 284         REAL end, end_alg;
 285         REAL time;
 286
 287         BOOLEAN tflag;
 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
 292
 293         static BOOLEAN called = FALSE;
 294         static INTEGER* list;
 295         static REAL* weight;
 296
 297         REAL* recon;
 298         REAL* recon_copy;
 299         INTEGER area; //size of picture
 300         INTEGER wx; //size of picture in x direct
 301         INTEGER wy; //size of picture in y direct
 302
 303         if (ProjFile.Open("prjfil") != 0)
 304         {
 305                 fprintf(output, "\n **** unable to open prjfil");
 306                 fprintf(output, "\n **** program aborted\n");
 307                 exit(-1);
 308         };
 309
 310         ProjFile.ReadHeader(head, &NoisePar, &Spctrm, &Anglst);
 311
 312         // PERFORM INITIALIZATION ONLY ON THE FIRST CALL
 313         if (!called)
 314         {
 315
 316                 Creacm.recon = new REAL[GeoPar.area];
 317
 318                 if (!Blob.pix_basis)
 319                         Blob.recon = new REAL[Blob.area];
 320
 321         };
 322         recon_copy = new REAL[GeoPar.area];
 323
 324         RaySel.effcnt = 0;
 325
 326         if (!RaySel.useray)
 327         {
 328                 RaySel.numray = GeoPar.snrays;
 329         }
 330         else
 331         {
 332                 RaySel.numray = GeoPar.usrays;
 333         }
 334
 335         if (!RaySel.useray)
 336         {
 337                 RaySel.initnr = GeoPar.fsnray;
 338                 RaySel.endray = GeoPar.lsnray;
 339         }
 340         else
 341         {
 342                 RaySel.initnr = GeoPar.fusray;
 343                 RaySel.endray = GeoPar.lusray;
 344         }
 345
 346         RaySel.begray = RaySel.initnr;
 347         RaySel.curnr = RaySel.initnr;
 348
 349         if (RaySel.stype == 0)
 350         {
 351                 RaySel.curnr = RaySel.initnr - RaySel.rayinc;
 352         }
 353
 354         RaySel.initnp = 0;
 355         RaySel.curnp = RaySel.initnp;
 356
 357         if (RaySel.stype == 1)
 358         {
 359                 RaySel.curnp = RaySel.initnp - RaySel.prjinc;
 360         }
 361
 362         RaySel.rcnt1 = 0;
 363         RaySel.rcnt2 = 0;
 364
 365         // INITIALIZE OPTION FLAGS
 366         optsw = 0;
 367         for (i = 0; i < 50; i++)
 368         {
 369                 optflg[i] = 0;
 370         }
 371
 372         tflag = FALSE;
 373
 374         word = InFile.getwrd(FALSE, &eol, exec_codes, executeOp1Cnt + algCnt);
 375         contin = FALSE;
 376
 377         val = 0.0;
 378
 379         if (word == exec_codes[0])
 380         {
 381                 val = GeoPar.aveden;
 382                 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
 383         }
 384         else if (word == exec_codes[1])
 385         {
 386                 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
 387         }
 388         else if (word == exec_codes[2])
 389         {
 390                 contin = TRUE;
 391                 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
 392         }
 393         else if (word == exec_codes[3])
 394         {
 395                 flagphan = TRUE;
 396                 if (Creacm.test == NULL)
 397                 {
 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");
 400                         return;
 401                 }
 402                 word = InFile.getwrd(FALSE, &eol, &(exec_codes[algSt]), algCnt);
 403         }
 404
 405         // LOOK UP ROUTINE NAME
 406         // check algorithm name
 407         if ((AlgIndex = Algorithms.GetIndex(word)) < 0)
 408         {
 409                 fprintf(output, "\n **** unknown algorithm %s", int2str(word));
 410                 fprintf(output, "\n **** algorithm execution aborted");
 411                 ProjFile.Close();
 412                 return;
 413         }
 414
 415         AlgName = word;
 416
 417         BOOLEAN pix_basis = Blob.pix_basis;
 418
 419         //check for algorithm that goes with blobs
 420         if (!Blob.pix_basis)
 421         {
 422                 // check if selected algorithm is iterative
 423                 if (BlobAlgorithms.GetIndex(AlgName) < 0)
 424                 {
 425                         fprintf(output, "\n **** blob basis is not appropriate for algorithm %s, the algorithm is being executed using pixel basis instead.", int2str(AlgName));
 426
 427                         pix_basis = true;
 428                         Blob.pix_basis = true;
 429                 }
 430         }
 431
 432         // ACCEPT AND DISPLAY COMMENT CARD
 433         InFile.GetComment(head);
 434
 435         fprintf(output, "\n         %s", head);
 436
 437         word = InFile.getwrd(FALSE, &eol, &(exec_codes[executeOp2St]),
 438                         executeOp2Cnt);
 439
 440         // CHECK TO SEE IF OPTION HAS BEEN SPECIFIED
 441         if (!eol)
 442         {
 443                 if (word == exec_codes[executeOp2St])
 444                 {
 445                         //smooth
 446                         optsw = 1;
 447                 }
 448
 449                 if (word == exec_codes[executeOp2St + 1])
 450                 {
 451                         //contour
 452                         optsw = 2;
 453                 }
 454
 455                 thresh = InFile.getnum(TRUE, &eol);
 456                 w1 = InFile.getnum(FALSE, &eol);
 457                 w2 = InFile.getnum(FALSE, &eol);
 458                 w3 = InFile.getnum(FALSE, &eol);
 459
 460                 // ECHO OPTION SELECTED
 461                 if (optsw == 1)
 462                 {
 463                         fprintf(output, "\n         reconstruction is smoothed after");
 464                 }
 465
 466                 if (optsw == 2)
 467                 {
 468                         fprintf(output, "\n         reconstruction is contoured after");
 469                 }
 470
 471                 InFile.listit(optflg);
 472
 473                 if (optsw == 1)
 474                 {
 475                         if (pix_basis)
 476                         {
 477                                 fprintf(output, "\n         threshold = %9.5f      w1 = %9.5f      w2 = %9.5f      w3 = %9.5f", thresh, w1, w2, w3);
 478                         }
 479                         else
 480                         {
 481                                 fprintf(output, "\n         threshold = %9.5f      w1 = %9.5f      w2 = %9.5f      w3 = %9.5f", thresh, w1, w2, w3);
 482
 483                         }
 484                 }
 485
 486                 if (optsw == 2)
 487                 {
 488                         fprintf(output, "\n         threshold = %9.5f      w1 = %9.5f      w2 = %9.5f      w3 = %9.5f", thresh, w1, w2, w3);
 489                 }
 490         }
 491
 492         if (contin)
 493         {
 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)
 497                         {
 498                                 fprintf(output, "\n **** error: current blob parameters must match those of the previous blob basis execution");
 499                                 fprintf(output, "\n **** algorithm execution aborted");
 500                                 ProjFile.Close();
 501                                 return;
 502                         }
 503                 }
 504
 505                 if (!pix_old_flag && pix_basis)
 506                 { // previously blob, now pixel
 507                         Blob.blob2pix(Blob.recon, Creacm.recon);
 508                         delete[] Blob.recon;
 509                 }
 510
 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);
 515                 }
 516         }
 517         else
 518         {
 519                 if (!pix_basis)
 520                         Blob.recon = new REAL[Blob.area];
 521
 522                 // initiate recon with val or phantom image
 523                 if (flagphan == TRUE)
 524                 {
 525                         for (i = 0; i < GeoPar.area; i++)
 526                         {
 527                                 Creacm.recon[i] = Creacm.test[i];
 528                         }
 529                         if (!pix_basis)
 530                         {
 531                                 Blob.pix2blob(Creacm.recon, Blob.recon);
 532                         }
 533                 }
 534                 else
 535                 {
 536                         if (pix_basis)
 537                         {
 538                                 for (i = 0; i < GeoPar.area; i++)
 539                                 {
 540                                         Creacm.recon[i] = val;
 541                                 }
 542                         }
 543                         else
 544                         {
 545                                 for (i = 0; i < Blob.area; i++)
 546                                 {
 547                                         Blob.recon[i] = val;
 548                                 }
 549                         }
 550                 }
 551         }
 552
 553         if (pix_basis)
 554         {
 555                 bldlst(&list, &weight);
 556                 recon = Creacm.recon;
 557                 wx = GeoPar.nelem;
 558                 wy = GeoPar.nelem;
 559                 area = GeoPar.area;
 560         }
 561         else
 562         {
 563                 Blob.bbldlst(&list, &weight);
 564                 recon = Blob.recon;
 565                 wx = Blob.M;
 566                 wy = Blob.H;
 567                 area = Blob.area;
 568         }
 569
 570         called = TRUE;
 571
 572         second(&begin);
 573         start = begin;
 574
 575         // initialize selected algorithm
 576         if (AlgClasses[AlgIndex]->Init() != 0)
 577         {
 578                 fprintf(output, "\n **** algorithm %s failed to initialize", int2str(algn[AlgIndex]));
 579                 fprintf(output, "\n **** algorithm execution aborted");
 580                 ProjFile.Close();
 581                 return;
 582         }
 583
 584         count = 0;
 585
 586         // run until the termination criterion is met or an error occurs
 587         while (TRUE)
 588         {
 589                 count++;
 590                 totalIters++;
 591                 iter = ((count - 1) % 50) + 1;
 592
 593                 // execute single iteration of the selected algorithm
 594                 if (SuperSet.superiorizationEnabled)
 595                 {
 596                         tflag = executeSuperiorization(recon, list, weight, count, AlgClasses[AlgIndex]);
 597                 }
 598                 else
 599                 {
 600                         tflag = AlgClasses[AlgIndex]->Run(recon, list, weight, count);
 601                 }
 602
 603                 second(&end_alg);
 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);
 607
 608                 // if error
 609                 if (tflag)
 610                 {
 611                         fprintf(output, "\n **** algorithm indicated failure to continue");
 612
 613                         // reset algorithm after failure
 614                         if (AlgClasses[AlgIndex]->Reset() != 0)
 615                         {
 616                                 fprintf(output, "\n **** algorithm %s failed to reset", int2str(algn[AlgIndex]));
 617                         }
 618
 619                         fprintf(output, "\n **** algorithm execution aborted");
 620                         break;
 621                 }
 622
 623                 if (Term.termInstance == NULL)
 624                 {
 625                         tflag = ((count - Term.iterations) >= 0);
 626                 }
 627                 else
 628                 {
 629                         tflag = Term.termInstance->Run(recon, list, weight, count);
 630                 }
 631
 632                 if (tflag)
 633                 {
 634                         // END THE ROUTINE IF THE TEST WAS SATISFIED
 635                         // CHECK IF OPTION SPECIFIED FOR LAST ITERATION
 636
 637                         iter = 0;
 638
 639                         if (!pix_basis)
 640                         { // if currently blob, pixelize it
 641
 642                                 Blob.recon = recon;
 643                                 Blob.blob2pix(Blob.recon, Creacm.recon);
 644                         };
 645
 646                         if ((optsw == 1) && (optflg[iter] != 0))
 647                         {
 648                                 smooth(Creacm.recon, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
 649
 650 #ifdef FFCOMPARE
 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++)
 654                                 {
 655                                         if ((i % 3) == 0) fprintf(output, "\n");
 656                                         fprintf(output, " %36.30f", (double) recon[i]);
 657                                 }
 658                                 //////////////////////////////////
 659 #endif
 660                         };
 661
 662                         if ((optsw == 2) && (optflg[iter] != 0))
 663                         {
 664                                 contur(Creacm.recon, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
 665
 666 #ifdef FFCOMPARE
 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++)
 670                                 {
 671                                         if ((i % 3) == 0) fprintf(output, "\n");
 672                                         fprintf(output, " %36.30f", (double) recon[i]);
 673                                 }
 674                                 //////////////////////////////////
 675 #endif
 676                         };
 677
 678                         RecFile.WriteRec(head, int2str(AlgName), count, iter, Creacm.recon);
 679
 680                         fprintf(output, "\n         reconstruction completed after iteration %4i", count);
 681
 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.
 686                         if (!pix_basis)
 687                         {
 688                                 support_old = Blob.support;
 689                                 shape_old = Blob.shape;
 690                                 delta_old = Blob.delta;
 691                         }
 692                         // reset algorithm after termination
 693                         if (AlgClasses[AlgIndex]->Reset() != 0)
 694                         {
 695                                 fprintf(output, "\n **** algorithm %s failed to reset", int2str(algn[AlgIndex]));
 696                         }
 697
 698                         break;
 699                 }
 700
 701                 if (!pix_basis)
 702                 { // if currently blob, pixelize it
 703                         Blob.recon = recon;
 704                         Blob.blob2pix(Blob.recon, Creacm.recon);
 705                 };
 706
 707                 for (i = 0; i < GeoPar.area; i++) recon_copy[i] = Creacm.recon[i];
 708
 709                 // EXECUTE SMOOTHING OR CONTOURING AS APPROPRIATE
 710                 if ((optsw == 1) && (optflg[iter] != 0))
 711                 {
 712                         smooth(recon_copy, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
 713
 714 #ifdef FFCOMPARE
 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++)
 718                         {
 719                                 if ((i % 3) == 0) fprintf(output, "\n"); fprintf(output, " %36.30f", (double) recon[i]);
 720                         }
 721                         //////////////////////////////////
 722 #endif
 723                 };
 724
 725                 if ((optsw == 2) && (optflg[iter] != 0))
 726                 {
 727                         contur(recon_copy, GeoPar.nelem, GeoPar.nelem, thresh, w1, w2, w3);
 728
 729                         ////////// dump recon ////////////
 730 #ifdef FFCOMPARE
 731                         fprintf(output, "\nexalg: algorithm %s result of iteration %i after contur", int2str(algn[AlgIndex]), count);
 732                         for (i = 0; i < area; i++)
 733                         {
 734                                 if ((i % 3) == 0) fprintf(output, "\n");
 735                                 fprintf(output, " %36.30f", (double) recon[i]);
 736                         }
 737 #endif
 738                         //////////////////////////////////
 739                 }
 740
 741                 RecFile.WriteRec(head, int2str(AlgName), count, iter, recon_copy);
 742
 743                 second(&end);
 744                 fprintf(output, "\n         iteration %4i completed", count);
 745                 time = end - begin;
 746                 fprintf(output, "\n         %10.3f seconds for this iteration", time);
 747                 begin = end;
 748
 749         }
 750
 751         second(&end);
 752         time = end - begin;
 753         fprintf(output, "\n         %10.3f seconds for this iteration", time);
 754
 755         if (count > 1)
 756         {
 757                 time = end - start;
 758                 fprintf(output, "\n         %10.3f seconds for all iterations", time);
 759         };
 760
 761         ProjFile.Close();
 762
 763         if (list != NULL) delete[] list;
 764         if (weight != NULL)     delete[] weight;
 765         if (recon_copy != NULL) delete[] recon_copy;
 766
 767         return;
 768 }