src/snark/wray.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/wray.cpp $
   5  $LastChangedRevision: 121 $
   6  $Date: 2014-07-09 17:52:58 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  DEFINITION OF VARIABLES
  11  NUMB      THE NUMBER OF PIXELS A RAY PASSES THROUGH
  12  SNORM     SUM OF SQUARES OF WEIGHTS
  13  WEIGHT    THE LENGTH OF A RAY IN A PIXEL (CORRESPONDING TO
  14  ENTRY IN THE ARRAY  LIST)
  15  THETA     THE ANGLE IN RADIANS OF THE SOURCE OF THE RAY
  16  WITH RESPECT TO THE PICTURE AREA
  17  SORSTOC   THE DISTANCE FROM THE SOURCE TO THE CENTER OF THE BOD
  18  THE SELECTED PIXELS
  19  D         THE DISTANCE (ON THE TARGET )  FROM THE PERPENDICULAR
  20  THROUGH THE SOURCE TO THE PLACE AT WHICH THE
  21  DENSITY VALUE WAS MEASURED
  22  LIST      THE LIST OF INDECIES INTO THE ARRAY PICT THAT CORRESP
  23  TO THE SILECTED PIXELS
  24
  25  05/21/2008      jk      added 4 default parameters
  26  this allows to pass sorx, sory, cf and sf and avoid
  27  a call to posit() - used when taking projections
  28  of phantoms with variability when the call to posit()
  29  is made right before the call wray()
  30  */
  31
  32 #include <cstdio>
  33 #include <cmath>
  34
  35 #include "blkdta.h"
  36 #include "geom.h"
  37 #include "uiod.h"
  38
  39 #include "posit.h"
  40
  41 #include "wray.h"
  42
  43 void wray(INTEGER np, INTEGER nr, INTEGER* list, REAL* weight, INTEGER* numb, REAL* snorm, REAL * _sorx, REAL * _sory, REAL * _cf, REAL * _sf)
  44 //_sorx, _sory, _cf and _sf are default parametes with values set to REAL * _sorx=0, REAL * _sory=0, REAL * _cf=0, REAL * _sf=0
  45 {
  46         INTEGER flip;
  47
  48         REAL scale;
  49         REAL elem;
  50         REAL sorx;
  51         REAL sory;
  52         REAL cf;
  53         REAL sf;
  54         REAL hx;
  55         REAL hy;
  56         REAL tanphi;
  57         REAL rtnphi;
  58         REAL hxhx;
  59         REAL hyhy;
  60         REAL c;
  61         INTEGER ix;
  62         INTEGER iy;
  63         REAL fp;
  64         INTEGER incx;
  65         REAL w;
  66         REAL yc;
  67         REAL xc;
  68         INTEGER index;
  69
  70         int next;
  71
  72         scale = GeoPar.pixsiz;
  73         elem = (REAL) GeoPar.nelem;
  74         *numb = 0;
  75         *snorm = 0;
  76
  77         //jk  5/21/2008
  78         //if the values for _sorx and _sory are provided, do not call posit()
  79         if (_sorx != 0 && _sory != 0 && _cf != 0 && _sf != 0)
  80         {
  81                 sorx = *_sorx;
  82                 sory = *_sory;
  83                 cf = *_cf;
  84                 sf = *_sf;
  85         }
  86         else
  87                 posit(np, nr, &sorx, &sory, &cf, &sf);
  88
  89         sorx /= GeoPar.pixsiz;
  90         sory /= GeoPar.pixsiz;
  91
  92         if (sf < 0.0)
  93         {
  94                 sf = -sf;
  95                 cf = -cf;
  96         }
  97
  98         flip = 0;
  99         if (sf * cf < 0.0)
 100         {
 101                 flip = 1;
 102                 cf = -cf;
 103         }
 104
 105         /*
 106          IF SIN(PHI) IS TOO CLOSE TO ZERO SET TO SOME SMALL VALUE
 107          WITH PROPER SIGN
 108          LIKEWISE WITH COS(PHI)  THIS AVOIDS TROUBLES WITH DIVISION
 109          AND INFINITE SLOPES
 110          */
 111
 112         if (fabs(sf) <= 0.0000001)
 113         {
 114                 if (sf < 0)
 115                 {
 116                         sf = (REAL) -0.0000001;
 117                 }
 118                 else
 119                 {
 120                         sf = (REAL) 0.0000001;
 121                 }
 122         }
 123
 124         hx = ((REAL) 1.0) / sf;
 125
 126         if (((REAL) fabs(cf)) <= (REAL) 0.0000001)
 127         {
 128                 if (cf < 0)
 129                 {
 130                         cf = (REAL) -0.0000001;
 131                 }
 132                 else
 133                 {
 134                         cf = (REAL) 0.0000001;
 135                 }
 136         }
 137
 138
 139         hy = ((REAL) 1.0) / cf;
 140         tanphi = sf / cf;
 141         rtnphi = cf / sf;
 142         hx *= scale;
 143         hxhx = hx * hx;
 144         hy *= scale;
 145         hyhy = hy * hy;
 146
 147         // CALCULATE THE INTERSECTION OF REFERANCE LINE AND THE RAY IN QUESTI
 148
 149         sory += elem / ((REAL) 2.0);
 150
 151         if (flip > 0)
 152         {
 153                 // REFLECT IN Y=0 IF NEEDED
 154                 sorx = -sorx;
 155         }
 156
 157         sorx += elem / ((REAL) 2.0);
 158         c = sory - sorx * tanphi;
 159
 160         // FLIP LINE 180 DEGREES
 161
 162         c = (((REAL) 1.0) - tanphi) * elem - c;
 163
 164         // EQUATION OF RAY IS NOW Y=X*TANPHI+C
 165         // C IS THE Y INTERCEPT
 166         // IF INTERCEPT IS ABOVE (N) THEN THE RAY HAS MISSED ALL PIXELS
 167
 168         if ((c - elem) < 0)
 169         {
 170                 if (c < 0)
 171                 {
 172
 173                         // IF Y INTERCEPT IS NEGATIVE THEN TAKE X INTERCEPT (WHICH MUST
 174                         // BE POSITIVE SINCE SLOPE IS ALWAYS POSITIVE
 175
 176                         //  RAY ENTERS FROM BOTTOM
 177
 178                         c = -c * rtnphi;
 179
 180                         // RECALCULATE C SO THAT EQUATION OF LINE IS NOW
 181                         // X=RTANPHI*Y+C
 182                         // IF X INTERCEPT IS GREATER TNAN N THE RAY HAS MISSED THE RECONSTRUC
 183                         // AREA
 184
 185                         if ((c - elem) >= 0)
 186                                 goto L360;
 187
 188                         // INITIAL GRID (IX,IY) COORDINATE IS NOW CALCULATED
 189
 190                         ix = (INTEGER) c;
 191                         iy = 0;
 192
 193                         // FP IS ACTUAL DISPLACEMENT OF INTERCEPT FROM THE NEXT GREATER
 194                         // GRID LINE
 195
 196                         fp = ((REAL) (ix + 1)) - c;
 197
 198                         // START AT LOW X VALUES AND INCREMENT ACROSS PICTURE UNLESS
 199                         // SLOPE IS NEGATIVE
 200
 201                         incx = 1;
 202
 203                         // FLIP INDICATOR IS EXACTLY WRONG FOR HERMANS MATRIX STORAGE
 204
 205                         if (flip <= 0)
 206                         {
 207
 208                                 // START AT HIGH X VALUES AND DECREMENT (FROM RIGHT TO LEFT)
 209
 210                                 ix = GeoPar.nelem - ix - 1;
 211                                 incx = -1;
 212                         }
 213
 214                         // SET INITIAL VALUES OF THE LENGTH OF THE RAY IN CURRENT CELL
 215
 216                         w = hy * fp;
 217                         yc = fp * tanphi;
 218                         xc = -fp;
 219
 220                         // INDEX IS THE LINEAR ARRAY SUBSCRIPT OF THE  PICTURE ARRAY
 221                         // INDEX INTO FORTRAN TYPE ARRAY
 222
 223                         index = iy * GeoPar.nelem + ix + 1;
 224
 225                         // INITIAL CONDITIONS SET FOR ENTRANCE TO LOOPS WE HAVE ENTERED
 226                         // THE PICTURE AREA FROM THE BOTTOM
 227
 228                         next = 3;
 229                 }
 230                 else
 231                 {
 232
 233                         // ENTER PICTURE FROM SIDE (LEFT)
 234                         // CALCULATE INITIAL GRID COORDINATES (IX,IY)
 235
 236                         ix = 0;
 237                         incx = 1;
 238
 239                         // START WITH LOW X VALUES AND GO UP UNLESS THE FLIP INDICATOR
 240                         // IS SET
 241                         // FLIP INDICATOR IS EXACTLY WRONG FOR HERMANS MATRIX STORAGE
 242
 243                         if (flip <= 0)
 244                         {
 245                                 // START WITH HIGH X VALUES AND DECREMENT
 246
 247                                 ix = GeoPar.nelem - 1;
 248                                 incx = -1;
 249                         }
 250
 251                         iy = (INTEGER) c;
 252
 253                         // FP IS DISTANCE FROM Y INTERCEPT TO NEXT GREATER GRID LINE
 254
 255                         fp = ((REAL) (iy + 1)) - c;
 256
 257                         // CALCULATE ARRAY SUBSCRIPT IN PICTURE AREA OF INITIAL CELL
 258                         // INDEX INTO FORTRAN TYPE ARRAY
 259
 260                         index = iy * GeoPar.nelem + ix + 1;
 261
 262                         // CALCULATE INITIAL DESCRIPTION OF RAY IN FIRST CELL FOR LOOP
 263
 264                         yc = -fp;
 265                         w = hx * fp;
 266                         xc = fp * rtnphi;
 267
 268                         // RAY ENTERS FROM LEFT SIDE
 269
 270                         next = 1;
 271                 }
 272
 273                 for (;;)
 274                 {
 275                         switch (next)
 276                         {
 277                         case 1:
 278                                 //  MAKE ROOM FOR NEW PICTURE ELEMENT
 279                                 (*numb)++;
 280                                 list[(*numb) - 1] = index - 1;
 281
 282                                 // STEP ONE IN X DIRECTION AND CORRESPONDING STEP IN Y DIRECTION
 283
 284                                 xc -= (REAL) 1.0;
 285                                 yc += tanphi;
 286
 287                                 // IF XC IS NEGATIVE THE Y GRID LINE HAS BEEN CROSSED THE RAY
 288                                 // HAS LEFT THROUGHTHE TOP
 289
 290                                 if (xc > (REAL) 0.0)
 291                                 {
 292                                         // ENTERED SIDE,LEFT THROUGH SIDE --ASSIGN MAXIMUM Y-TYPE WEIGHT
 293                                         weight[(*numb) - 1] = hy;
 294                                         (*snorm) += hyhy;
 295                                         w -= hy;
 296                                         next = 4;
 297                                 }
 298                                 else
 299                                 {
 300                                         weight[(*numb) - 1] = w;
 301                                         (*snorm) += w * w;
 302
 303                                         if ((w - (REAL) 0.0000001) < 0)
 304                                         {
 305                                                 w = hy;
 306                                                 (*numb)--;
 307                                         }
 308                                         else
 309                                         {
 310                                                 w = hy - w;
 311                                         }
 312
 313                                         next = 2;
 314                                 }
 315                                 break;
 316                                 //////////////////////////////////////////////////////////////
 317                         case 2:
 318                                 iy++;
 319
 320                                 // RAY ENTERS FROM BOTTOM
 321                                 // STEP ONE CELL IN Y DIRECTION  AND TEST FOR COMPLETION
 322
 323                                 if ((iy - GeoPar.nelem) >= 0)
 324                                         goto L360;
 325
 326                                 index += GeoPar.nelem;
 327
 328                                 next = 3;
 329                                 break;
 330
 331                                 //////////////////////////////////////////////////////////////
 332
 333                                 // NEW LIST ELEMENT MEANS ANOTHER PIXEL SELECTED
 334                         case 3:
 335                                 (*numb)++;
 336                                 list[(*numb) - 1] = index - 1;
 337
 338                                 // STEP ONE IN Y DIRECTION AND CORRESPONDING X INCREMENT
 339
 340                                 yc -= 1.0;
 341                                 xc += rtnphi;
 342
 343                                 // IF YC IS NEGATIVE WH HAVE CROSSED A X-GRID LINE AND THAT
 344                                 // WE LEAVE THROUGE THE SIDE
 345
 346                                 if (yc > (REAL) 0.0)
 347                                 {
 348                                         w -= hx;
 349
 350                                         // ENTERED BOTTOM LEFT TOP ASSIGN MAXIMUM X-TYPE WEIGHT
 351
 352                                         weight[(*numb) - 1] = hx;
 353                                         (*snorm) += hxhx;
 354                                         next = 2;
 355                                 }
 356                                 else
 357                                 {
 358                                         weight[(*numb) - 1] = w;
 359                                         (*snorm) += w * w;
 360
 361                                         // IF THE WEIGHT SO CALCULATED IS TOO SMALL DELETE THIS PIXEL
 362                                         // THIS CAUSED BY CROSSING GRID LINES AT A CORNER
 363
 364                                         if ((w - (REAL) 0.000001) <= 0.0)
 365                                         {
 366                                                 w = hx;
 367                                                 (*numb)--;
 368                                         }
 369                                         else
 370                                         {
 371                                                 w = hx - w;
 372                                         }
 373                                         next = 4;
 374                                 }
 375                                 break;
 376
 377                                 /////////////////////////////////////////////////////////////
 378                         case 4:
 379                                 // ENTER PIXEL FROM SIDE
 380                                 ix += incx;
 381
 382                                 // STEP ONE IN X DIRECTION
 383                                 // AND TEST FOR TERMINATION
 384
 385                                 if (ix < 0)
 386                                         goto L360;
 387
 388                                 if ((ix - GeoPar.nelem) >= 0)
 389                                         goto L360;
 390                                 index += incx;
 391
 392                                 next = 1;
 393                                 break;
 394                         }
 395                 }
 396         }
 397
 398         L360:
 399
 400         if (trace > 6)
 401         {
 402                 fprintf(output, "\n          wray    np = %5i  nr = %5i  numb = %5i  snorm = %9.4f", np, nr, *numb, *snorm);
 403         }
 404         if ((*numb > 0) && (trace > 9))
 405         {
 406                 for (int i = 0; i < *numb; i++)
 407                 {
 408                         fprintf(output, "\n                  %5i     %9.4f", list[i], weight[i]);
 409                 }
 410         }
 411 }