src/snark/back.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/back.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  COMBINED BACK-PROJECTION RECONSTRUCTION ALGORITHM.
  11
  12  THIS SUBROUTINE COMPUTES A BACK-PROJECTION RECONSTRUCTION
  13  USING EITHER THE EQUATIONS OF THE LINES AS GIVEN BY RAY OR WRAY
  14  OR BY INTERPOLATING BETWEEN THE GIVEN PROJECTION DATA TO
  15  APPROXIMATE THE RADON TRANSFORM AT A PARTICULAR POINT.  THE
  16  FORMER IS SPECIFIED BY THE 'DISCRETE' PARAMETER, THE LATTER BY
  17  'CONTINUOUS'.
  18
  19  IN THE DISCRETE CASE, THE USER HAS THE ADDITIONAL OPTIONS OF
  20  SPECIFYING IF THE RAYSUM IS TO BE DIVIDED BY 'SNORM' BEFORE
  21  BEING BACK-PROJECTED ('DISTRIBUTED') AND IF PIXELS WHICH ARE
  22  INTERCEPTED BY A ZERO OR NEGATIVE RAY SHOULD BE REMOVED FROM
  23  THE COMPUTATIONS ('ZEROIZING').
  24
  25  IN THE CONTINUOUS CASE, THE DESIRED INTERPOLATION METHOD MUST
  26  BE SPECIFIED.  WARNING- THE CONTINUOUS CASE WEIGHTS PROJECTION
  27  "NP" BY "(THETA(NP+1)-THETA(NP-1))/2.0".  IN ORDER TO ACCOMPLISH
  28  THIS IT IS ASSUMED THAT THE PROJECTIONS ARE ARRANGED IN INCREASING
  29  ORDER ON THETA AND THAT THETA(PRJNUM)-THETA(1) IS LESS THAN
  30  PI IN THE PARALLEL CASE OR TWOPI IN THE DIVERGENT CASE.
  31
  32  THE RESULT OF THE BACK-PROJECTION WILL BE NORMALIZIED TO
  33  THE CORRECT AVERAGE DENSITY IF EITHER 'ADDITIVE' OR 'MULTIPLICATIV
  34  HAS BEEN SPECIFIED.  THE FORMER PERFORMS THE NORMALIZATION BY
  35  THE ADDITION OF A CONSTANT TO THE ENTIRE PICTURE, THE LATTER BY
  36  MULTIPLICATION.
  37  */
  38
  39 #include <cstdio>
  40 #include <cmath>
  41
  42 #include "blkdta.h"
  43 #include "geom.h"
  44 #include "spctrm.h"
  45 #include "fourie.h"
  46 #include "consts.h"
  47 #include "uiod.h"
  48
  49 #include "bckprj.h"
  50 #include "anglst.h"
  51 #include "ray.h"
  52 #include "wray.h"
  53
  54 #include "projfile.h"
  55 #include "infile.h"
  56
  57 #include "back.h"
  58
  59 BOOLEAN back_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  60 {
  61         REAL neginf, zflag;
  62         INTEGER word;
  63         REAL* g = NULL;       //wei 3/2005
  64         BOOLEAN disc, zeros, dist, add, mult;
  65         BOOLEAN eol;
  66         static const INTEGER hcont = CHAR2INT('c', 'o', 'n', 't');
  67         static const INTEGER hdisc = CHAR2INT('d', 'i', 's', 'c');
  68         static const INTEGER hdist = CHAR2INT('d', 'i', 's', 't');
  69         static const INTEGER hzero = CHAR2INT('z', 'e', 'r', 'o');
  70         static const INTEGER haddi = CHAR2INT('a', 'd', 'd', 'i');
  71         static const INTEGER hmult = CHAR2INT('m', 'u', 'l', 't');
  72
  73         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
  74         static const INTEGER back_codes[6] =
  75         {
  76                         CHAR2INT('c', 'o', 'n', 't'),
  77                         CHAR2INT('d', 'i', 's', 'c'),
  78                         CHAR2INT('z', 'e', 'r', 'o'),
  79                         CHAR2INT('d', 'i', 's', 't'),
  80                         CHAR2INT('a', 'd', 'd', 'i'),
  81                         CHAR2INT('m', 'u', 'l', 't') };
  82
  83         INTEGER nsize;
  84         REAL thmod;
  85         REAL thfac;
  86         REAL theta;
  87         REAL sinth;
  88         REAL costh;
  89         REAL theta0;
  90         INTEGER np;
  91         REAL theta1;
  92         INTEGER np2;
  93         REAL theta2;
  94         REAL sinth2;
  95         REAL costh2;
  96         REAL w;
  97         INTEGER nr;
  98         INTEGER numb;
  99         REAL snorm;
 100         INTEGER nb;
 101         INTEGER k;
 102         REAL amt;
 103         REAL sum;
 104         INTEGER knt;
 105         REAL corr;
 106
 107         if (iter > 1)
 108                 return FALSE;
 109         neginf = -Consts.infin;
 110         zflag = neginf + neginf;
 111         dist = FALSE;
 112         zeros = FALSE;
 113         add = FALSE;
 114         mult = FALSE;
 115
 116         // PROCESS CONTROL CARD
 117
 118         word = InFile.getwrd(TRUE, &eol, back_codes, 2);
 119         if (word != hcont)
 120         {
 121                 if (word != hdisc)
 122                 {
 123                         fprintf(output, "\n **** neither continuous or discrete specified");
 124                         return TRUE;
 125                 }
 126
 127                 disc = TRUE;
 128
 129                 word = InFile.getwrd(FALSE, &eol, &(back_codes[2]), 4);
 130
 131                 if (word == hzero)
 132                 {
 133                         zeros = TRUE;
 134                         word = InFile.getwrd(FALSE, &eol, &(back_codes[3]), 3);
 135                 }
 136
 137                 if (word == hdist)
 138                 {
 139                         dist = TRUE;
 140                         word = InFile.getwrd(FALSE, &eol, &(back_codes[4]), 2);
 141                 }
 142
 143         }
 144         else
 145         {
 146                 disc = FALSE;
 147                 Fourie.interp = InFile.getint(FALSE, &eol);
 148                 if ((eol) || (Fourie.interp < -1) || (Fourie.interp > 6))
 149                 {
 150                         // ERROR CONDITIONS
 151
 152                         fprintf(output, "\n **** invalid interpolation method");
 153
 154                         return TRUE;
 155                 }
 156
 157                 word = InFile.getwrd(FALSE, &eol, &(back_codes[4]), 2);
 158         }
 159
 160         if (!eol)
 161         {
 162                 if (word == haddi)
 163                         add = TRUE;
 164                 if (word == hmult)
 165                         mult = TRUE;
 166         }
 167
 168         if (disc)
 169         {
 170                 fprintf(output, "\n          discrete back-projection");
 171         }
 172         else
 173         {
 174
 175                 fprintf(output, "\n          continuous back-projection");
 176
 177                 if (Fourie.interp == -1)
 178                 {
 179                         fprintf(output, "\n          modified cubic spline interpolation");
 180                 }
 181
 182                 if (Fourie.interp == 0)
 183                 {
 184                         fprintf(output, "\n          band limiting sinc interpolation");
 185                 }
 186
 187                 if (Fourie.interp > 0)
 188                 {
 189                         fprintf(output, "\n          %2i order lagrange interpolation",
 190                                         Fourie.interp);
 191                 }
 192         }
 193
 194         if (dist)
 195         {
 196                 fprintf(output, "\n          distributed back-projection");
 197         }
 198
 199         if (zeros)
 200         {
 201                 fprintf(output, "\n          zeroizing is in effect");
 202         }
 203
 204         if (add)
 205         {
 206                 fprintf(output, "\n          additive normalization");
 207         }
 208
 209         if (mult)
 210         {
 211                 fprintf(output, "\n          multiplicative normalization");
 212         }
 213
 214         if (!disc)
 215         {
 216
 217                 // CONTINUOUS CASE
 218
 219                 nsize = GeoPar.snrays + 2 * Fourie.interp;
 220                 if (Fourie.interp == -1)
 221                         nsize = GeoPar.snrays + 6;
 222
 223                 g = new REAL[nsize];
 224
 225                 thmod = Consts.pi;
 226                 thfac = 2.0;
 227                 if (!GeoPar.par)
 228                 {
 229                         thmod = Consts.twopi;
 230                         thfac = 4.0;
 231                 }
 232                 Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
 233                 theta0 = theta - thmod;
 234
 235                 for (np = 0; np < GeoPar.prjnum; np++)
 236                 {
 237
 238                         ProjFile.ReadProj(np, g, nsize);
 239
 240                         Anglst.getang(np, &theta1, &sinth, &costh);
 241
 242                         np2 = np + 1;
 243
 244                         if (np == (GeoPar.prjnum - 1))
 245                                 np2 = 0;
 246
 247                         Anglst.getang(np2, &theta2, &sinth2, &costh2);
 248
 249                         if (np == (GeoPar.prjnum - 1))
 250                                 theta2 += thmod;
 251
 252                         w = (theta2 - theta0) / thfac;
 253
 254                         if (w <= 0.0)
 255                         {
 256                                 fprintf(output,
 257                                                 "\n **** projection angles not in increasing order");
 258
 259                                 if (g != NULL)
 260                                         delete[] g;   //wei 3/2005
 261                                 return TRUE;
 262                         }
 263
 264                         Fourie.rinc = GeoPar.pinc;
 265
 266                         if (GeoPar.div)
 267                                 Fourie.rinc *= GeoPar.radius / GeoPar.stod;
 268                         if (GeoPar.vri)
 269                                 Fourie.rinc *= (REAL) MAX0(fabs(sinth), fabs(costh));
 270                         if (GeoPar.strip)
 271                                 w /= Fourie.rinc;
 272
 273                         theta0 = theta1;
 274
 275                         for (nr = 0; nr < nsize; nr++)
 276                         {
 277                                 g[nr] *= w;
 278                         }
 279
 280                         bckprj(recon, GeoPar.nelem, g, nsize, sinth, costh,
 281                                         GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
 282
 283                 }
 284
 285         }
 286         else
 287         {
 288                 // DISCRETE CASE
 289
 290                 g = new REAL[GeoPar.nrays];
 291
 292                 if (zeros)
 293                 {
 294
 295                         // FLAG PIXELS INTERCEPTED WITH ZERO OR NEGATIVE RAYS
 296
 297                         for (np = 0; np < GeoPar.prjnum; np++)
 298                         {
 299
 300                                 ProjFile.ReadProj(np, g, GeoPar.nrays);
 301
 302                                 for (nr = GeoPar.fsnray; nr <= GeoPar.lsnray; nr++)
 303                                 { // changed < to <=. lajos Nov 17, 2004
 304                                         if (g[nr] <= Consts.zero)
 305                                         {
 306
 307                                                 if (GeoPar.strip)
 308                                                         ray(np, nr, list, weight, &numb, &snorm);
 309                                                 if (GeoPar.line)
 310                                                         wray(np, nr, list, weight, &numb, &snorm);
 311
 312                                                 if (numb != 0)
 313                                                 {
 314
 315                                                         for (nb = 0; nb < numb; nb++)
 316                                                         {
 317                                                                 k = list[nb];
 318                                                                 recon[k] = zflag;
 319                                                         }
 320                                                 }
 321                                         }
 322                                 }
 323                         }
 324
 325                         // PERFORM DISCRETE BACK-PROJECTION
 326                 }
 327
 328                 for (np = 0; np < GeoPar.prjnum; np++)
 329                 {
 330
 331                         ProjFile.ReadProj(np, g, GeoPar.nrays);
 332                         for (nr = GeoPar.fsnray; nr <= GeoPar.lsnray; nr++)
 333                         { // changed < to <=. lajos Nov 17, 2004
 334
 335                                 if (GeoPar.strip)
 336                                         ray(np, nr, list, weight, &numb, &snorm);
 337                                 if (GeoPar.line)
 338                                         wray(np, nr, list, weight, &numb, &snorm);
 339
 340                                 if (numb != 0)
 341                                 {
 342                                         amt = g[nr];
 343                                         if (dist)
 344                                         {
 345                                                 if (zeros)
 346                                                 {
 347
 348                                                         // CORRECT SNORM BY REMOVING FLAGGED PIXELS
 349                                                         snorm = 0.0;
 350                                                         for (nb = 0; nb < numb; nb++)
 351                                                         {
 352                                                                 k = list[nb];
 353                                                                 if (recon[k] > neginf)
 354                                                                         snorm += weight[nb] * weight[nb];
 355                                                         }
 356                                                 }
 357                                                 amt /= snorm;
 358                                         }
 359                                         for (nb = 0; nb < numb; nb++)
 360                                         {
 361                                                 k = list[nb];
 362                                                 recon[k] += amt * weight[nb];
 363                                         }
 364                                 }
 365                         }
 366
 367                 }
 368         }
 369
 370         // PERFORM NORMALIZATION
 371
 372         if (!(add || mult))
 373         {
 374                 if (!zeros)
 375                 {
 376                         if (g != NULL)
 377                                 delete[] g;   //wei 3/2005
 378                         return FALSE;
 379                 }
 380                 // NO NORMALIZATION REQUESTED - JUST CLEAR FLAGGED PIXELS
 381
 382                 for (k = 0; k < GeoPar.area; k++)
 383                 {
 384                         if (recon[k] <= neginf)
 385                                 recon[k] = 0.0;
 386                 }
 387
 388                 if (g != NULL)
 389                         delete[] g;   //wei 3/2005
 390                 return FALSE;
 391
 392                 // COMPUTE TOTAL DENSITY AND NUMBER OF UNFLAGGED PIXELS
 393         }
 394
 395         sum = 0.0;
 396         knt = 0;
 397
 398         for (k = 0; k < GeoPar.area; k++)
 399         {
 400
 401                 if (recon[k] > neginf)
 402                 {
 403                         sum += recon[k];
 404                         knt++;
 405                 }
 406
 407         }
 408
 409         if (!mult)
 410         {
 411
 412                 // ADDITIVE CORRECTION TO AVERAGE DENSITY
 413
 414                 corr = GeoPar.aveden - sum / ((REAL) (knt));
 415
 416                 for (k = 0; k < GeoPar.area; k++)
 417                 {
 418                         recon[k] += corr;
 419                         if (recon[k] <= neginf)
 420                                 recon[k] = 0.0;
 421                 }
 422                 if (g != NULL)
 423                         delete[] g;   //wei 3/2005
 424                 return FALSE;
 425
 426                 // MULTIPLICATIVE CORRECTION TO AVERAGE DENSITY
 427         }
 428
 429         corr = ((REAL) (knt)) * GeoPar.aveden / sum;
 430
 431         for (k = 0; k < GeoPar.area; k++)
 432         {
 433                 if (recon[k] <= neginf)
 434                         recon[k] = 0.0;
 435                 recon[k] *= corr;
 436         }
 437         if (g != NULL)
 438                 delete[] g;   //wei 3/2005
 439         return FALSE;
 440 }
 441