src/snark/quad.cpp

   1 /*
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                           *
   4  *                              S N A R K   1 4                              *
   5  *                                                                           *
   6  *                     A PICTURE RECONSTRUCTION PROGRAM                      *
   7  *                                                                           *
   8  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   9
  10  quad.cpp,v 1.6 2009/06/01 03:33:26 jklukowska Exp
  11
  12  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **
  13  *                                                                *
  14  *      QUADRATIC ITERATIVE RECONSTRUCTION ALGORITHMS             *
  15  *                                                                *
  16  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **
  17
  18
  19     PROGRAM BY EHUD ARTZY 1977
  20
  21      INTRODUCTION
  22      ============
  23
  24      THIS PROGRAM IS SET TO MINIMIZE THE FUNCTION
  25
  26      K(X)=SA*(P-M*X)#*A*(P-M*X)+(X-XBAR)#*(SB*B**L+SC*C**(-1)*(X-XBAR)
  27
  28      WHERE
  29
  30      M IS AN MM BY NN MATRIX
  31      P AN MM-DIMENSIONAL VECTOR
  32      A IS A NONNEGATIVE DIAGONAL MATRIX OF RANK MM
  33      B,C ARE POSITIVE DEFINITE SYMETRIC NN BY NN MATRICES
  34      SA,SB,SC ARE NONNEGATIVE REAL NUMBERS
  35      XBAR IS AN NN-DIMENSIONAL VECTOR
  36      L IS A NONNEGATIVE INTEGER
  37      # STANDS FOR TRANSPOSE
  38
  39      IF SB = SC = 0 THEN AN ADDITIONAL DIAGONAL MATRIX D OF RANK NN IS
  40      INTRODUCED WHICH, FOR ALL X WHICH MINIMIZE K(X), WILL CHOOSE THE
  41      ONE WHICH MINIMIZES THE NORM OF D**(-1)*(X-XBAR)
  42
  43      IF SB + SC GREATER THAN 0 THEN
  44
  45         LET
  46
  47    (1)  Q = C**(T/2)*((SA*M#*A*M +SB*B**L)*C**(T/2) +(SC*C**(T/2 - 1))
  48
  49         F = C**(T/2)*SA*M#*A*(P - M*XBAR)
  50
  51         IF Y IS A SOLUTION OF Q*Y = F THEN  X = C**(T/2)*Y + XBAR
  52         MINIMIZES K(X)
  53
  54      ELSE LET
  55
  56    (2)   Q = D*M#*A*M*D
  57
  58          F = D*M#*A*P
  59
  60          IF Y IS A SOLUTION OF Q*Y = F THEN X = D*Y IS A
  61          MINIMUM POINT OF K(X)
  62
  63      FOR FURTHER DETAILS SEE:
  64            QUADRATIC OPTIMIZATION FOR IMAGE RECONSTRUCTION.
  65            GABOR T. HERMAN AND ARNOLD LENT
  66            COMPUTER GRAPHICS AND IMAGE PROCESSING 5,319-332 (1976)
  67
  68       THE FOLLOWING ALGORITHM IS USED TO FIND THE SOLUTION OF Q*Y = F
  69
  70       G(0) = R(0) = Q*Y(0) - F
  71       Y(K+1) = Y(K) - DELTA(K)*G(K)
  72       R(K+1) = Q*Y(K+1) - F
  73       G(K+1) = R(K+1) + GMA(K)*G(K)
  74
  75       AN ALTERNATIVE WAY TO COMPUTE R(K+1) IS GIVEN BY
  76       R(K+1) = R(K) - DELTA(K)*Q*G(K)
  77
  78       THE FORMER IS REFFERED TO AS DIRECT METHOD AND THE LATER INDIRECT
  79       METHOD
  80
  81       THE CHOICE OF DELTA(K) AND GMA(K) DETERMINE THE PARTICULAR KIND
  82       OF ALGORITHM USED. WHEN GMA(K) = 0 FOR EVERY K WE HAVE A ONE STEP
  83       METHOD (OSM) OTHERWISE WE HAVE TWO STEPS METHOD (TSS)
  84
  85       FOR USING THE PROGRAM THE USER MUST SPECIFY THE FOLLOWING
  86       PARAMETERS:
  87
  88      CARD 1 - AOPT,BOPT,COPT,MINZ,PERIOD,TOLER,DELTA,GMA
  89
  90       AOPT  : 1 - OSM WITH DIRECT COMPUTATION OF GRADIENT
  91               2 - OSM WITH INDIRECT COMPUTATION OF GRADIENT
  92               3 - TSS WITH DIRECT COMPUTATION OF GRADIENT
  93               4 - TSS WITH INDIRECT COMPUTATION OF GRADIENT
  94
  95       BOPT  : 0 - THIS IS A STATIONARY METHOD FOR OSM USER SUPPLIES
  96                   DELTA FOR TSS USER SUPPLIES DELTA AND GMA
  97               1 - THIS IS A STATIONARY METHOD WITH DELTA COMPUTED BY
  98                  THE PROGRAM TO BE
  99
 100                          DELTA = 2.0/(LMD + ESMIN)
 101                  WHERE
 102                    LMD AND ESMIN ARE ESTIMATES OF LARGEST AND SMALLEST
 103                    POSITIVE EIGENVALUES OF THE MATRIX Q.
 104              -1 -  USER SUPPLY LMD AND ESMIN FOR BOPT 1
 105               2 -  FOR OSM THIS IS RICHARDSON METHOD OF LENGTH PERIOD
 106                    FOR TSS THIS IS SEMI ITERATIVE METHOD BASED ON
 107                    RICHARDSON STATIONARY METHOD. LMD AND ESMIN(SEE 1)
 108                    ARE COMPUTED BY THE PROGRAM FOR REFERENCE SEE:
 109                    DAVID M. YOUNG
 110                    ITERATIVE SOLUTION OF LARGE LINEAR SYSTEMS
 111                    ACADEMIC PRESS 1971
 112              -2 -  USER SUPPLY LMD AND ESMIN FOR BOPT 2
 113               3 -  FOR OSM THIS IS THE STEEPEST DESCENT METHOD WHERE
 114                    DELTA(K) = R(K)#*R(K)/(R(K)#*Q*R(K))
 115                    FOR TSS THIS IS THE CONJUGATE GRADIENT METHOD WHERE
 116                    DELTA(K) = R(K)#*R(K)/(G(K)#*Q*G(K))
 117                    GMA(K) = R(K)#*R(K)/(R(K-1)#*R(K-1))
 118
 119       COPT  : 1 -  EQUATION 2 IS USED WITH GOITEIN D MATRIX
 120                    FOR REFERENCE SEE:
 121                    M. GOITEIN
 122                    THREE-DIMENSIONAL DENSITY RECONSTRUCTION FROM A
 123                    SERIES OF TWO-DIMENSIONAL PROJECTIONS
 124                    NUCL. INSTR. METHODS 101, 509-518 (1972)
 125               2 -  EQUATION 2 IS USED WITH LSIRT D MATRIX
 126                    FOR REFERENCE SEE:
 127                    A. V. LAKSHMINARAYNAN AND A. LENT
 128                    THE SIMULTANEOUS ITERATIVE RECONSTRUCTION TECHNIQUE
 129                    AS LEAST SQUARES METHOD
 130                    SPIE. VOL. 96 OPTICAL INSTRUMENTATION IN
 131                    MEDICINE V (1976)
 132               3 -  THE USER SUPPLY HIS OWN DIAGONAL D MATRIX. THE
 133                    VALUES OF THE DIAGONAL ELEMENTS ARE COMPUTED BY
 134                    DSET SQUARED AND STORED IN A VECTOR OF LENGTH
 135                    AREA (SEE SNARK MANUAL). THE CALLING SEQUENCE
 136                    FOR DSET IS
 137                    CT(D,LIST,WEIGHT)
 138                    WHERE
 139                    D     -CALCULATED VECTOR
 140                    LIST  -SEE THE SNARK MANUAL
 141                    WEIGHT-SEE THE SNARK MANUAL
 142               4 -  EQUATION 1 IS USED CARD 3-5 MUST BE PRESENTED
 143
 144       MINZ  : 0 -  FOR BOPT 1 OR 2 ESMIN IS ASSUMED TO BE 0
 145               1 -  FOR BOPT 1 OR 2 ESMIN IS CALCULATED
 146
 147       PERIOD: N -  N IS THE LENGTH OF THE PERIOD FOR THE RICHARDSON OSM
 148                    N MUST BE IN THE RANGE 1-32
 149
 150       TOLER :   -  THIS IS THE TOLERANCE FOR WHICH LMD AND ESMIN IS
 151                    COMPUTED FOR REFERENCE SEE:
 152                    A COMPUTER INPLEMENTATION OF BAYESIAN ANALYSIS
 153                    OF IMAGE RECONSTRUCTION
 154                    G. T. HERMAN AND A. LENT
 155                    INFORMATION AND CONTROL 31,364-384 (1976)
 156
 157       DELTA :   -  FOR BOPT 0 THIS IS CONSTANT DELTA
 158                    FOR BOPT -1,-2 THIS IS USER SUPPLY ESMIN
 159
 160       GMA   :   -  FOR BOPT 0 AND TSS THIS IS CONSTANT GMA
 161                    FOR BOPT -1,-2 THIS IS USER SUPPLY LMD
 162
 163      CARD 2 - DOPT,EOPT,FOPT,GOPT,HOPT,POPT,ERROR,MINPRD,INTERP
 164
 165       DOPT  : 1 -  NATURAL ORDER FOR PICKING DELTA FOR RICHARDSON OSM
 166                    (AOPT LE 2 BOPT 2)
 167               2 -  USER SUPPLY ORDER OF PICKING DELTA CARD 6 MUST BE
 168                    PRESENTED
 169               3 - LEBEDTV-FINOGENOV ORDERING FOR PICKING DELTA
 170                   FOR REFERENCE SEE:
 171                   R.S. ANDERSSEN AND G.H. GOLUB
 172                   RICHARDSONS NON-STATIONARY MATRIX ITERATIVE PROCEDURE
 173                   STAN-CS-72-307 (1972)
 174                   STANFORD UNIVERSITY
 175       EOPT  : 1 - A IS THE IDENTITY MATRIX TIMES THE CONSTANT ERROR
 176               2 - THE VALUE OF THE MATRIX A IS COMPUTED BY ESTIMATING
 177                   THE VARIANCE IN THE MEASUREMENTS. THIS OPTION IS
 178                   AVAILABLE ONLY FOR NOISY DATA (SEE SNARK MANUAL FOR
 179                   NOISE ESTIMATION).
 180               3 - USER SUPPLY AN ERROR ROUTINE. THE CALLING SEQUENCE
 181                   IS
 182                     PUNCTION UERROR(NP,NR,PD)
 183                     WHERE
 184                     NP - IS A PROJECTION NUMBER
 185                     NR - IS A RAY NUMBER IN PROJECTION NP
 186                     PD - IS THE VALUE OF THE (NP,NR) PROJECTION
 187
 188       FOPT  : 1 - XBAR (EXPECTED VALUE) IS THE ZERO VECTOR
 189               2 - XBAR IS AVERAGE DENSITY AS APROXIMATE BY SNARK
 190               3 - XBAR IS NORMALIZE BACK PROJECTION
 191               4 - XBAR IS CONTINUOUS BACK PROJECTION
 192               5 - XBAR IS STARTING POINT
 193       GOPT  : 0 - NO SPLITTING
 194               1 - DIAGONAL SPLITTING
 195               2 - USER SUPPLY SPLITTING MATRIX
 196                     CALL SSET(S,LIST,WEIGHT)
 197             WHERE
 198             S IS A DIAGONAL SPLITTING MATRIX
 199             LIST SEE THE SNARK05 MANUAL
 200             WEIGHT SEE THE SNARK05 MANUAL
 201
 202       HOPT  : 1 - RECONSTRUCTED PICTURE IS CORRECTED SO THAT AVERAGE
 203                   DENSITY WILL BE EQUAL TO THAT OF AVERAGE DENSITY
 204                   COMPUTED BY SNARK
 205               0 - NO CORRECTION
 206
 207       POPT  : 0 - NO REPORTS
 208               1 - THE FOLLOWING VALUES ARE REPORTED
 209                   GRADIENT (ONLY FOR AOPT GE 2 OR BOPT 3)
 210                   NORM OF CORRECTION,DELTA,NORM OF CORRECTING VECTOR
 211               2 - IN ADDITION COST FUNCTION IF POSSIBLE
 212
 213       ERROR :   - FOR EOPT 1 THE VALUE OF THE DIAGONAL ELEMENT
 214       MINPRD:   - IF MULTIPLICATIVE NOISE IS PRESENTED AND EOPT 2
 215                   A CONSTANT TO USE WHEN MEASURMENT IS ZERO MUST BE
 216                   SPECIFIED
 217
 218       INTERP:   - FOR FOPT 4 THIS IS THE INTERPOLATION INDEX
 219                (SEE THE SNARK05 MANUAL)
 220
 221       THE FOLLOWING CARD(S) IS(ARE) READ ONLY IF BOPT = 2
 222
 223             ORDER(I), I = 1,PERIOD
 224       ORDER :  PERIOD NUMBERS IN THE RANGE 1 TO PERIOD
 225
 226       THE FOLLOWING 3 CARDS ARE READ ONLY IF COPT GE 4
 227
 228      CARD 3 - SA,SB,SC
 229
 230       SA    :   CONSTANT (SEE EQUATION 1)
 231       SB    :   CONSTANT (SEE EQUATION 1)
 232       SC    :   CONSTANT (SEE EQUATION 1)
 233
 234      CARD 4 - L,BW1,BW2,BBCON          *** MATRIX B ***
 235
 236       L     : POWER OF MATRIX B
 237       BW1   : VARIANCE OF DENSITY IN AN ARBITRARY PIXEL
 238       BW2   : COVARIANCE OF DENSITIES IN NEIGHBORING PIXELS
 239       BW3   : COVARIANCE OF DENSITIES IN DIAGONALY NEIGHBORING PIXELS
 240       BBCON : KAMI - GENERALIZED KASHYAP-MITTAL SMOOTING (L MUST BE 1)
 241                      FOR REFERENCE SEE:
 242                      R.L. KASHYAP AND M.C. MITTLE
 243                      PICTURE RECONSTRUCTION FROM PROJECTIONS
 244                      1EEE TRANS. COMPUTERS, C-27,915-923 (1975)
 245               SAJZ - SMOTH AS IN SNAR WITHOUT NORMALIZATION
 246               SAJC - SAME AS SAJZ EXCEPT THAT BOARDER PIXELS ASSUME
 247                      THE SAME DENSITY OUTSIDE THE PICTURE
 248               SNAR - SMOTH AS IN SNARK(SEE SNARK MANUAL) FOR SNAR
 249                      B IS NOT A SYMETRIC MATRIX
 250
 251      CARD 5 - T,CW1,CW2,CW3,CBCON      *** MATRIX C ***
 252
 253       T     : T/2 IS THE POWER OF MATRIX C IN EQUATION 1
 254       CW1   : SEE BW1
 255       CW2   : SEE BW2
 256       CW3   : SEE BW3
 257       CBCON : SEE BBCON
 258
 259
 260
 261 ***********************************************************************
 262 */
 263
 264 #include <cstdlib>
 265 #include <cstdio>
 266 #include <cmath>
 267
 268 #include "blkdta.h"
 269 #include "geom.h"
 270 #include "spctrm.h"
 271 #include "fourie.h"
 272 #include "anglst.h"
 273 #include "bckprj.h"
 274 #include "int2str.h"
 275 #include "modefl.h"
 276 #include "consts.h"
 277 #include "uiod.h"
 278 #include "projfile.h"
 279 #include "infile.h"
 280 #include "errpar.h"
 281 #include "err.h"
 282 #include "blob.h"
 283
 284 #include "quad.h"
 285
 286 INTEGER quad_class::Init()
 287 {
 288
 289   cg_ = NULL;
 290   g_ = NULL;
 291   v_ = NULL;
 292   w_ = NULL;
 293   r_ = NULL;
 294   f_ = NULL;
 295   fb_ = NULL;
 296   s_ = NULL;
 297   d_ = NULL;
 298   m = NULL;
 299   n = NULL;
 300
 301   del = NULL;               //bug 174, wei, 10/2005
 302   order = NULL;
 303
 304   return 0;                      //wei, 3/2005
 305 }
 306
 307 INTEGER quad_class::Reset()
 308 {
 309   if(m != NULL&& m != n) delete[] m;
 310   if(n != NULL&& n != w_) delete[] n;                 //bug 191, wei, 12/2005
 311   if(cg_ != NULL&&cg_!=g_) delete[] cg_;
 312   if(g_ != NULL&&g_!=r_) delete[] g_;
 313   if(v_ != NULL&&v_!=r_) delete[] v_;
 314   if(w_ != NULL&&w_!=r_) delete[] w_;
 315   if(r_ != NULL) delete[] r_;
 316   if(f_ != NULL) delete[] f_;
 317   if(fb_ != NULL) delete[] fb_;
 318   if(s_ != NULL) delete[] s_;
 319   if(d_ != NULL) delete[] d_;
 320   if(del != NULL) delete[] del;   //bug 174, wei, 10/2005
 321   if(order != NULL) delete[] order;
 322
 323
 324
 325   return 0;                    //wei, 3/2005
 326 }
 327
 328 BOOLEAN quad_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
 329 {
 330   static INTEGER minz;
 331
 332   static REAL delta;
 333   REAL toler;
 334   static REAL gma;//set static. hstau see deltac()  below
 335   //INTEGER popt;  poot is global, why is it declared here? hstau
 336   INTEGER ir;
 337   //REAL ratio;
 338   //INTEGER iratio;
 339   INTEGER i;
 340
 341   ///INTEGER* order;
 342   //REAL error;
 343   //REAL sa;
 344   //REAL sb;
 345   //REAL sc;
 346   //INTEGER l;
 347   //INTEGER t;
 348   //REAL bw1;
 349   //REAL bw2;
 350   //REAL bw3;
 351   //INTEGER bbcon;
 352   //REAL cw1;
 353   //REAL cw2;
 354   //REAL cw3;
 355   //INTEGER cbcon;
 356   //REAL xbar;
 357   ///INTEGER iterx;
 358   //INTEGER aopt;
 359   //INTEGER bopt;
 360   //INTEGER copt;
 361   //INTEGER period;
 362   //INTEGER dopt;
 363   //INTEGER eopt;
 364   //INTEGER fopt;
 365   //INTEGER gopt;
 366   //INTEGER  hopt;
 367
 368   INTEGER k;
 369   REAL dumy;
 370   REAL dum = 0.0; // this is dummy variable
 371   //INTEGER js;
 372   static REAL esmin; // needed to be set static? hstau
 373   INTEGER j;
 374
 375   INTEGER np;
 376   REAL theta;
 377   REAL sinth;
 378   REAL costh;
 379   REAL spac;
 380
 381
 382   BOOLEAN xalg8_tmp;
 383
 384 /*
 385   static REAL* cg_;//set static. hstau see deltad()  below
 386   static REAL* g_;//set static. hstau see deltad()  below
 387   static REAL* v_;
 388   static REAL* w_;  //set static. hstau
 389   static REAL* r_;
 390   static REAL* f_;  //set static. hstau see resedu() and prpict()  below
 391   static REAL* fb_;  //for blob case;
 392   static REAL* s_; //set static. hstau
 393   static REAL* d_; //set static. hstau                     //moved to class quad,wei,3/2005
 394 */
 395
 396   static REAL lmd; // needed to be set static? hstau
 397   ///dimension init(20),dmart(6),let(8)
 398   INTEGER headng;
 399   static BOOLEAN alg;//set static. hstau
 400   BOOLEAN eol;
 401
 402   /*
 403   //data let / 1ha,1hb,1hc,1hd,1he,1hf,1hg,1hh /
 404    static const INTEGER let[] = {// needed to be set static? hstau
 405      CHAR2INT(' ',' ',' ',' '),   //added to match Fortran indexing .hstau
 406      CHAR2INT('a',' ',' ',' '),
 407      CHAR2INT('b',' ',' ',' '),
 408      CHAR2INT('c',' ',' ',' '),
 409      CHAR2INT('d',' ',' ',' '),
 410      CHAR2INT('e',' ',' ',' '),
 411      CHAR2INT('f',' ',' ',' '),
 412      CHAR2INT('g',' ',' ',' '),
 413      CHAR2INT('h',' ',' ',' ')
 414    };
 415    */
 416   const char *let[] = {
 417     "algorithm",
 418     "delta_flag",
 419     "d_flag",
 420     "order",
 421     "error_type",
 422     "exp_value",
 423     "split",
 424     "normalize",
 425     "period",
 426     "toler",
 427     "print",
 428     "error",
 429     "minv",
 430     "interp"};
 431
 432   //data dmart /4hgoit,4hein ,4hsirt,4h    ,4huser,4h    /
 433 //  static const INTEGER dmart[] = {// needed to be set static? hstau
 434 //    CHAR2INT(' ',' ',' ',' '),   //added to match Fortran indexing .hstau
 435 //    CHAR2INT('g','o','i','t'),
 436 //    CHAR2INT('e','i','n',' '),
 437 //    CHAR2INT('s','i','r','t'),
 438 //    CHAR2INT(' ',' ',' ',' '),
 439 //    CHAR2INT('u','s','e','r'),
 440 //    CHAR2INT(' ',' ',' ',' ')
 441 //  };
 442
 443   ///data init / 4hzero,4h vec,4htor ,4h    ,4haver,4hage ,4hdens,
 444   ///            4hity ,4hback,4h pro,4hject,4hion , 4hcont,4h bac,
 445   ///            4hk pr,4hojs ,4hinit,4hial ,4hvect,4hor  /
 446
 447   static const INTEGER init[] = {// needed to be set static? hstau
 448     CHAR2INT(' ',' ',' ',' '),   //added to match Fortran indexing .hstau
 449     CHAR2INT('z','e','r','o'),
 450     CHAR2INT(' ','v','e','c'),
 451     CHAR2INT('t','o','r',' '),
 452     CHAR2INT(' ',' ',' ',' '),
 453     CHAR2INT('a','v','e','r'),
 454     CHAR2INT('a','g','e',' '),
 455     CHAR2INT('d','e','n','s'),
 456     CHAR2INT('i','t','y',' '),
 457     CHAR2INT('b','a','c','k'),
 458     CHAR2INT(' ','p','r','o'),
 459     CHAR2INT('j','e','c','t'),
 460     CHAR2INT('i','o','n',' '),
 461     CHAR2INT('c','o','n','t'),
 462     CHAR2INT(' ','b','a','c'),
 463     CHAR2INT('k',' ','p','r'),
 464     CHAR2INT('o','j','s',' '),
 465     CHAR2INT('i','n','i','t'),
 466     CHAR2INT('i','a','l',' '),
 467     CHAR2INT('v','e','c','t'),
 468     CHAR2INT('o','r',' ',' ')
 469   };
 470
 471   static const INTEGER hbb = CHAR2INT(' ',' ',' ',' ');// needed to be set static? hstau
 472   static const INTEGER hin = CHAR2INT(' ',' ','i','n');// needed to be set static? corr. move i n to the last two positions hstau
 473
 474   // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 16, 2004
 475   static const INTEGER quad_blurck_codes[4] =
 476   {
 477     CHAR2INT('s','n','a','r'),
 478     CHAR2INT('k','a','m','i'),
 479     CHAR2INT('s','a','j','z'),
 480     CHAR2INT('s','a','j','c')
 481   };
 482
 483
 484   INTEGER area;
 485
 486
 487   iterx = iter;
 488
 489
 490   if(Blob.pix_basis) {
 491     area = GeoPar.area;
 492   }
 493   else {
 494     area = Blob.area;
 495   }
 496
 497  ///tracer
 498   //   for(i = 0; i < area; i++) {
 499   //  fprintf(output,"\nrecon=%f", recon[i]);
 500   //}
 501       ////tracer
 502
 503 // THE FOLLOWING SECTION IS EXECUTED ONLY AT FIRST ITERATION
 504
 505   if(!(iter != 1)) {
 506
 507     // INITILIZE DUMMY VARIABLE
 508
 509     xalg8_tmp = FALSE;
 510     alg = FALSE;
 511     ///  d = 1;
 512     ///  f = 1;
 513     ///  s = 1;
 514
 515     // READ INPUT
 516
 517     ///aopt = getint(input, &eol);
 518     ///bopt = getint(NULL, &eol);
 519     ///copt = getint(NULL, &eol);
 520     ///minz = getint(NULL, &eol);
 521     ///period = getint(NULL, &eol);
 522     ///toler = getnum(NULL, &eol);
 523     ///delta = getnum(NULL, &eol);
 524     ///gma = getnum(NULL, &eol);
 525     ///dopt = getint(input, &eol);
 526     ///eopt = getint(NULL, &eol);
 527     ///fopt = getint(NULL, &eol);
 528     ///gopt = getint(NULL, &eol);
 529     ///hopt = getint(NULL, &eol);
 530     ///popt = getint(NULL, &eol);
 531     ///error_1 = getnum(NULL, &eol);
 532     ///Err.minprd = getnum(NULL, &eol);
 533     ///Fourie.interp = getint(NULL, &eol);
 534
 535     aopt = InFile.getint(TRUE, &eol);
 536     bopt = InFile.getint(FALSE, &eol);
 537     copt = InFile.getint(FALSE, &eol);
 538     minz = InFile.getint(FALSE, &eol);
 539     period = InFile.getint(FALSE, &eol);
 540     toler = InFile.getnum(FALSE, &eol);
 541     delta = InFile.getnum(FALSE, &eol);
 542     gma = InFile.getnum(FALSE, &eol);
 543
 544     if(eol) {
 545       fprintf(output, "\n         *** insufficient data item - end of line encountered ***");
 546       xalg8_tmp = TRUE;
 547       return xalg8_tmp;
 548     }
 549
 550     dopt = InFile.getint(TRUE, &eol);
 551     eopt = InFile.getint(FALSE, &eol);
 552     fopt = InFile.getint(FALSE, &eol);
 553     gopt = InFile.getint(FALSE, &eol);
 554     hopt = InFile.getint(FALSE, &eol);
 555     popt = InFile.getint(FALSE, &eol);
 556     error_1 = InFile.getnum(FALSE, &eol);
 557     Err.minprd = InFile.getnum(FALSE, &eol);
 558     Fourie.interp = InFile.getint(FALSE, &eol);
 559     //tracer
 560     //fprintf(output,"dopt=%d eopt=%d fopt=%d  gopt=%d hopt=%d popt=%d error=%f minv=%f interp=%d", dopt, eopt, fopt, gopt, hopt, popt, error_1, Err.minprd, Fourie.interp);
 561     //tracer
 562     if(eol) {
 563       fprintf(output, "\n         *** insufficient data item - end of line encountered ***");
 564       xalg8_tmp = TRUE;
 565       return xalg8_tmp;
 566     }
 567
 568     // CHECK RANGE OF OPTIONS
 569
 570     // Add a general checking regardless whether the parameter will be used or not,
 571     // reflecting the table summrizing all the restrictions on the QUAD variables in the manual
 572     // bug 191, wei, 12/2006
 573
 574     if(!((aopt>=1 ) &&(aopt<=4)))
 575         { fprintf(output, "\n         *** option algorithm is out of range ***");
 576       xalg8_tmp = TRUE;
 577       return xalg8_tmp;}
 578
 579     if(!((bopt>= -2) &&(bopt<=3)))
 580         {fprintf(output, "\n         *** option delta-flga is out of range ***");
 581           xalg8_tmp = TRUE;
 582       return xalg8_tmp;}
 583
 584     if(!((copt>= 1) &&(copt<=4)))
 585         {fprintf(output, "\n         *** option d-flag is out of range ***");
 586           xalg8_tmp = TRUE;
 587       return xalg8_tmp;}
 588
 589     if(!(period>=1 ))
 590         {fprintf(output, "\n         *** option period is out of range ***");
 591           xalg8_tmp = TRUE;
 592       return xalg8_tmp;}
 593
 594     if(!(toler> Consts.zero))
 595         {fprintf(output, "\n         *** option toler is out of range ***");
 596           xalg8_tmp = TRUE;
 597       return xalg8_tmp;}
 598
 599     if(!((dopt>=1 ) &&(dopt<=3)))
 600         {fprintf(output, "\n         *** option order is out of range ***");
 601           xalg8_tmp = TRUE;
 602       return xalg8_tmp;}
 603
 604     if(!((eopt>=1 ) &&(eopt<=3)))
 605         {fprintf(output, "\n         *** option error-type is out of range ***");
 606           xalg8_tmp = TRUE;
 607       return xalg8_tmp;}
 608
 609     if(!((fopt>= 1) &&(fopt<=5)))
 610         {fprintf(output, "\n         *** option exp-valueis out of range ***");
 611           xalg8_tmp = TRUE;
 612       return xalg8_tmp;}
 613
 614     if(!((gopt>=0 ) &&(gopt<=2)))
 615         {fprintf(output, "\n         *** option split is out of range ***");
 616           xalg8_tmp = TRUE;
 617       return xalg8_tmp;}
 618
 619     if(!((hopt>=0 ) &&(hopt<=1)))
 620         {fprintf(output, "\n         *** option normalize is out of range ***");
 621           xalg8_tmp = TRUE;
 622       return xalg8_tmp;}
 623
 624     if(!((popt>=0 ) &&(popt<=2)))
 625         {fprintf(output, "\n         *** option print is out of range ***");
 626           xalg8_tmp = TRUE;
 627       return xalg8_tmp;}
 628
 629     if(!(error_1> Consts.zero ))
 630         {fprintf(output, "\n         *** option error is out of range ***");
 631           xalg8_tmp = TRUE;
 632       return xalg8_tmp;}
 633
 634     if(!(Err.minprd> Consts.zero ))
 635         {fprintf(output, "\n         *** option minv is out of range ***");
 636           xalg8_tmp = TRUE;
 637       return xalg8_tmp;}
 638
 639     if(!((Fourie.interp>= -1) &&(Fourie.interp<=6)))
 640         {fprintf(output, "\n         *** option interp is out of range ***");
 641           xalg8_tmp = TRUE;
 642       return xalg8_tmp;}
 643
 644
 645     ir = 0;
 646     if(((aopt == 3) || (aopt == 4)) && ((bopt == 1) || (bopt == -1))) ir = 3;   // corr. and bug in Fortran hstau
 647     if((aopt < 1) || (aopt > 4)) ir = 1;
 648     if((bopt < -2) || (bopt > 3)) ir += 2;  //bug in Fortran hstau
 649     if((copt < 1) || (copt > 4)) ir += 4;
 650
 651     if(!((aopt > 2) || (abs(bopt))  !=  2)) {
 652
 653       if((dopt < 1) || (dopt > 3)) ir += 8;
 654     }
 655
 656     if((eopt < 1) || (eopt > 3)) ir += 16;
 657     if(copt  >  3)  {
 658
 659       if((fopt < 1) || (fopt > 5)) ir += 32;//bug in Fortran hstau
 660     }
 661     if((gopt < 0) || (gopt > 2)) ir += 64;
 662     if((hopt < 0) || (hopt > 1)) ir += 128;
 663     if(period < 1)               ir += 256;  // added hstau
 664     if(toler < Consts.zero)      ir += 512;  // added hstau
 665     if( popt < 0  ||  popt > 2 ) ir += 1024;  // added hstau
 666     if(error_1 < Consts.zero)    ir += 2048;  // added hstau
 667
 668     if(eopt == 2) {
 669       if(Err.minprd < Consts.zero) ir += 4096;  // added hstau
 670     }
 671
 672     if(fopt == 4) {
 673       if(Fourie.interp < -1 ||  Fourie.interp > 6 ) ir += 8192; // added hstau
 674     }
 675
 676     if(aopt > 2 && abs(bopt) == 2) {  //eleminate this quad algorithm from the old manual. hstau 1/04
 677       fprintf(output, "\n\n      ***delta_flag cannot be 2 or -2 for algorithm = 3 or 4***");
 678       xalg8_tmp = TRUE;
 679       return xalg8_tmp;
 680     }
 681
 682       // OPTION OUT OF RANGE
 683     if(ir > 0) {
 684       for(i = 0; i < 14; i++) {
 685         if(ir != ((ir/2)*2)) {
 686           //fprintf(output, "\n\n      *** %c opt out of range, program stops***", let[i]);
 687           fprintf(output, "\n\n      ***option %s out of range, program stops***", let[i]); //hstau
 688         }
 689
 690         ir = ir/2;
 691       }
 692       xalg8_tmp = TRUE;
 693       return xalg8_tmp;
 694
 695     }
 696
 697     // PRINT HEADING
 698
 699     if((aopt >= 3) && (bopt == 3)) {
 700       fprintf(output, "\n           the conjugate gradient algorithm");
 701     }
 702
 703     if((aopt <= 2) && (bopt == 3)) {
 704       fprintf(output, "\n           the steepest descent algorithm");
 705     }
 706
 707     if((aopt >= 3) && (abs(bopt) == 2)) {
 708       fprintf(output, "\n           richardson semi iterative method");
 709     }
 710
 711     if((aopt <= 2) && (abs(bopt) == 2)) {
 712       fprintf(output, "\n           richardson one step method");
 713     }
 714
 715     if((aopt >= 3) && (abs(bopt)  <=  1)) {
 716       fprintf(output, "\n           a stationary two steps method");
 717     }
 718
 719     if((aopt <= 2) && (abs(bopt)  <=  1)) {
 720       fprintf(output, "\n           a stationary one steps method");
 721     }
 722
 723     headng = hbb;
 724     if(aopt == (aopt/2)*2) headng = hin;
 725     fprintf(output, "\n\n            %2sdirect gradient computation", int2str(headng)); //corr. reduce space betw %2s and direct hstau
 726     fprintf(output, "\n          =================================");
 727
 728     if(gopt == 1) {
 729       fprintf(output, "\n\n              diagonal scaling is used");
 730     }
 731
 732     if(gopt == 2) {
 733       fprintf(output, "\n\n              with user supplied scalling");
 734     }
 735
 736     //bug 154 this checking is ignored hstau. 8 2005
 737     /*
 738     // CHECK RATIO OF PIXSIZ AND RINC (APPLY ONLY FOR STRIP MODE) affect blobs? hstau
 739
 740     ratio = GeoPar.pixsiz/GeoPar.pinc;
 741     iratio = (INTEGER)(ratio + 0.5);
 742     if((GeoPar.strip) && ((iratio < 1) || (fabs(iratio-ratio) >= 0.05))) {
 743       fprintf(output, "\n\n      ***warning pinc is not an integer multiplication of pixsiz");
 744     }
 745     */
 746
 747     // READ USER SUPLIED ORDER IF CALLED FOR
 748
 749     if(!((aopt > 2) || (abs(bopt) != 2))) {
 750       order = new INTEGER[period];
 751
 752       if(dopt == 2) {
 753
 754         order[0] = InFile.getint(TRUE, &eol);
 755
 756       if(eol || order[0] < 1 || order[0] > period){
 757         fprintf(output, "\n\n       insufficient data item, user specified order required \n   or illegal order number, order should be integer in [1, period]");
 758         xalg8_tmp = TRUE;
 759         return xalg8_tmp;}
 760
 761         for(i = 1; i < period; i++) {   //corr.hstau
 762           //alg = eol;
 763           ///if(alg) order[i] = getint(input, &eol);
 764           //if(alg) order[i] = InFile.getint(TRUE, &eol);
 765           ///if(!alg) order[i] = getint(NULL, &eol);
 766           //if(!alg) order[i] = InFile.getint(FALSE, &eol);
 767            order[i] = InFile.getint(FALSE, &eol);
 768
 769       if(eol || order[0] < 1 || order[0] > period){
 770         fprintf(output, "\n\n       insufficient data item, user specified order required \n   or illegal order number, order should be integer in [1, period] \n");
 771         xalg8_tmp = TRUE;
 772         return xalg8_tmp;}                                //wei, bug 174, 10/2005
 773         }
 774         //alg = FALSE;
 775       }
 776     }
 777
 778     // FOR  COPT .GT. 3 READ NEXT CARDS
 779
 780     if(copt <= 3) {
 781
 782       // SET UNUSED VARIABLES
 783
 784       sa = 1.0;
 785       sb = 0.0;
 786       sc = 0.0;
 787       l = 0;
 788       t = 0;
 789     }
 790     else {
 791       ///sa = getnum(input, &eol);
 792       ///sb = getnum(NULL, &eol);
 793       ///sc = getnum(NULL, &eol);
 794       ///l = getint(input, &eol);
 795       ///bw1 = getnum(NULL, &eol);
 796       ///bw2 = getnum(NULL, &eol);
 797       ///bw3 = getnum(NULL, &eol);
 798       ///bbcon = getwrd(NULL, &eol);
 799       ///t = getint(input, &eol);
 800       ///cw1 = getnum(NULL, &eol);
 801       ///cw2 = getnum(NULL, &eol);
 802       ///cw3 = getnum(NULL, &eol);
 803       ///cbcon = getwrd(NULL, &eol);
 804
 805       sa = InFile.getnum(TRUE, &eol);
 806       sb = InFile.getnum(FALSE, &eol);
 807       sc = InFile.getnum(FALSE, &eol);
 808
 809       if(eol){
 810         fprintf(output, "\n\n       insufficient data item, a, b, c missing, check the manual \n");
 811         xalg8_tmp = TRUE;
 812         return xalg8_tmp;}
 813
 814       l = InFile.getint(TRUE, &eol);
 815       bw1 = InFile.getnum(FALSE, &eol);
 816       bw2 = InFile.getnum(FALSE, &eol);
 817       bw3 = InFile.getnum(FALSE, &eol);
 818       bbcon = InFile.getwrd(FALSE, &eol, quad_blurck_codes, 4);
 819
 820      if(eol){
 821         fprintf(output, "\n\n       insufficient data item, matrix B missing, check the manual \n");
 822         xalg8_tmp = TRUE;
 823         return xalg8_tmp;}
 824
 825       t = InFile.getint(TRUE, &eol);
 826       cw1 = InFile.getnum(FALSE, &eol);
 827       cw2 = InFile.getnum(FALSE, &eol);
 828       cw3 = InFile.getnum(FALSE, &eol);
 829       cbcon = InFile.getwrd(FALSE, &eol, quad_blurck_codes, 4);
 830
 831       if(eol){
 832         fprintf(output, "\n\n       insufficient data item \n, matrix C missing, check the manual \n");
 833         xalg8_tmp = TRUE;
 834         return xalg8_tmp;}
 835
 836       //tracer
 837       //      fprintf(output,"t=%d", t);
 838       //tracer
 839
 840       // CHECK INPUT
 841
 842       if(t == 0) copt = 5;  // should it not be a 4?? hstau
 843       fprintf(output, "\n\n      the weights are: sa= %8.4f sb= %8.4f sc= %8.4f", sa, sb, sc);
 844
 845       if(!((sa >= 0.0) && (sb >= 0.0) && (sc >= 0.0))) { // there was a bug in this line (see Fortran version). hstau 6/2/2003
 846         fprintf(output, "\n\n      ***negative weight, program stops***");
 847
 848         xalg8_tmp = TRUE;
 849       }
 850
 851       if(!((sb == 0.0) || (l == 0))) {
 852         fprintf(output, "\n\n      matrix b is raised to the power %4i", l);
 853
 854
 855         if(l <= 0) {
 856           fprintf(output, "\n\n      ***negative power, program stops***");
 857           xalg8_tmp = TRUE;
 858         }
 859         blurck(bw1, bw2, bw3, bbcon, l, &alg);
 860
 861         xalg8_tmp = alg;
 862         if((Blob.pix_basis && GeoPar.nelem < 3) || (!Blob.pix_basis && Blob.H < 5) ) {
 863           fprintf(output, "\n\n      ***nelem (or the number of blobs in one direction) must be greater than 2 (4) for smoothing");
 864           fprintf(output, "\n      program stops***");
 865           xalg8_tmp = TRUE;
 866         }
 867       }
 868       k = t;
 869       if(copt != 5) {
 870
 871         fprintf(output, "\n\n      t value for the matrix c is %4i", t);
 872
 873         if(t <= 0) {
 874           //fprintf(output, "\n\n      ***t value %4i is not allowed*** program stops", t);
 875           fprintf(output, "\n\n      ***t value %4i is not allowed, program stops***", t); //hstau
 876
 877           xalg8_tmp = TRUE;
 878         }
 879         else {
 880           if(!((aopt == 1) && (bopt  <=  2) && (gopt == 0))) {
 881
 882             // SET T ACORDING TO METHOD
 883
 884             t = k / 2;
 885             if(k != t*2) {
 886               //fprintf(output, "\n\n      ***with either aopt .gt. 1 .or. bopt .gt. 2 t .or. splitting must be even, program stops***");  //Fortran code like output?? hstau
 887               fprintf(output, "\n\n      ***with either aopt > 1 or bopt > 2, t or splitting must be even, program stops***");  //Fortran code like output?? hstau
 888
 889               xalg8_tmp = TRUE;
 890               goto L7;
 891             }
 892           }
 893           blurck(cw1, cw2, cw3, cbcon, t, &alg);
 894           xalg8_tmp = alg;
 895           if((Blob.pix_basis && GeoPar.nelem < 3) || (!Blob.pix_basis && Blob.H < 5) ) {
 896             fprintf(output, "\n\n      ***nelem (or the number of blobs in one direction) must be greater than 2 (4) for smoothing");
 897             fprintf(output, "\n      program stops***");
 898
 899             xalg8_tmp = TRUE;
 900           }
 901         }
 902       }
 903     }
 904     // ERROR CHECK
 905
 906 L7:
 907     if((eopt - 2) <= 0) {
 908       if((eopt - 2) < 0) {
 909
 910        //bug 154 this report is removed hstau. 8 2005
 911         /*
 912         fprintf(output, "\n\n      error is identity matrix times the constant %9.4f", error_1);
 913         */
 914
 915         if(fabs(error_1) > Consts.zero) goto L11;
 916         fprintf(output, "\n\n      ***abs(error) must be greater than ZERO*** program stops");
 917
 918         xalg8_tmp = TRUE;
 919         goto L11;
 920       }
 921       else {
 922         fprintf(output, "\n\n      error is calculated by program");
 923
 924         fprintf(output, "\n\n      minimum measurement is %9.4f", Err.minprd);
 925
 926         if(Err.minprd <= Consts.zero) {
 927           xalg8_tmp = TRUE;
 928           //fprintf(output, "\n\n      *** min variance is less than zero, program stops***");
 929           fprintf(output, "\n\n      *** minv is less than ZERO, program stops***");  //corr. hstau
 930
 931           goto L12;
 932         }
 933         errpar(alg);
 934          if(Err.errtyp == 0) error_1 = Err.minprd;
 935          if(Err.errtyp == 0) eopt = 1;  //suspicious logic. hstau
 936         xalg8_tmp = alg;
 937
 938         goto L11;
 939       }
 940     }
 941
 942     fprintf(output, "\n\n      user error routine");
 943
 944     dumy = uerror(1, 1, dum, &alg);
 945     xalg8_tmp = alg;
 946
 947
 948
 949 L11:
 950     if(!xalg8_tmp) {
 951       if(copt <= 3) {
 952
 953          //bug 154 this report is removed hstau. 8 2005
 954
 955         /*
 956         fprintf(output, "\n\n      second criterion matrix is %4s", int2str(dmart[2*copt-1]));
 957         //fprintf(output, "%4s matrix", int2str(dmart[2*copt-1]));
 958         fprintf(output, "%4s matrix", int2str(dmart[2*copt])); //corr.hstau
 959         */
 960
 961         d_ = new REAL[area];
 962
 963         if(copt  <=  2) deset(d_, list, weight, &alg);  // check this.hstau ok
 964         if(copt == 3) dset(d_, list, weight, &alg);
 965         xalg8_tmp = alg;
 966
 967         if(Blob.pix_basis) {
 968           prpict(GeoPar.nelem, d_, 1, "dset", "d");
 969         }
 970         else {
 971           bprpict(Blob.M, Blob.H, d_, 1, "dset", "d");
 972         }
 973       }
 974       // SPLIT IF CALLED FOR
 975
 976       if(!xalg8_tmp) {
 977         if(gopt > 0) {
 978           s_ = new REAL[area];
 979         }
 980         if(gopt == 1) seset(s_, list, weight, &alg);
 981         if(gopt == 2) sset(s_, list, weight, &alg);
 982         xalg8_tmp = alg;
 983         if(gopt > 0) {
 984           if(Blob.pix_basis) {
 985             prpict(GeoPar.nelem, s_, 1, "set ", "s");
 986           }
 987           else {
 988             bprpict(Blob.M, Blob.H, s_, 1, "set", "s");
 989           }
 990         }
 991         // FOR COPT .LE. 3 AND SPLITTING COMBINE THE VECTORS
 992
 993         if(!((copt > 3) || (gopt == 0))) {
 994           //js = s - 1;
 995           for(i = 0; i < area; i++) {
 996             d_[i] *= s_[i];
 997           }
 998           delete[] s_;  // bug 92 - Lajos - 03/02/2005
 999           s_ = NULL;    //bug 174, wei, 10/2005
1000           ///  s = 1; !!!              //check this.hstau
1001         }
1002       }
1003     }
1004 L12:
1005
1006     if(bopt != 3) {
1007       if((bopt == 0) && (aopt >= 3)) {
1008         fprintf(output, "\n\n      user supplied delta is %9.4f", delta);
1009         fprintf(output, "\n      user supplied gma   is %9.4f", gma);
1010       }
1011       if((bopt == 0) && (aopt <= 2)) {
1012         fprintf(output, "\n\n      user supplied delta is %9.4f", delta);
1013       }
1014       if(bopt != 0) {
1015         if(bopt <= 0) {
1016           //fprintf(output, "\n\n      user supplied smallest eigenvalue %14.6e", delta);
1017           fprintf(output, "\n\n      user supplied smallest eigenvalue %14.6e", gma);  // bug in Fortran .hstau
1018           //fprintf(output, "\n      user supplied largest eigenvalue  %14.6e", gma);
1019           fprintf(output, "\n      user supplied largest eigenvalue  %14.6e", delta);   //bug in Fortran .hstau
1020           if((delta+gma) <= Consts.zero) {
1021             fprintf(output, "\n\n      ***largest eigenvalue must be positive number");
1022             fprintf(output, "\n      greater than smallest eigenvalue, program stops***");
1023
1024             xalg8_tmp = TRUE;
1025             goto L13;
1026           }
1027
1028           //esmin = delta;
1029           //lmd = gma;
1030           esmin = gma; // bug in Fortran.hstau
1031           lmd = delta;  // bug in Fortran. hstau
1032           bopt = - bopt;
1033         }
1034         else {
1035           if(minz == 0) {
1036             fprintf(output, "\n\n      smallest eigenvalue is set to 0");
1037           }
1038           fprintf(output, "\n\n      eigenvalue is calculated to tolerance %9.4f", toler);
1039
1040           if(toler <= Consts.zero) {
1041             fprintf(output, "\n\n      ***tolerance too small, program stops***");
1042
1043             xalg8_tmp = TRUE;
1044             goto L13;
1045           }
1046
1047           if(xalg8_tmp) {
1048
1049             goto L13;
1050           }
1051
1052           r_ = new REAL[area];
1053           //eigval(r, (minz == 0), d, list, weight, toler, esmin, lmd, s);
1054           eigval(r_, (minz == 0), d_, list, weight, toler, &esmin, &lmd, s_); //corr. hstau
1055           delete[] r_;  // bug 92 - Lajos - 03/02/2005
1056           r_ = NULL; //bug 147, wei, 2005/10
1057
1058           if((esmin+lmd) <= Consts.zero) {
1059             xalg8_tmp = TRUE;
1060             fprintf(output, "\n\n      ***largest eigenvalue must be positive number");
1061             fprintf(output, "\n\n      greater than smallest eigenvalue, program stops***");
1062
1063             goto L13;
1064           }
1065         }
1066  //tracer
1067         //      fprintf(output,"\n quad: bopt=%d esmin=%f lmd=%f ",bopt,esmin, lmd);
1068
1069       //tracer
1070
1071         if(bopt == 1) delta = ((REAL) 2.0) / (esmin + lmd);
1072         if(bopt == 1) {
1073           fprintf(output, "\n\n      delta computed by program is %14.4e", delta);
1074         }
1075         if((bopt == 2) && (aopt  <=  2)) {
1076           fprintf(output, "\n\n      the period for richardson method is %2i", period);
1077         }
1078         if(bopt == 2) prsemi(esmin, lmd, &alg, &delta);//affected by blobs? hstau. no
1079         xalg8_tmp = alg;
1080
1081       }
1082     }
1083
1084 L13:
1085     j = 4*(fopt - 1);
1086     if(copt > 3) {
1087       fprintf(output, "\n\n      expected value is %4s", int2str(init[j+1]));
1088       fprintf(output, "%4s", int2str(init[j+2]));  // corr. hstau
1089       fprintf(output, "%4s", int2str(init[j+3]));
1090       fprintf(output, "%4s", int2str(init[j+4]));
1091     }
1092
1093     if(Modefl.lofl) {
1094       fprintf(output, "\n\n      lower value of pixels is set to %9.4f", Modefl.lower);
1095     }
1096     if(Modefl.upfl) {
1097       fprintf(output, "\n\n      upper value of pixels is set to %9.4f", Modefl.upper);
1098     }
1099     if(hopt >= 1) {
1100       fprintf(output, "\n\n      average density is corrected after each iteration");
1101     }
1102
1103     // RETURN IF THERE ARE ERRORS IN INPUT
1104
1105     if(xalg8_tmp) {
1106
1107
1108       return xalg8_tmp;
1109     }
1110
1111     // INITILIZE EXPECTED VECTOR
1112
1113     xbar = 0.0;
1114     if(fopt == 2) xbar = GeoPar.aveden;
1115     if(fopt > 2) {
1116       f_ = new REAL[GeoPar.area];   // f always in pixel basis. hstau
1117       if(fopt == 4) {
1118         if(GeoPar.div) {
1119           fprintf(output, "\n\n      warning   continuous back projection is not intended for divergent data");
1120         }
1121         Fourie.rinc = GeoPar.pinc;
1122         r_ = new REAL[GeoPar.nrays];
1123
1124         for(i = 0; i < GeoPar.area; i++) {
1125           f_[i] = 0.0;
1126         }
1127
1128         for(np = 0; np < GeoPar.prjnum; np++) {
1129           ProjFile.ReadProj(np, r_, GeoPar.nrays);
1130           Anglst.getang(np, &theta, &sinth, &costh);
1131           if(GeoPar.vri) Fourie.rinc = (REAL) (GeoPar.pinc * MAX0(fabs(costh), fabs(sinth)));
1132           spac = GeoPar.pixsiz/Fourie.rinc;
1133           bckprj(f_, GeoPar.nelem, r_, GeoPar.nrays, sinth, costh, spac, Fourie.interp);
1134         }
1135
1136         delete[] r_;  // bug 92 - Lajos - 03/02/2005
1137         r_ = NULL;//bug 147, wei, 2005/10
1138
1139         for(i = 0; i < GeoPar.area; i++) {
1140           f_[i] /= ((REAL)(GeoPar.prjnum));
1141         }
1142       }
1143       else {
1144         if(fopt == 3) back(f_, list, weight);
1145         if(fopt == 5) {
1146            if(Blob.pix_basis) {
1147               for(i = 0; i < GeoPar.area; i++) {
1148                  f_[i] = recon[i];
1149               }
1150            }
1151            else {
1152               fb_ = new REAL[Blob.area];
1153               for(i = 0; i < Blob.area; i++) {
1154                  fb_[i] = recon[i];
1155               }
1156            }
1157         }
1158       }
1159     }
1160
1161     if(!Blob.pix_basis && fopt > 2 && fopt !=5) {
1162        fb_ = new REAL[Blob.area];
1163        //tracer
1164        //fprintf(output,"\nf=%f %f %f area=%d", f[0],f[1],f[2],Blob.area);fflush(output);
1165        //tracer
1166        Blob.pix2blob(f_,fb_);
1167        delete[] f_;  // bug 92 - Lajos - 03/02/2005
1168        f_ = NULL; //bug 147, wei, 2005/10
1169
1170        //f = fb;
1171     }
1172
1173     // ALLOCATE AREA
1174
1175     r_ = new REAL[area];
1176     g_ = r_;
1177     v_ = r_;
1178     w_ = r_;
1179     if((aopt == 2) || (aopt == 4)) {
1180       v_ = new REAL[area];
1181     }
1182     if(aopt >= 3) {
1183       g_ = new REAL[area];
1184     }
1185
1186     cg_ = g_;
1187 //tracer
1188     //fprintf(output,"\n quad: cg=%f",cg[0]);
1189         ///tracer
1190
1191     if((sc != 0.0) || ((sb !=  0.0) && (l > 0))) {
1192       w_ = new REAL[area];
1193     }
1194
1195
1196     if((((copt !=  5) || (gopt > 0))  &&   // corr. why so many () hstau
1197        ((aopt !=  1) || (bopt > 2))) ||
1198        ((aopt == 1) && (bopt > 2))
1199        ) {
1200
1201       cg_ = new REAL[area];
1202     }
1203     // FIRST TIME AROUND CALL RESEDU FOR ALL METHODS
1204
1205   }  // if iter == 1
1206   else {
1207
1208   // FLOW CONTROL
1209
1210
1211   // AOPT EQUAL 2 OR AOPT EQUAL 4
1212
1213     if(aopt == (aopt/2)*2) {
1214       nextr(r_, v_, delta, area);
1215       if(Blob.pix_basis) {
1216         prpict(GeoPar.nelem, r_, iter, "next", "r");
1217       }
1218       else {
1219         bprpict(Blob.M, Blob.H, r_, iter, "next", "r");
1220       }
1221       goto L22;
1222     }
1223   }
1224
1225   // AOPT EQUAL 1 OR AOPT EQUAL 3
1226
1227   if(Blob.pix_basis) {
1228     resedu(recon, r_, w_, d_, list, weight, f_, s_);
1229 //tracer
1230     //fprintf(output,"\n quad2: r=%f, cg=%f",r[0],cg[0]);
1231         ///tracer
1232
1233     prpict(GeoPar.nelem, r_, iter, "rese", "r");
1234   }
1235   else {
1236     resedu(recon, r_, w_, d_, list, weight, fb_, s_);
1237
1238     bprpict(Blob.M, Blob.H, r_, iter, "rese", "r");
1239   }
1240
1241   // FOR BOPT = 2 GET DELTA AND GMA
1242
1243 L22:
1244   if(bopt == 2) semi(&delta, &gma, &alg);
1245
1246   xalg8_tmp = alg;
1247   if(xalg8_tmp) return xalg8_tmp;
1248
1249   // FOR CONSTANT DELTA AND AOPT EQUAL 1 CALCULATE NEXT PICTURE
1250
1251   if(!((bopt  <=  2) && (aopt == 1))) {
1252
1253     // SECOND CRITERIA
1254
1255     if(copt <= 3) {
1256
1257       //deltad(r, g, v, cg, d, delta, gma, list, weight, &alg);
1258       deltad(r_, g_, v_, cg_, d_, &delta, &gma, list, weight, &alg); //corr.hstau
1259
1260       if(Blob.pix_basis) {
1261         prpict(GeoPar.nelem, r_, iter, "deld", "r");
1262         prpict(GeoPar.nelem, cg_,iter, "deld", "cg");
1263       }
1264       else {
1265         bprpict(Blob.M, Blob.H, r_, iter, "deld", "r");
1266         bprpict(Blob.M, Blob.H, cg_,iter, "deld", "cg");
1267       }
1268       xalg8_tmp = alg;
1269       if(alg) return xalg8_tmp;
1270     }
1271     else {
1272
1273       //deltac(r, g, v, cg, w, delta, gma, list, weight, &alg, s);
1274       deltac(r_, g_, v_, cg_, w_, &delta, &gma, list, weight, &alg, s_);// corr.hstau
1275
1276
1277       if(Blob.pix_basis) {
1278         prpict(GeoPar.nelem, r_, iter, "delc", "r");
1279         prpict(GeoPar.nelem, g_, iter, "delc", "g");
1280         prpict(GeoPar.nelem, v_, iter, "delc", "v");
1281         prpict(GeoPar.nelem, cg_, iter, "delc", "cg");
1282       }
1283       else {
1284         bprpict(Blob.M, Blob.H, r_, iter, "delc", "r");
1285         bprpict(Blob.M, Blob.H, g_, iter, "delc", "g");
1286         bprpict(Blob.M, Blob.H, v_, iter, "delc", "v");
1287         bprpict(Blob.M, Blob.H, cg_, iter, "delc", "cg");
1288       }
1289
1290       xalg8_tmp = alg;
1291       if(alg) return xalg8_tmp;
1292     }
1293   }
1294
1295   // COMPUTE NEXT PICTURE
1296
1297   nxtpct(recon, cg_, delta);
1298
1299   // IF POPT .EQ. 2 AND POSSIBLE COMPUTE THE COST FUNCTION
1300
1301   if(popt >= 2) {
1302     if(Blob.pix_basis) {
1303       costfn(recon, r_, g_, cg_, w_, list, weight, f_);
1304     }
1305     else {
1306       costfn(recon, r_, g_, cg_, w_, list, weight, fb_);
1307     }
1308   }
1309   return xalg8_tmp;
1310 }
1311