src/snark/dcon.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/dcon.cpp $
   5 $LastChangedRevision: 85 $
   6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7 $Author: agulati $
   8 **********************************************************************
   9
  10  PURPOSE :
  11      THIS IS A CONVOLUTION TYPE RECONSTRUCTION ALGORITHM FOR
  12      DIVERGENT BEAM PROJECTION DATA.  VARIOUS KINDS OF REGULARIZING
  13      FUNCTIONS ARE IMPLEMENTED.
  14  INPUT FILE :
  15      TWO (OR MORE) DATA CARDS HAVE TO BE SUPPLIED BY THE USER
  16      FOLLOWING EXECUTION COMMAND. THE INFORMATION ON AND AFTER
  17      THE THIRD CARD ARE FILTER PARAMETERS, THEREFORE THEY ARE
  18      FILTER DEPENDENT.
  19       CARD            LIST
  20       ====            ====
  21        1      INTERP,MISSP,NCONR,WCONV,METHOD
  22        2      NF
  23        3      C,ALPHA          -- IF(NF.NE.4HLINE
  24                                   .AND.NF.NE.4HBAND)
  25        3      X1,Y1,X2,Y2,...,XL,YL
  26                                -- IF(NF.EQ.4HLINE)
  27      DESCRIPTION OF INPUT PARAMETERS :
  28        INTERP - PARAMETER WHICH SPECIFIES INTERPOLATION METHOD,
  29                 THE VALUE OF INTERP MUST LIE BETWEEN -1 AND 6 AND
  30                 HAS THE FOLLOWING MEANING:
  31                 -1 -- MODIFIED CUBIC SPLINE INTERPOLATION,
  32                  0 -- BAND-LIMITING (SINC) INTERPOLATION,
  33                 POSITIVE -- INTERP-POINTS LAGRANGE INTERPOLATION.
  34        MISSP  - EVERY MISSP'TH PROJECTION IS TO BE USED IN
  35                 RECONSTRUCTION.
  36        NCONR  - NUMBER OF RAYS USED ON EACH SIDE OF CONVOLUTION,
  37                 FULL CONVOLUTION USED IF NCONR IS NEGATIVE.
  38        WCONV  - WEIGHT FOR CENTER RAY IN PRE-SMOOTHING.
  39        METHOD - PARAMETER WHICH DETERMINES THE METHOD IN FINDING
  40                 THE RAY NUMBER AND WEIGHT FOR BACK PROJECTION.
  41                  0 -- ACCURATE BACK PROJECTION,
  42                  1 -- WEIGHT FOR BACKPROJECTION APPROXIMATED,
  43                 -1 -- BOTH WIEGHT AND RAY NUMBER APPROXIMATED.
  44        NF     - FILTER NAME, CAN BE ONE OF LINE, HAMMING, SINC,
  45                 EXPONENTIAL, PARABOLIC, COSINE, SHEPP, OR BAND.
  46        C      - RATIO OF FILTER BANDWIDTH TO NYQUIST CUTOFF FREQ.
  47        ALPHA  - PARAMETER USED IN HAMMING FILTER ONLY.
  48        (XN,YN)- COORDINATES WHICH DEFINE THE FILTER FUNCTION
  49                   SEE SNARK05 MANUAL FOR DETAIL
  50      EXAMPLES OF CONTROL CARD SEQUENCE :
  51        1)     4    1   -1  1.0    0
  52           HAMMING
  53                  1.0      0.54
  54        2)     2    3  100 0.75    1
  55           LINE
  56                  0.3       1.0       0.5       0.7
  57                  0.8       0.4       1.0       0.0
  58        3)     0    1   -1  0.9   -1
  59           BAND-LIMITING
  60  WARNING :
  61      IT IS ASSUMED THAT PROJECTION ANGLES ARE GIVEN IN INCREASING
  62      ORDER AND THE LAST ANGLE IS LESS THAN THE FIRST ANGLE PLUS
  63      TWOPI.
  64      IT IS THE RESPONSBILITY OF THE USER TO CHECK THAT USRAYS
  65      IS LARGE ENOUGH FOR THE INTERPOLATION METHOD TO WORK.
  66  REFERENCES :
  67      1) FAST IMAGE RECONSTRUCTION BASED ON A RADON INVERSION
  68         FORMULA APPROPRIATE FOR RAPIDLY COLLECTED DATA
  69         G. T. HERMAN AND A. NAPARSTEK
  70         SIAM JOURNAL ON APPLIED MATHEMATICS V.33 (1977)
  71      2) FILTER SELECTION FOR FAN-BEAM CONVOLUTION RECONSTRUCTION
  72         ALGORITHM
  73         T. CHANG AND G. T. HERMAN
  74         J. COMPUTER ASSISTED TOMOGRAPHY, TO APPEAR
  75  WRITTEN BY T. CHANG, MAY 1977
  76  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  77 */
  78
  79 #include <cmath>
  80 #include <cstdio>
  81
  82 #include "blkdta.h"
  83 #include "geom.h"
  84 #include "spctrm.h"
  85 #include "consts.h"
  86 #include "uiod.h"
  87 #include "anglst.h"
  88 #include "qintp.h"
  89 #include "second.h"
  90 #include "projfile.h"
  91 #include "infile.h"
  92
  93 #include "dcon.h"
  94
  95 BOOLEAN dcon_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  96 {
  97   //fprintf(output,"\n xalg4");
  98
  99   static const INTEGER band = CHAR2INT('b','a','n','d');
 100
 101   // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
 102   // added missing options from "dcon_dconft.cpp". Lajos, Jan 13, 2005
 103   static const INTEGER dcon_codes[8] =
 104   {
 105     CHAR2INT('b','a','n','d'),
 106     CHAR2INT('l','i','n','e'),
 107     CHAR2INT('h','a','m','m'),
 108     CHAR2INT('p','a','r','a'),
 109     CHAR2INT('e','x','p','o'),
 110     CHAR2INT('s','i','n','c'),
 111     CHAR2INT('c','o','s','i'),
 112     CHAR2INT('s','h','e','p')
 113   };
 114
 115
 116   BOOLEAN eol;
 117   REAL theta,sinth,costh;
 118
 119   REAL kpr;
 120   REAL* g = NULL;
 121   REAL* g1= NULL;
 122
 123   REAL* q1 = NULL;
 124   REAL* q2 = NULL;
 125   REAL*  detecn = NULL;
 126
 127
 128   ////////////////////
 129
 130   INTEGER interp;
 131   int i;
 132   INTEGER missp;
 133   INTEGER nconr;
 134   REAL wconv;
 135   INTEGER method;
 136   REAL wconw;
 137   INTEGER nf;
 138   INTEGER nraym;
 139   REAL fmdray;
 140
 141   INTEGER mdraym;
 142   INTEGER mraym2;
 143   REAL halfl;
 144
 145   REAL rcube;
 146
 147   REAL alfa;
 148
 149   REAL pcovst;
 150
 151   REAL sigma;
 152
 153   REAL ctime;
 154
 155   REAL btime;
 156
 157   REAL beg;
 158
 159   REAL end;
 160   REAL ftime;
 161   REAL* ncos = NULL;
 162
 163   REAL limmid;
 164
 165   REAL* mdmod;
 166   REAL* modi = NULL;
 167   REAL* mdg;
 168   REAL* mdg1;
 169   REAL anginy;
 170
 171   REAL angiin;
 172
 173   REAL sbtrin;
 174
 175   REAL factor;
 176
 177   INTEGER lpj;
 178   REAL beta0;
 179
 180   int np;
 181   REAL dtc;
 182
 183   REAL datm;
 184   REAL dath;
 185   REAL datl;
 186
 187   REAL cbeg;
 188   int k;
 189   REAL* mr1;
 190   REAL* ml1;
 191   REAL* mr;
 192   REAL* ml;
 193   REAL pr1;
 194   REAL pr2;
 195   REAL pl1;
 196   REAL pl2;
 197   REAL* mq1;
 198   REAL* mq2;
 199   int i1;
 200   REAL cend;
 201   REAL beta1;
 202
 203   REAL sbeta;
 204
 205   REAL cbeta;
 206
 207   REAL w;
 208
 209   REAL b;
 210   REAL cbhl;
 211   REAL sbhl;
 212   REAL cbpxs;
 213
 214   REAL sbpxs;
 215
 216   REAL usini;
 217   REAL ucosi;
 218
 219   REAL bovrc;
 220   REAL wbkpjy;
 221
 222   REAL wincy;
 223
 224   REAL wincx;
 225
 226   REAL angsin;
 227   REAL angcos;
 228
 229   REAL anginc;
 230
 231   REAL subtr;
 232
 233   INTEGER ind;
 234
 235   int iy;
 236   REAL usin;
 237   REAL ucos;
 238   REAL wbkpj;
 239   REAL ang;
 240   int ix;
 241
 242   REAL bend;
 243
 244   BOOLEAN ret_xalg4;
 245
 246   ////////////////////
 247
 248   if(GeoPar.nrays <= 1) {
 249
 250     fprintf(output, "\n         ***** nrays too small - job aborted *****");
 251
 252     return TRUE;
 253   }
 254
 255   // ALL BUT THE FIRST ITERATION ARE FOR SMOOTHING ONLY
 256
 257   ret_xalg4 = FALSE;
 258   if(iter > 1) return false;
 259
 260   // CHECK DATA COLLECTING GEOMETRY
 261
 262   if(!GeoPar.div) {
 263     fprintf(output, "\n      *****geometry error*****");
 264     fprintf(output, "\n\n     this algorithm is written for fan beam only \n\n");
 265     return TRUE;
 266   }
 267
 268   // READ CHECK AND ECHO USER SUPPLIED INFORMATION
 269
 270   interp = InFile.getint(TRUE, &eol);
 271
 272   if((interp < -1) || (interp > 6)) {
 273
 274     fprintf(output, "\n          ***** illegal value for interp *****");
 275     ret_xalg4 = TRUE;
 276   }
 277
 278   i = GeoPar.nrays - GeoPar.snrays;
 279   if(((interp == -1) && (i < 2)) || ((interp > 2) && (i < (interp-2)))) {
 280
 281     fprintf(output, "\n          ***** size of modified projection data too small");
 282     fprintf(output, "\n               for interpolation method specified *****");
 283
 284     ret_xalg4 = TRUE;
 285   }
 286
 287   missp = InFile.getint(FALSE, &eol);
 288
 289   if((missp <= 0) || (missp >= GeoPar.prjnum)) {
 290     fprintf(output, "\n         ***** illegal value for missp *****");
 291
 292     ret_xalg4 = TRUE;
 293   }
 294
 295
 296   nconr = InFile.getint(FALSE, &eol);
 297   wconv = InFile.getnum(FALSE, &eol);
 298   if((wconv <= 0) || (wconv > 1)) {
 299     fprintf(output, "\n         ***** illegal value for wconv *****");
 300
 301         ret_xalg4 = TRUE;
 302   }
 303
 304
 305   method = InFile.getint(FALSE, &eol);
 306
 307   if(eol) {
 308     fprintf(output, "\n         ***** insufficient data item - end of line encountered *****");
 309
 310     return TRUE;
 311   }
 312
 313   if(ret_xalg4) return ret_xalg4;
 314
 315   if(interp > 0) {
 316     fprintf(output, "\n          %1i point lagrange interpolation", interp);
 317   }
 318
 319   if(interp < 0) {
 320     fprintf(output, "\n          cubic spline interpolation");
 321   }
 322
 323   if(interp == 0) {
 324     fprintf(output, "\n          band limiting (sinc) interpolation");
 325   }
 326
 327   fprintf(output, "\n          one out of every %5i projection(s) used for reconstruction", missp);
 328
 329
 330   if((nconr < 0) || (nconr >= GeoPar.nrays)) nconr = GeoPar.nrays - 1;
 331
 332   fprintf(output, "\n          %3i points used on each side of convolution", nconr);
 333
 334   wconw = (REAL) 0.5 * ((REAL) 1.0 - wconv);
 335   fprintf(output, "\n          pre-smoothing weight on center ray = %6.4f", wconv);
 336   fprintf(output, "\n          pre-smoothing weight on side rays  = %6.4f", wconw);
 337
 338
 339   if(method > 0) {
 340     fprintf(output, "\n          weight for back-projection approximated");
 341   }
 342
 343   if(method == 0) {
 344     fprintf(output, "\n          accurate back-projection");
 345   }
 346
 347   if(method < 0) {
 348     fprintf(output, "\n          both weight and ray number for back-projection approximated");
 349   }
 350
 351   nf = InFile.getwrd(TRUE, &eol, dcon_codes, 8);
 352   if(eol) {
 353     fprintf(output, "\n         ***** insufficient data item - end of line encountered *****");
 354
 355     return TRUE;
 356   }
 357
 358   // USE HERMAN, LAKSHMINARAYANAN AND NAPARSTEK'S ALGORITHM FOR
 359   // BAND-LIMITING FILTER
 360
 361   if(nf == band) {
 362     fprintf(output, "\n          band-limiting filter");
 363
 364     dconbl(recon, interp, missp, nconr, wconv, method);
 365     return ret_xalg4;
 366   }
 367
 368   // INITIALIZE CONSTANTS
 369
 370   nraym = GeoPar.nrays - 1;
 371   fmdray = (REAL) GeoPar.midray + 1;
 372   mdraym = GeoPar.midray;
 373   mraym2 = GeoPar.midray - 1;
 374   halfl = (REAL) GeoPar.midpix * GeoPar.pixsiz;
 375   rcube = GeoPar.radius * GeoPar.radius * GeoPar.radius;
 376   alfa = GeoPar.pinc / GeoPar.stod;
 377
 378   // INITIALIZE CONSTANTS FOR TANGENT DETECTOR CASE
 379   // DETECN IS THE BASE FOR THE ARRAY WHICH STORES THE
 380   // ARC DETECTOR LOCATIONS IN TANGENT DETECTOR UNIT
 381
 382   if(!(GeoPar.arc || (GeoPar.nrays <= 3))) {
 383     pcovst = alfa;
 384     alfa = (REAL) atan((REAL) mdraym * alfa) / (REAL) mdraym;
 385
 386
 387     detecn = new REAL[mraym2];
 388
 389     sigma = alfa;
 390     for(i = 0; i < mraym2; i++) {
 391       detecn[i] = (REAL)( (REAL) sin(sigma) / (REAL) cos(sigma)) / pcovst;
 392       sigma += alfa;
 393     }
 394   }
 395
 396   // INITIALIZE TIMERS
 397
 398   ctime = 0.;
 399   btime = 0.;
 400   second(&beg);
 401
 402   // ALLOCATE WORK AREAS FOR CONVOLUTING FUNCTION
 403   // Q1 IS THE INDEX FOR THE FIRST ELEMENT IN Q1(M)
 404   // Q2 IS THE INDEX FOR THE FIRST ELEMENT IN Q2(M)
 405
 406   q1 = new REAL[GeoPar.nrays+1];        //(JD 11/06)
 407   q2 = new REAL[GeoPar.nrays+1];        //(JD 11/06)
 408
 409   // CONVOLUTING FUNCTION
 410
 411   if(dconft(q1, q2, nconr+1, alfa, nf)) {
 412
 413
 414   if(detecn != NULL) delete[] detecn ;
 415   if(q1 != NULL) delete[] q1 ;
 416   if(q2 != NULL) delete[] q2 ;
 417   if(ncos != NULL) delete[] ncos ;
 418   if(g != NULL) delete[] g ;
 419   if(g1 != NULL) delete[] g1;
 420   if(modi != NULL) delete[] modi ;           //wei 3/2005
 421
 422     return TRUE;
 423   }
 424   second(&end);
 425   ftime = end - beg;
 426
 427   // WEIGHT ARRAY FOR CONVOLUTION
 428   // NCOS IS THE BASE FOR CONVOLUTION WEIGHT ARRAY
 429
 430   ncos = new REAL[mdraym];
 431
 432   sigma = alfa;
 433
 434   for(i = 0; i < mdraym; i++) {
 435     ncos[i] = (REAL) cos(sigma);
 436     sigma += alfa;
 437   }
 438
 439   // ALLOCATE AREAS FOR PROJECTION DATA
 440   // G IS THE BASE FOR PROJECTION DATA FOR ONE PROJECTION
 441   // G1 IS THE BASE FOR WEIGHTED PROJECTION DATA
 442   // MODI IS THE BASE FOR MODIFIED PROJECTION DATA
 443
 444   g = new REAL[GeoPar.nrays];
 445   g1 = new REAL[GeoPar.nrays];
 446
 447   modi = new REAL[GeoPar.nrays];
 448
 449   limmid = (REAL) MIN0(nconr, mdraym);
 450
 451   //mdmod = modi + GeoPar.midray + 1;
 452   //mdg = g + GeoPar.midray + 1;
 453   //mdg1 = g1 + GeoPar.midray + 1;
 454   mdmod = modi + GeoPar.midray;        //(JD 11/06)
 455   mdg = g + GeoPar.midray;        //(JD 11/06)
 456   mdg1 = g1 + GeoPar.midray;        //(JD 11/06)
 457
 458   anginy = SQR(GeoPar.pixsiz / GeoPar.radius);
 459   angiin = (REAL) GeoPar.midpix * anginy;
 460   sbtrin = (REAL) GeoPar.midpix * angiin;
 461   factor = (REAL) 0.25 * alfa * GeoPar.radius / Consts.pisq;
 462   //lpj = ((GeoPar.prjnum-1) / missp) * missp + 1;
 463   lpj = ((GeoPar.prjnum-1) / missp) * missp;    //new indexing (from zero) (JD)
 464
 465   Anglst.getang(lpj, &beta0, &sinth, &costh);
 466   //Anglst.getang(1, &theta, &sinth, &costh);
 467   Anglst.getang(0, &theta, &sinth, &costh);    //new indexing (from zero) (JD)
 468
 469   if((theta + Consts.twopi) <=  beta0) {
 470
 471     fprintf(output, "\n         ***** illegal arrangement of projection angles *****");
 472
 473
 474
 475   if(detecn != NULL) delete[] detecn ;
 476   if(q1 != NULL) delete[] q1 ;
 477   if(q2 != NULL) delete[] q2 ;
 478   if(ncos != NULL) delete[] ncos ;
 479   if(g != NULL) delete[] g ;
 480   if(g1 != NULL) delete[] g1;
 481   if(modi != NULL) delete[] modi ;           //wei 3/2005
 482
 483     return TRUE;
 484   }
 485
 486   // THIS IS WHERE THE RECONSTRUCTION STARTS
 487   // BUILD UP RECONSTRUCTED PICTURE ONE PROJECTION AT A TIME
 488
 489   for(np = 0; np < GeoPar.prjnum; np += missp) {
 490     ProjFile.ReadProj(np, g, GeoPar.nrays);
 491
 492     // CONVERT THE PROJECTION DATA TO ARC GEOMETRY BY LINEAR
 493     // INTERPOLATION IF NECESSARY
 494
 495     if(!(GeoPar.arc || GeoPar.nrays <= 3)) {
 496
 497       *mdg1 = *mdg;
 498
 499       //for(i = 0; i < mraym2; i++){
 500       for(i = 1; i <= mraym2; i++){        //(JD 11/06)
 501         dtc = detecn[i-1];
 502         *(mdg1+i) = qintp(fmdray + dtc, g, GeoPar.nrays, 2);
 503         *(mdg1-i) = qintp(fmdray - dtc, g, GeoPar.nrays, 2);
 504       }
 505
 506       for(i = 1; i < nraym; i++){
 507         g[i] = g1[i];
 508       }
 509
 510     }
 511
 512     // PRE-SMOOTHING OF PROJECTION DATA
 513
 514     if(wconv != 1.0) {
 515       datm = g[0];
 516       dath = g[1];
 517       g[0] = wconv * datm + wconw * dath;
 518
 519       //for(i = 1; i < nraym; i++) {
 520       for(i = 1; i < nraym; i++) {         //(JD 11/06)
 521         datl = datm;
 522         datm = dath;
 523         dath = g[i+1];
 524         g[i] = wconv * datm + wconw * (dath + datl);
 525       }
 526
 527       g[GeoPar.nrays-1] = wconv * dath + wconw * datm;
 528     }
 529
 530     // WEIGHTED PROJECTION DATA
 531
 532     *(mdg1) = *(mdg);
 533
 534     //for(i = 0; i < mdraym; i++) {
 535     for(i = 1; i <= mdraym; i++) {        //(JD 11/06)
 536       //*(mdg1+i) = *(mdg+i) * ncos[i];
 537       //*(mdg1-i) = *(mdg-i) * ncos[i];
 538       *(mdg1+i) = *(mdg+i) * ncos[i-1];      //(JD 11/06)
 539       *(mdg1-i) = *(mdg-i) * ncos[i-1];      //(JD 11/06)
 540     }
 541
 542     // CONVOLUTION
 543
 544     second(&cbeg);
 545     mdmod[0] = *q1 * (*mdg1) + (*q2) * (*mdg);
 546     if(nconr != 0) {
 547
 548       //for(i = 0; i < limmid; i++) {
 549       for(i = 1; i <= limmid; i++) {      //(JD 11/06)
 550         mdmod[0] += *(q1+i) * ( *(mdg1+i) + *(mdg1-i) )+ *(q2+i) * ( *(mdg+i) + *(mdg-i));
 551       }
 552     }
 553
 554     //for(k = 0; k < mdraym; k++) {
 555     for(k = 1; k <= mdraym; k++) {                 //(JD 11/06)
 556       mr1 = mdg1+k;
 557       ml1 = mdg1-k;
 558       mr = mdg+k;
 559       ml = mdg-k;
 560       pr1 = *q1 * (*mr1);
 561       pr2 = *q2 * (*mr);
 562       pl1 = *q1 * (*ml1);
 563       pl2 = *q2 * (*ml);
 564       if(nconr != 0) {
 565
 566         //for(i = 0; i < nconr; i++) {
 567         for(i = 1; i <= nconr; i++) {              //(JD 11/06)
 568           if((k+i) > mdraym) {
 569
 570             //for(i1 = i; i1 < nconr; i1++) {
 571             for(i1 = i; i1 <= nconr; i1++) {       //(JD 11/06)
 572
 573               //if(!((mr1-i1) <= g1)) {
 574               if(!((mr1-i1) <= (g1-1))) { //i.e. if (mr1-i1) <g1
 575                 mq1 = q1 + i1;
 576                 mq2 = q2 + i1;
 577                 pr1 += *mq1 * *(mr1-i1);
 578                 pr2 += *mq2 * *(mr -i1);
 579                 pl1 += *mq1 * *(ml1+i1);
 580                 pl2 += *mq2 * *(ml +i1);
 581                 //break;        //(JD 11/06)
 582               }
 583             }
 584             break;
 585           }
 586           mq1 = q1 + i;
 587           mq2 = q2 + i;
 588           pr1 += *mq1 * (*(mr1+i) + *(mr1-i));
 589           pr2 += *mq2 * (*(mr +i) + *(mr -i));
 590           pl1 += *mq1 * (*(ml1+i) + *(ml1-i));
 591           pl2 += *mq2 * (*(ml +i) + *(ml -i));
 592         }
 593
 594       }
 595
 596       //*(mdmod + k) = pr1 + pr2 * ncos[k];
 597       //*(mdmod - k) = pl1 + pl2 * ncos[k];
 598       *(mdmod + k) = pr1 + pr2 * ncos[k-1];   // (JD 11/07/03)
 599       *(mdmod - k) = pl1 + pl2 * ncos[k-1];   // (JD 11/07/03)
 600
 601     }
 602
 603 #ifdef FFCOMPARE
 604     fprintf(output,"\n xalg4: modi");
 605     for(i = 0; i < GeoPar.nrays; i++) {
 606       if((i % 3) == 0) fprintf(output,"\n");
 607       fprintf(output," %36.30f", (double) modi[i]);
 608     }
 609 #endif
 610
 611     second(&cend);
 612     ctime += cend-cbeg;
 613
 614     // COMPUTE THE TRAPEZOIDAL WIDTH FOR BETA INTEGRATION
 615
 616     beta1 = theta;
 617     sbeta = -costh;
 618     cbeta = sinth;
 619     if(np < lpj) Anglst.getang(np + missp, &theta, &sinth, &costh);
 620     if((np < lpj) && (theta <= beta1)) {
 621       fprintf(output, "\n         ***** illegal arrangement of projection angles *****");
 622
 623
 624   if(detecn != NULL) delete[] detecn ;
 625   if(q1 != NULL) delete[] q1 ;
 626   if(q2 != NULL) delete[] q2 ;
 627   if(ncos != NULL) delete[] ncos ;
 628   if(g != NULL) delete[] g ;
 629   if(g1 != NULL) delete[] g1;
 630   if(modi != NULL) delete[] modi ;           //wei 3/2005
 631
 632       return TRUE;
 633     }
 634
 635     //if(np >= lpj) Anglst.getang(1,&theta, &sinth, &costh);
 636     if(np >= lpj) Anglst.getang(0,&theta, &sinth, &costh);         // (JD 11/07/03)
 637     w = (REAL) 0.5 * (REAL) fmod(Consts.twopi + theta - beta0, Consts.twopi);
 638
 639     beta0 = beta1;
 640     b = w * factor;
 641
 642     // BACK PROJECTION
 643
 644     cbhl = halfl * cbeta;
 645     sbhl = halfl * sbeta;
 646     cbpxs = GeoPar.pixsiz * cbeta;
 647     sbpxs = GeoPar.pixsiz * sbeta;
 648     usini = sbhl - cbhl;
 649     ucosi = GeoPar.radius - sbhl - cbhl;
 650
 651     if(method != 0) {
 652       bovrc = b / rcube;
 653       wbkpjy = (GeoPar.radius + (REAL) 2.0 * (cbhl + sbhl)) * bovrc;
 654       wincy = (REAL) -2.0 * cbpxs * bovrc;
 655       wincx = (REAL) -2.0 * sbpxs * bovrc;
 656
 657       if(method <= 0) {
 658         angsin = sbhl;
 659         angcos = GeoPar.radius - cbhl;
 660         anginc = cbpxs / GeoPar.radius - angiin;
 661         subtr = cbhl / GeoPar.radius - sbtrin;
 662       }
 663     }
 664
 665     ind = 0;
 666     for(iy = 0; iy < GeoPar.nelem; iy++) {
 667       if(method >= 0) {
 668         usini -= sbpxs;
 669         usin = usini;
 670         ucosi += cbpxs;
 671         ucos = ucosi;
 672       }
 673       if(method != 0) {
 674         wbkpjy = wbkpjy + wincy;
 675         wbkpj = wbkpjy;
 676         if(method <= 0) {
 677           angsin -= sbpxs;
 678           angcos += cbpxs;
 679           anginc += anginy;
 680           subtr += angiin;
 681           ang = (REAL) atan2(angsin, angcos) - subtr;
 682         }
 683       }
 684       for(ix = 0; ix < GeoPar.nelem; ix++) {
 685
 686         if(method >= 0) {
 687
 688           usin += cbpxs;
 689           ucos += sbpxs;
 690           ang = (REAL) atan2(usin, ucos);
 691         }
 692
 693         if(method == 0) {
 694
 695           wbkpj = b / (usin * usin + ucos * ucos);
 696         }
 697         else {
 698           wbkpj += wincx;
 699           if(method < 0) ang += anginc;
 700         }
 701         kpr = fmdray - ang / alfa;
 702         recon[ind] += qintp(kpr, modi, GeoPar.nrays, interp) * wbkpj;
 703         ind++;
 704       }
 705
 706     }
 707
 708     second(&bend);
 709     btime += bend - cend;
 710   }
 711   fprintf(output, "\n          total time for obtaining convoluting function was %10.3f seconds", ftime); // changed precision to three digits - swr 1/21/06
 712   fprintf(output, "\n          total time for convolution was                    %10.3f seconds", ctime); // ditto
 713   fprintf(output, "\n          total time for back-projection was                %10.3f seconds", btime); // ditto
 714
 715   if(detecn != NULL) delete[] detecn ;
 716   if(q1 != NULL) delete[] q1 ;
 717   if(q2 != NULL) delete[] q2 ;
 718   if(ncos != NULL) delete[] ncos ;
 719   if(g != NULL) delete[] g ;
 720   if(g1 != NULL) delete[] g1;
 721   if(modi != NULL) delete[] modi ;           //wei 3/2005
 722
 723   return ret_xalg4;
 724 }
 725