src/snark/dcon_dconbl.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_dconbl.cpp $
   5 $LastChangedRevision: 85 $
   6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7 $Author: agulati $
   8 **********************************************************************
   9
  10  CONVOLUTION-TYPE METHOD FOR DIVERGENT DATA WITH BAND-LIMITING
  11  FILTER.
  12      BASED ON
  13         CONVOLUTION RECONSTRUCTION TECHNIQUES FOR DIVERGENT BEAMS
  14         G. T. HERMAN, A. V. LAKSHMINARAYANAN AND A. NAPARSTEK
  15         COMPUTERS IN BIOLOGY AND MEDICINE
  16         TO APPEAR
  17      WRITTEN BY G. T. HERMAN, MARCH 1975
  18      REVISED FOR THE CDC3500, MARCH 1976
  19      REVISED FOR SNARK77 BY T. CHANG, MAY 1977
  20  ***********************************************************
  21 */
  22
  23 #include <cstdio>
  24 #include <cmath>
  25
  26 #include "blkdta.h"
  27 #include "geom.h"
  28 #include "consts.h"
  29 #include "uiod.h"
  30 #include "anglst.h"
  31 #include "qintp.h"
  32 #include "second.h"
  33 #include "projfile.h"
  34
  35 #include "dcon.h"
  36
  37 #ifdef FFCOMPARE
  38
  39 REAL Sin(REAL x)
  40 {
  41   REAL sum = 0;
  42
  43   REAL y = x;
  44   REAL z = 1;
  45   REAL a = 1;
  46
  47   for(int i = 1; i < 10; i++) {
  48     if((i % 2) == 0) {
  49       sum -= z/a;
  50     }
  51     else {
  52       sum += z/a;
  53     }
  54
  55     a *= (2 * i + 1) * 2 * i; // !
  56     z *= x * x; // x ^ (2 *  i)
  57   }
  58
  59   return x * sum;
  60 }
  61
  62 #else
  63
  64 #define Sin(X) sin((X))
  65
  66 #endif
  67
  68 void dcon_class::dconbl(REAL* recon, INTEGER interp, INTEGER missp, INTEGER nconr, REAL wconv, INTEGER method)
  69 {
  70 ////////////////////////////////////////////////////////
  71
  72   INTEGER limup;
  73   REAL wconw;
  74   REAL* j = NULL;
  75   INTEGER mdrym1;
  76   REAL* mq = NULL;
  77   REAL alfa;
  78   REAL alfasq;
  79   REAL a;
  80   int i;
  81   int m;
  82   int n;
  83
  84   INTEGER nraym1;
  85   REAL concen;
  86   INTEGER limmid;
  87
  88   REAL halfl;
  89   REAL twodr;
  90   REAL fmidry;
  91   INTEGER lastpr;
  92   REAL beta0;
  93   REAL sinth;
  94   REAL costh;
  95   REAL theta;
  96   REAL stime;
  97   REAL ctime;
  98   REAL btime;
  99   REAL bend;
 100   int np;
 101   REAL dtam;
 102   REAL dtar;
 103   int nr;
 104   REAL dtal;
 105   REAL send;
 106   INTEGER mqm;
 107   int nm;
 108   INTEGER inp;
 109
 110   INTEGER inm;
 111   INTEGER inr;
 112   REAL pr;
 113   INTEGER inl;
 114   REAL pl;
 115   int mm;
 116   REAL cend;
 117   REAL cbeta;
 118   REAL sbeta;
 119   REAL beta1;
 120   REAL w;
 121   REAL rcosi;
 122   REAL rcosy;
 123   REAL rcosx;
 124   REAL rsini;
 125   REAL rsiny;
 126   REAL rsinx;
 127   REAL wtwodr;
 128   REAL wjrfi;
 129   REAL wjrfx;
 130   REAL wjrfy;
 131   REAL anginy;
 132   REAL angiin;
 133   REAL fkjrfx;
 134   REAL subtr;
 135   REAL sob;
 136   REAL cob;
 137   REAL wrsq;
 138   REAL prini;
 139   INTEGER ind;
 140   REAL rcos;
 141   REAL rsin;
 142   REAL wjrf;
 143   REAL fract;
 144   REAL fkjrf;
 145   REAL* mod;
 146   int k;
 147   REAL usq;
 148   REAL frynum;
 149   REAL p;
 150
 151 ////////////////////////////////////////////////////////
 152   REAL* g;
 153
 154   // ASSIGN SPACE FOR SMOOTHED PROJECTION DATA, ECHO SMOOTHING WEIGHTS
 155
 156   g = new REAL[GeoPar.nrays];
 157   limup = GeoPar.nrays - 1; // index of last ray
 158
 159   if(!((wconv > 0.0) && (wconv <= 1.0))) wconv = (REAL) 1.0;
 160   wconw = ((REAL) 1.0 - wconv) / (REAL) 2.0;
 161
 162   // ALLOCATE SPACE FOR CONVOLUTION INFORMATION
 163
 164   j = new REAL[GeoPar.nrays];
 165
 166   mdrym1 = GeoPar.midray - 1;
 167
 168   //mq = new REAL[mdrym1];
 169   mq = new REAL[mdrym1 + 1];  //(JD 1/20/04)
 170
 171
 172   // WEIGHT ARRAY J(I)
 173
 174   alfa = GeoPar.pinc / GeoPar.stod;
 175   alfasq = alfa * alfa;
 176   a = alfa * GeoPar.radius;
 177   j[GeoPar.midray] = a;
 178
 179   for(i = 1; i <= GeoPar.midray; i++)
 180   {
 181     if(GeoPar.tang) j[GeoPar.midray + i] = a / (REAL) sqrt((REAL) 1.0 + alfasq * SQR(((REAL)(i))));
 182     if(GeoPar.arc) j[GeoPar.midray + i]  = a * (REAL) cos( ((REAL)(i)) * alfa);
 183     j[GeoPar.midray - i] = j[GeoPar.midray + i];
 184   }
 185
 186   ////////////////////////////////////////////////////
 187 #ifdef FFCOMPARE
 188   fprintf(output,"\n***** test data 1 *****************************");
 189   for(i = 0; i < GeoPar.nrays; i++)
 190   {
 191     if((i % 3) == 0) fprintf(output,"\n");
 192     fprintf(output," %36.30f", j[i]);
 193   }
 194   fprintf(output,"\n");
 195 #endif
 196   ////////////////////////////////////////////////////
 197
 198
 199   // CONVOLUTING FUNCTION Q(M)
 200
 201   for(m = 0; m < GeoPar.midray; m++)
 202   {
 203     n = 2 * m + 1; // 1, 3, 5, ...
 204     if(GeoPar.tang) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * n * a));
 205 #ifndef FFCOMPARE
 206     if(GeoPar.arc) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * (REAL) sin((REAL) (n * alfa))));
 207
 208 #else
 209     REAL xxx = n * alfa;
 210     REAL yyy = Sin(xxx);
 211     REAL zzz = sin(xxx);
 212
 213     if(GeoPar.arc) mq[m] = (REAL) -1.0 / ((REAL) 2.0 * SQR(Consts.pi * yyy));
 214
 215 //    fprintf(output,"\n %36.30f %36.30f %36.30f", (double) mq[m], (double) xxx, (double) alfa);
 216     fprintf(output,"\n %36.30f", (double) mq[m]);
 217 #endif
 218   }
 219
 220   ////////////////////////////////////////////////////
 221 #ifdef FFCOMPARE
 222   fprintf(output,"\n***** test data 2*****************************");
 223   for(m = 0; m < GeoPar.midray; m++)
 224   {
 225     if((m % 3) == 0) fprintf(output,"\n");
 226     fprintf(output," %36.30f", mq[m]);
 227   }
 228   fprintf(output,"\n");
 229 #endif
 230   ////////////////////////////////////////////////////
 231
 232
 233
 234   // CONSTANTS NEEDED FOR CONVOLUTION AND BACKPROJECTION
 235
 236   nraym1 = GeoPar.nrays - 1;
 237   if(GeoPar.tang) concen = (REAL) 1.0 / ((REAL) 8.0 * SQR(a));
 238   if(GeoPar.arc) concen = (REAL) 1.0 / ((REAL) 8.0 * SQR(alfa));
 239   limmid = MIN0(nconr, mdrym1);
 240   mod = new REAL[GeoPar.nrays];
 241
 242   ///midmod = mod + midray;
 243   halfl = ((REAL)(GeoPar.midpix)) * GeoPar.pixsiz;
 244   twodr = (REAL) 2.0 / GeoPar.radius;
 245   fmidry = (REAL) GeoPar.midray + 1;
 246
 247   // INITIALIZE BETA0 TO BE USED IN CALCULATING WEIGHTS FOR
 248   // BETA INTEGRATION
 249
 250   lastpr = ((GeoPar.prjnum - 1) / missp) * missp;
 251   Anglst.getang(lastpr, &beta0, &sinth, &costh);
 252   Anglst.getang(0, &theta, &sinth, &costh);
 253
 254
 255   // INITIALIZE TIMERS
 256
 257   stime = 0.0;
 258   ctime = 0.0;
 259   btime = 0.0;
 260   second(&bend);
 261
 262
 263   // THIS IS WHERE THE RECONSTRUCTION STARTS
 264   // BUILD UP RECONSTRUCTED PICTURE ONE PROJECTION AT A TIME
 265
 266   for(np = 0; np < GeoPar.prjnum; np += missp)
 267   {
 268
 269     ProjFile.ReadProj(np, g, GeoPar.nrays);
 270
 271     // debug ////////////////////////////////////////////////////
 272     // input verification
 273 #ifdef FFCOMPARE
 274     fprintf(output,"\n**********************************");
 275     for(nr = 0; nr < GeoPar.nrays; nr++)
 276         {
 277       if((nr % 3) == 0) fprintf(output,"\n");
 278       fprintf(output," %36.30f", g[nr]);
 279     }
 280 #endif
 281     //////////////////////////////////////////////////////
 282
 283
 284     dtam = g[0];
 285     dtar = g[1];
 286     g[0] = (wconv * dtam + wconw * dtar) * j[0];
 287
 288     for(nr = 1; nr < nraym1; nr++)
 289         {
 290       dtal = dtam;
 291       dtam = dtar;
 292       dtar = g[nr + 1];
 293       g[nr] = (wconv * dtam + wconw * (dtar + dtal)) * j[nr];
 294     }
 295
 296     dtal = dtam;
 297     dtam = dtar;
 298     g[nraym1] = (wconv * dtam + wconw * dtal) * j[nraym1];
 299
 300     second(&send);
 301     stime += send - bend;
 302
 303
 304 #ifdef FFCOMPARE
 305     fprintf(output,"\n***** test data 3 *****************************");
 306     for(nr = 0; nr < nraym1; nr++) {
 307       if((nr % 3) == 0) fprintf(output,"\n");
 308       fprintf(output," %36.30f", g[nr]);
 309     }
 310     fprintf(output,"\n");
 311 #endif
 312
 313
 314     // COMPUTE CONVOLUTED PROJECTION DATA
 315
 316     p = concen * g[GeoPar.midray];
 317
 318     mqm = 0;
 319     for(nm = 1; nm <= limmid; nm += 2) {
 320       p += (g[GeoPar.midray + nm] + g[GeoPar.midray - nm]) * mq[mqm];
 321       mqm++;
 322     }
 323
 324     mod[GeoPar.midray] = p;
 325
 326     // WE HAVE NOW TAKEN CARE OF THE CENTER RAY
 327
 328     for(m = 1; m <= GeoPar.midray; m++) {
 329       inr = GeoPar.midray + m;
 330       pr = concen * g[inr];
 331       inl = GeoPar.midray - m;
 332       pl = concen * g[inl];
 333
 334       if(!(nconr == 0)) {
 335
 336         mqm = 0;
 337
 338         for(nm = 1; nm <= nconr; nm += 2) {
 339           inp = inr + nm;
 340           if(inp > limup) {
 341             for(mm = nm; mm <= nconr; mm += 2) {
 342               inm = inr - mm;
 343               if(inm < 0) break;
 344               pr += g[inm] * mq[mqm];
 345               inp = inl + mm;
 346               pl += g[inp] * mq[mqm];
 347               mqm++;
 348             }
 349             break;
 350           }
 351
 352           inm = inr - nm;
 353
 354           pr += ( g[inp] + g[inm] ) * mq[mqm];
 355           inp = inl + nm;
 356           inm = inl - nm;
 357           pl += ( g[inp] + g[inm] ) * mq[mqm];
 358
 359
 360
 361           mqm++;
 362         }
 363       }
 364
 365       mod[GeoPar.midray + m] = pr;
 366       mod[GeoPar.midray - m] = pl;
 367
 368     }
 369
 370 #ifdef FFCOMPARE
 371
 372   fprintf(output,"\n***** test data mod *****************************");
 373   for(m = 0; m < GeoPar.nrays; m++) {
 374     if((m % 3) == 0) fprintf(output,"\n");
 375     fprintf(output," %36.30f", mod[m]);
 376   }
 377   fprintf(output,"\n");
 378
 379 #endif
 380   ////////////////////////////////////////////////////
 381
 382
 383     second(&cend);
 384     ctime += cend - send;
 385
 386     // BACK-PROJECT CONVOLUTED PROJECTION DATA
 387
 388     cbeta = sinth;
 389     sbeta = -costh;
 390
 391     // CALCULATE WEIGHTS FOR BETA INTEGRATION
 392
 393     beta1 = theta;
 394     if(np >= lastpr) Anglst.getang(0, &theta, &sinth, &costh);
 395     if(np < lastpr) Anglst.getang(np+missp, &theta, &sinth, &costh);
 396
 397     w = (REAL) fmod(Consts.twopi + theta - beta0, Consts.twopi) / (REAL) 2.0;
 398
 399     beta0 = beta1;
 400
 401     // CALCULATE BACK-PROJECTION INFORMATION
 402
 403     rcosi = halfl * (sbeta - cbeta);
 404     rcosy = - GeoPar.pixsiz * sbeta;
 405     rcosx = GeoPar.pixsiz * cbeta;
 406     rsini = -halfl * (sbeta + cbeta);
 407     rsiny = rcosx;
 408     rsinx = -rcosy;
 409     if(GeoPar.arc) w /= SQR(GeoPar.radius);
 410     wtwodr = w * twodr;
 411     wjrfi = w - wtwodr * rsini;
 412     wjrfx = -wtwodr * rsinx;
 413     wjrfy = -wtwodr * rsiny;
 414     anginy = (SQR(GeoPar.pixsiz / GeoPar.radius)) / alfa;
 415     angiin = ((REAL)(GeoPar.midpix)) * anginy;
 416     fkjrfx = rcosx / a - angiin;
 417     subtr = fkjrfx * GeoPar.midpix;
 418     sob = halfl * sbeta;
 419     cob = GeoPar.radius - halfl * cbeta;
 420     wrsq = w * SQR(GeoPar.radius);
 421     prini = GeoPar.radius + rsini;
 422
 423     // BACK-PROJECTION
 424
 425     ind = 0;
 426     for(i = 0; i < GeoPar.nelem; i++)
 427         {
 428                 if(method >= 0)
 429                 {
 430                         rcosi += rcosy;
 431                         rcos = rcosi;
 432                         prini += rsiny;
 433                         rsin = prini;
 434                 }
 435                 if(method != 0)
 436                 {
 437
 438                         wjrfi += wjrfy;
 439                         wjrf = wjrfi;
 440                         if(method <= 0) {
 441                                 sob += rcosy;
 442                                 cob += rsiny;
 443                                 fract = sob / cob;
 444                                 fkjrfx += anginy;
 445                                 subtr += angiin;
 446                                 if(GeoPar.arc) fract = (REAL) atan(fract);
 447                                 fkjrf = fract / alfa - subtr;
 448                         }
 449                 }
 450
 451                 for(k = 0; k < GeoPar.nelem; k++)
 452                 {
 453                         if(method >= 0)
 454                         {
 455                                 rcos += rcosx;
 456                                 rsin += rsinx;
 457
 458                                 if(method <= 0)
 459                                 {
 460                                         usq = SQR(rsin);
 461                                         if(GeoPar.arc) usq += SQR(rcos);
 462                                         wjrf = wrsq / usq;
 463                                 }
 464                                 fract = rcos / rsin;
 465                                 if(GeoPar.arc) fract = (REAL) atan(fract);
 466                                 fkjrf = fract / alfa;
 467                         }
 468                         if(method != 0)
 469                         {
 470                                 wjrf += wjrfx;
 471                                 if(method < 0) fkjrf += fkjrfx;
 472                         }
 473                         frynum = fmidry - fkjrf;
 474
 475
 476 #ifdef FFCOMPARE
 477                         fprintf(output,"\n frynum = %36.30f wjrf = %36.30f", frynum, wjrf);
 478 #endif
 479
 480
 481                         recon[ind] += wjrf * qintp(frynum, mod, GeoPar.nrays, interp);
 482
 483
 484                         ind++;
 485                 }
 486     }
 487
 488     second(&bend);
 489     btime += bend - cend;
 490   }
 491
 492   fprintf(output, "\n          total time for obtaining smoothed projection data was %9.3f", stime);   // changed precision to two digits - swr 1/21/06
 493   fprintf(output, "\n          total time for convoluting projection data was %16.3f", ctime);         // ditto
 494   fprintf(output, "\n          total time for back-projecting was %28.3f", btime);                     // ditto
 495
 496
 497   if(j != NULL) delete[] j;
 498   if(mq != NULL) delete[] mq ;           //wei 3/2005
 499
 500   delete[] mod;  // bug 92 - Lajos - 03/02/2005
 501   delete[] g;  // bug 92 - Lajos - 03/02/2005
 502   return;
 503 }