2 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
6 * A PICTURE RECONSTRUCTION PROGRAM *
8 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
10 quad.cpp,v 1.6 2009/06/01 03:33:26 jklukowska Exp
12 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **
14 * QUADRATIC ITERATIVE RECONSTRUCTION ALGORITHMS *
16 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **
19 PROGRAM BY EHUD ARTZY 1977
24 THIS PROGRAM IS SET TO MINIMIZE THE FUNCTION
26 K(X)=SA*(P-M*X)#*A*(P-M*X)+(X-XBAR)#*(SB*B**L+SC*C**(-1)*(X-XBAR)
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
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)
43 IF SB + SC GREATER THAN 0 THEN
47 (1) Q = C**(T/2)*((SA*M#*A*M +SB*B**L)*C**(T/2) +(SC*C**(T/2 - 1))
49 F = C**(T/2)*SA*M#*A*(P - M*XBAR)
51 IF Y IS A SOLUTION OF Q*Y = F THEN X = C**(T/2)*Y + XBAR
60 IF Y IS A SOLUTION OF Q*Y = F THEN X = D*Y IS A
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)
68 THE FOLLOWING ALGORITHM IS USED TO FIND THE SOLUTION OF Q*Y = F
70 G(0) = R(0) = Q*Y(0) - F
71 Y(K+1) = Y(K) - DELTA(K)*G(K)
73 G(K+1) = R(K+1) + GMA(K)*G(K)
75 AN ALTERNATIVE WAY TO COMPUTE R(K+1) IS GIVEN BY
76 R(K+1) = R(K) - DELTA(K)*Q*G(K)
78 THE FORMER IS REFFERED TO AS DIRECT METHOD AND THE LATER INDIRECT
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)
85 FOR USING THE PROGRAM THE USER MUST SPECIFY THE FOLLOWING
88 CARD 1 - AOPT,BOPT,COPT,MINZ,PERIOD,TOLER,DELTA,GMA
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
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
100 DELTA = 2.0/(LMD + ESMIN)
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:
110 ITERATIVE SOLUTION OF LARGE LINEAR SYSTEMS
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))
119 COPT : 1 - EQUATION 2 IS USED WITH GOITEIN D MATRIX
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
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
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
140 LIST -SEE THE SNARK MANUAL
141 WEIGHT-SEE THE SNARK MANUAL
142 4 - EQUATION 1 IS USED CARD 3-5 MUST BE PRESENTED
144 MINZ : 0 - FOR BOPT 1 OR 2 ESMIN IS ASSUMED TO BE 0
145 1 - FOR BOPT 1 OR 2 ESMIN IS CALCULATED
147 PERIOD: N - N IS THE LENGTH OF THE PERIOD FOR THE RICHARDSON OSM
148 N MUST BE IN THE RANGE 1-32
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)
157 DELTA : - FOR BOPT 0 THIS IS CONSTANT DELTA
158 FOR BOPT -1,-2 THIS IS USER SUPPLY ESMIN
160 GMA : - FOR BOPT 0 AND TSS THIS IS CONSTANT GMA
161 FOR BOPT -1,-2 THIS IS USER SUPPLY LMD
163 CARD 2 - DOPT,EOPT,FOPT,GOPT,HOPT,POPT,ERROR,MINPRD,INTERP
165 DOPT : 1 - NATURAL ORDER FOR PICKING DELTA FOR RICHARDSON OSM
167 2 - USER SUPPLY ORDER OF PICKING DELTA CARD 6 MUST BE
169 3 - LEBEDTV-FINOGENOV ORDERING FOR PICKING DELTA
171 R.S. ANDERSSEN AND G.H. GOLUB
172 RICHARDSONS NON-STATIONARY MATRIX ITERATIVE PROCEDURE
173 STAN-CS-72-307 (1972)
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
180 3 - USER SUPPLY AN ERROR ROUTINE. THE CALLING SEQUENCE
182 PUNCTION UERROR(NP,NR,PD)
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
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)
198 S IS A DIAGONAL SPLITTING MATRIX
199 LIST SEE THE SNARK05 MANUAL
200 WEIGHT SEE THE SNARK05 MANUAL
202 HOPT : 1 - RECONSTRUCTED PICTURE IS CORRECTED SO THAT AVERAGE
203 DENSITY WILL BE EQUAL TO THAT OF AVERAGE DENSITY
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
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
218 INTERP: - FOR FOPT 4 THIS IS THE INTERPOLATION INDEX
219 (SEE THE SNARK05 MANUAL)
221 THE FOLLOWING CARD(S) IS(ARE) READ ONLY IF BOPT = 2
223 ORDER(I), I = 1,PERIOD
224 ORDER : PERIOD NUMBERS IN THE RANGE 1 TO PERIOD
226 THE FOLLOWING 3 CARDS ARE READ ONLY IF COPT GE 4
230 SA : CONSTANT (SEE EQUATION 1)
231 SB : CONSTANT (SEE EQUATION 1)
232 SC : CONSTANT (SEE EQUATION 1)
234 CARD 4 - L,BW1,BW2,BBCON *** MATRIX B ***
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)
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
251 CARD 5 - T,CW1,CW2,CW3,CBCON *** MATRIX C ***
253 T : T/2 IS THE POWER OF MATRIX C IN EQUATION 1
261 ***********************************************************************
278 #include "projfile.h"
286 INTEGER quad_class::Init()
301 del = NULL; //bug 174, wei, 10/2005
304 return 0; //wei, 3/2005
307 INTEGER quad_class::Reset()
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;
325 return 0; //wei, 3/2005
328 BOOLEAN quad_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
334 static REAL gma;//set static. hstau see deltac() below
335 //INTEGER popt; poot is global, why is it declared here? hstau
370 REAL dum = 0.0; // this is dummy variable
372 static REAL esmin; // needed to be set static? hstau
385 static REAL* cg_;//set static. hstau see deltad() below
386 static REAL* g_;//set static. hstau see deltad() below
388 static REAL* w_; //set static. hstau
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
396 static REAL lmd; // needed to be set static? hstau
397 ///dimension init(20),dmart(6),let(8)
399 static BOOLEAN alg;//set static. hstau
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',' ',' ',' ')
416 const char *let[] = {
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(' ',' ',' ',' ')
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 /
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',' ',' ')
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
474 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 16, 2004
475 static const INTEGER quad_blurck_codes[4] =
477 CHAR2INT('s','n','a','r'),
478 CHAR2INT('k','a','m','i'),
479 CHAR2INT('s','a','j','z'),
480 CHAR2INT('s','a','j','c')
498 // for(i = 0; i < area; i++) {
499 // fprintf(output,"\nrecon=%f", recon[i]);
503 // THE FOLLOWING SECTION IS EXECUTED ONLY AT FIRST ITERATION
507 // INITILIZE DUMMY VARIABLE
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);
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);
545 fprintf(output, "\n *** insufficient data item - end of line encountered ***");
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);
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);
563 fprintf(output, "\n *** insufficient data item - end of line encountered ***");
568 // CHECK RANGE OF OPTIONS
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
574 if(!((aopt>=1 ) &&(aopt<=4)))
575 { fprintf(output, "\n *** option algorithm is out of range ***");
579 if(!((bopt>= -2) &&(bopt<=3)))
580 {fprintf(output, "\n *** option delta-flga is out of range ***");
584 if(!((copt>= 1) &&(copt<=4)))
585 {fprintf(output, "\n *** option d-flag is out of range ***");
590 {fprintf(output, "\n *** option period is out of range ***");
594 if(!(toler> Consts.zero))
595 {fprintf(output, "\n *** option toler is out of range ***");
599 if(!((dopt>=1 ) &&(dopt<=3)))
600 {fprintf(output, "\n *** option order is out of range ***");
604 if(!((eopt>=1 ) &&(eopt<=3)))
605 {fprintf(output, "\n *** option error-type is out of range ***");
609 if(!((fopt>= 1) &&(fopt<=5)))
610 {fprintf(output, "\n *** option exp-valueis out of range ***");
614 if(!((gopt>=0 ) &&(gopt<=2)))
615 {fprintf(output, "\n *** option split is out of range ***");
619 if(!((hopt>=0 ) &&(hopt<=1)))
620 {fprintf(output, "\n *** option normalize is out of range ***");
624 if(!((popt>=0 ) &&(popt<=2)))
625 {fprintf(output, "\n *** option print is out of range ***");
629 if(!(error_1> Consts.zero ))
630 {fprintf(output, "\n *** option error is out of range ***");
634 if(!(Err.minprd> Consts.zero ))
635 {fprintf(output, "\n *** option minv is out of range ***");
639 if(!((Fourie.interp>= -1) &&(Fourie.interp<=6)))
640 {fprintf(output, "\n *** option interp is out of range ***");
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;
651 if(!((aopt > 2) || (abs(bopt)) != 2)) {
653 if((dopt < 1) || (dopt > 3)) ir += 8;
656 if((eopt < 1) || (eopt > 3)) ir += 16;
659 if((fopt < 1) || (fopt > 5)) ir += 32;//bug in Fortran hstau
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
669 if(Err.minprd < Consts.zero) ir += 4096; // added hstau
673 if(Fourie.interp < -1 || Fourie.interp > 6 ) ir += 8192; // added hstau
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***");
682 // OPTION OUT OF RANGE
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
699 if((aopt >= 3) && (bopt == 3)) {
700 fprintf(output, "\n the conjugate gradient algorithm");
703 if((aopt <= 2) && (bopt == 3)) {
704 fprintf(output, "\n the steepest descent algorithm");
707 if((aopt >= 3) && (abs(bopt) == 2)) {
708 fprintf(output, "\n richardson semi iterative method");
711 if((aopt <= 2) && (abs(bopt) == 2)) {
712 fprintf(output, "\n richardson one step method");
715 if((aopt >= 3) && (abs(bopt) <= 1)) {
716 fprintf(output, "\n a stationary two steps method");
719 if((aopt <= 2) && (abs(bopt) <= 1)) {
720 fprintf(output, "\n a stationary one steps method");
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 =================================");
729 fprintf(output, "\n\n diagonal scaling is used");
733 fprintf(output, "\n\n with user supplied scalling");
736 //bug 154 this checking is ignored hstau. 8 2005
738 // CHECK RATIO OF PIXSIZ AND RINC (APPLY ONLY FOR STRIP MODE) affect blobs? hstau
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");
747 // READ USER SUPLIED ORDER IF CALLED FOR
749 if(!((aopt > 2) || (abs(bopt) != 2))) {
750 order = new INTEGER[period];
754 order[0] = InFile.getint(TRUE, &eol);
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]");
761 for(i = 1; i < period; i++) { //corr.hstau
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);
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");
772 return xalg8_tmp;} //wei, bug 174, 10/2005
778 // FOR COPT .GT. 3 READ NEXT CARDS
782 // SET UNUSED VARIABLES
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);
805 sa = InFile.getnum(TRUE, &eol);
806 sb = InFile.getnum(FALSE, &eol);
807 sc = InFile.getnum(FALSE, &eol);
810 fprintf(output, "\n\n insufficient data item, a, b, c missing, check the manual \n");
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);
821 fprintf(output, "\n\n insufficient data item, matrix B missing, check the manual \n");
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);
832 fprintf(output, "\n\n insufficient data item \n, matrix C missing, check the manual \n");
837 // fprintf(output,"t=%d", t);
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);
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***");
851 if(!((sb == 0.0) || (l == 0))) {
852 fprintf(output, "\n\n matrix b is raised to the power %4i", l);
856 fprintf(output, "\n\n ***negative power, program stops***");
859 blurck(bw1, bw2, bw3, bbcon, l, &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***");
871 fprintf(output, "\n\n t value for the matrix c is %4i", t);
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
880 if(!((aopt == 1) && (bopt <= 2) && (gopt == 0))) {
882 // SET T ACORDING TO METHOD
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
893 blurck(cw1, cw2, cw3, cbcon, t, &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***");
907 if((eopt - 2) <= 0) {
910 //bug 154 this report is removed hstau. 8 2005
912 fprintf(output, "\n\n error is identity matrix times the constant %9.4f", error_1);
915 if(fabs(error_1) > Consts.zero) goto L11;
916 fprintf(output, "\n\n ***abs(error) must be greater than ZERO*** program stops");
922 fprintf(output, "\n\n error is calculated by program");
924 fprintf(output, "\n\n minimum measurement is %9.4f", Err.minprd);
926 if(Err.minprd <= Consts.zero) {
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
934 if(Err.errtyp == 0) error_1 = Err.minprd;
935 if(Err.errtyp == 0) eopt = 1; //suspicious logic. hstau
942 fprintf(output, "\n\n user error routine");
944 dumy = uerror(1, 1, dum, &alg);
953 //bug 154 this report is removed hstau. 8 2005
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
963 if(copt <= 2) deset(d_, list, weight, &alg); // check this.hstau ok
964 if(copt == 3) dset(d_, list, weight, &alg);
968 prpict(GeoPar.nelem, d_, 1, "dset", "d");
971 bprpict(Blob.M, Blob.H, d_, 1, "dset", "d");
974 // SPLIT IF CALLED FOR
980 if(gopt == 1) seset(s_, list, weight, &alg);
981 if(gopt == 2) sset(s_, list, weight, &alg);
985 prpict(GeoPar.nelem, s_, 1, "set ", "s");
988 bprpict(Blob.M, Blob.H, s_, 1, "set", "s");
991 // FOR COPT .LE. 3 AND SPLITTING COMBINE THE VECTORS
993 if(!((copt > 3) || (gopt == 0))) {
995 for(i = 0; i < area; i++) {
998 delete[] s_; // bug 92 - Lajos - 03/02/2005
999 s_ = NULL; //bug 174, wei, 10/2005
1000 /// s = 1; !!! //check this.hstau
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);
1011 if((bopt == 0) && (aopt <= 2)) {
1012 fprintf(output, "\n\n user supplied delta is %9.4f", delta);
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***");
1030 esmin = gma; // bug in Fortran.hstau
1031 lmd = delta; // bug in Fortran. hstau
1036 fprintf(output, "\n\n smallest eigenvalue is set to 0");
1038 fprintf(output, "\n\n eigenvalue is calculated to tolerance %9.4f", toler);
1040 if(toler <= Consts.zero) {
1041 fprintf(output, "\n\n ***tolerance too small, program stops***");
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
1058 if((esmin+lmd) <= Consts.zero) {
1060 fprintf(output, "\n\n ***largest eigenvalue must be positive number");
1061 fprintf(output, "\n\n greater than smallest eigenvalue, program stops***");
1067 // fprintf(output,"\n quad: bopt=%d esmin=%f lmd=%f ",bopt,esmin, lmd);
1071 if(bopt == 1) delta = ((REAL) 2.0) / (esmin + lmd);
1073 fprintf(output, "\n\n delta computed by program is %14.4e", delta);
1075 if((bopt == 2) && (aopt <= 2)) {
1076 fprintf(output, "\n\n the period for richardson method is %2i", period);
1078 if(bopt == 2) prsemi(esmin, lmd, &alg, &delta);//affected by blobs? hstau. no
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]));
1094 fprintf(output, "\n\n lower value of pixels is set to %9.4f", Modefl.lower);
1097 fprintf(output, "\n\n upper value of pixels is set to %9.4f", Modefl.upper);
1100 fprintf(output, "\n\n average density is corrected after each iteration");
1103 // RETURN IF THERE ARE ERRORS IN INPUT
1111 // INITILIZE EXPECTED VECTOR
1114 if(fopt == 2) xbar = GeoPar.aveden;
1116 f_ = new REAL[GeoPar.area]; // f always in pixel basis. hstau
1119 fprintf(output, "\n\n warning continuous back projection is not intended for divergent data");
1121 Fourie.rinc = GeoPar.pinc;
1122 r_ = new REAL[GeoPar.nrays];
1124 for(i = 0; i < GeoPar.area; i++) {
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);
1136 delete[] r_; // bug 92 - Lajos - 03/02/2005
1137 r_ = NULL;//bug 147, wei, 2005/10
1139 for(i = 0; i < GeoPar.area; i++) {
1140 f_[i] /= ((REAL)(GeoPar.prjnum));
1144 if(fopt == 3) back(f_, list, weight);
1146 if(Blob.pix_basis) {
1147 for(i = 0; i < GeoPar.area; i++) {
1152 fb_ = new REAL[Blob.area];
1153 for(i = 0; i < Blob.area; i++) {
1161 if(!Blob.pix_basis && fopt > 2 && fopt !=5) {
1162 fb_ = new REAL[Blob.area];
1164 //fprintf(output,"\nf=%f %f %f area=%d", f[0],f[1],f[2],Blob.area);fflush(output);
1166 Blob.pix2blob(f_,fb_);
1167 delete[] f_; // bug 92 - Lajos - 03/02/2005
1168 f_ = NULL; //bug 147, wei, 2005/10
1175 r_ = new REAL[area];
1179 if((aopt == 2) || (aopt == 4)) {
1180 v_ = new REAL[area];
1183 g_ = new REAL[area];
1188 //fprintf(output,"\n quad: cg=%f",cg[0]);
1191 if((sc != 0.0) || ((sb != 0.0) && (l > 0))) {
1192 w_ = new REAL[area];
1196 if((((copt != 5) || (gopt > 0)) && // corr. why so many () hstau
1197 ((aopt != 1) || (bopt > 2))) ||
1198 ((aopt == 1) && (bopt > 2))
1201 cg_ = new REAL[area];
1203 // FIRST TIME AROUND CALL RESEDU FOR ALL METHODS
1211 // AOPT EQUAL 2 OR AOPT EQUAL 4
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");
1219 bprpict(Blob.M, Blob.H, r_, iter, "next", "r");
1225 // AOPT EQUAL 1 OR AOPT EQUAL 3
1227 if(Blob.pix_basis) {
1228 resedu(recon, r_, w_, d_, list, weight, f_, s_);
1230 //fprintf(output,"\n quad2: r=%f, cg=%f",r[0],cg[0]);
1233 prpict(GeoPar.nelem, r_, iter, "rese", "r");
1236 resedu(recon, r_, w_, d_, list, weight, fb_, s_);
1238 bprpict(Blob.M, Blob.H, r_, iter, "rese", "r");
1241 // FOR BOPT = 2 GET DELTA AND GMA
1244 if(bopt == 2) semi(&delta, &gma, &alg);
1247 if(xalg8_tmp) return xalg8_tmp;
1249 // FOR CONSTANT DELTA AND AOPT EQUAL 1 CALCULATE NEXT PICTURE
1251 if(!((bopt <= 2) && (aopt == 1))) {
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
1260 if(Blob.pix_basis) {
1261 prpict(GeoPar.nelem, r_, iter, "deld", "r");
1262 prpict(GeoPar.nelem, cg_,iter, "deld", "cg");
1265 bprpict(Blob.M, Blob.H, r_, iter, "deld", "r");
1266 bprpict(Blob.M, Blob.H, cg_,iter, "deld", "cg");
1269 if(alg) return xalg8_tmp;
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
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");
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");
1291 if(alg) return xalg8_tmp;
1295 // COMPUTE NEXT PICTURE
1297 nxtpct(recon, cg_, delta);
1299 // IF POPT .EQ. 2 AND POSSIBLE COMPUTE THE COST FUNCTION
1302 if(Blob.pix_basis) {
1303 costfn(recon, r_, g_, cg_, w_, list, weight, f_);
1306 costfn(recon, r_, g_, cg_, w_, list, weight, fb_);