src/snark/ray.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/ray.cpp $
   5  $LastChangedRevision: 91 $
   6  $Date: 2014-07-02 17:34:27 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  RAY RETURNS A LIST OF POINTS THAT LIE IN THE N RAY-TH RAY TAKEN
  11  AT ANGLE THETA AND WITH RAY WIDTHS DX AND DY MEASURED IN
  12  THE X AND Y DIRECTIONS.
  13  NUM IS THE NUMBER OF POINTS IN THE RAY.
  14  LIST IS THE SET OF POINTS RETURNED.  FOR A POINT WITH COORDINATES
  15  (I,J), LIST = (J-1)*WIDTH + I
  16
  17  Algorithm rewritten and simplified as part of bug 82 - swr 2/20/05
  18  */
  19
  20 #include <cstdio>
  21 #include <cmath>
  22 #include <stdlib.h>
  23
  24 #include "blkdta.h"
  25 #include "geom.h"
  26 #include "uiod.h"
  27
  28 #include "posit.h"
  29
  30 #include "ray.h"
  31
  32 void ray(INTEGER nproj, INTEGER nray, INTEGER* list, REAL* weight,
  33                 INTEGER* numb, REAL* snorm)
  34 {
  35         int i, j;
  36
  37         //INTEGER lo,hi;    // bug82 - swr 2/20/05
  38
  39         static INTEGER lastp = -1; // made static to match the intent of the Fortran code. Lajos, Feb 10, 2005
  40
  41         REAL ax;
  42         REAL ay;
  43         REAL costh;
  44         REAL sinth;
  45         static REAL dee;
  46         static REAL pmid; // made static to match the intent of the Fortran code. Lajos, Feb 10, 2005
  47         static REAL sin_1;    // bug82 start - swr 2/20/05
  48         static REAL cos_1;
  49         static INTEGER inci; // made static to match the intent of the Fortran code. Lajos, Feb 10, 2005
  50         static INTEGER incj; // made static to match the intent of the Fortran code. Lajos, Feb 10, 2005
  51         //INTEGER iscrpt;    // bug82 end - swr 2/20/05
  52         REAL pix2;
  53         static INTEGER midpix;    // bug82 - swr 2/20/05
  54         static INTEGER nrlo, nrhi;
  55
  56         // bug 168 - swr - 9/24/05
  57         if (GeoPar.div)
  58         {
  59                 fprintf(output,
  60                                 "\n         ***** ray is not usable with divergent geometry\n");
  61                 exit(1);
  62         }
  63
  64         *numb = 0;
  65         posit(nproj, nray, &ax, &ay, &costh, &sinth);
  66         //nr = (REAL) (nray - GeoPar.midray);    // bug82 - swr 2/20/05
  67
  68         if (nproj != lastp)
  69         {
  70
  71                 // INITIALIZATION FOR A NEW PROJECTION ANGLE
  72
  73                 lastp = nproj;
  74                 if (GeoPar.uni)
  75                         dee = GeoPar.pinc;
  76                 if (GeoPar.vri)
  77                         dee = GeoPar.pinc * ((REAL) MAX0(fabs(costh), fabs(sinth)));
  78                 pmid = (GeoPar.nelem + 1) / ((REAL) 2.0);
  79
  80                 if (fabs(sinth) < fabs(costh))
  81                 {
  82                         sin_1 = sinth;
  83                         cos_1 = costh;
  84                         inci = 1;
  85                         incj = GeoPar.nelem;
  86                 }
  87                 // FOR ANGLES THAT ARE CLOSER TO THE Y-AXIS THAN TO THE X-AXIS,
  88                 // THE FOLLOWING EFFECTIVELY PERFORMS A ROTATION ABOUT THE LINE X=-Y.
  89                 else
  90                 {
  91                         sin_1 = costh;
  92                         cos_1 = sinth;
  93                         inci = GeoPar.nelem;
  94                         incj = 1;
  95                 }
  96                 midpix = GeoPar.midpix - 1;    // bug82 - swr 2/20/05
  97                 // find limits of the rays that can intersect the reconstruction region
  98                 REAL halfw = midpix * GeoPar.pixsiz;
  99                 INTEGER nr = (INTEGER) ceil(halfw * (fabs(costh) + fabs(sinth)) / dee);
 100                 nrlo = GeoPar.midray - nr;
 101                 nrhi = GeoPar.midray + nr;
 102         }
 103
 104         if (nray >= nrlo && nray <= nrhi)
 105         {
 106
 107                 // DETERMINE THE POINTS THAT LIE IN THE RAY
 108                 // The equation of a line making an angle atan2(sin_1, cos_1)
 109                 // with the x axis and a perpendicular distance r from the origin
 110                 // is given by y*cos_1 - x * sin_1 = r.  swr - 2/20/05
 111
 112                 // find the distance the two edges of the strip are from the origin
 113                 REAL r1 = ((REAL) (nray - GeoPar.midray) - 0.5) * dee;
 114                 REAL r2 = ((REAL) (nray - GeoPar.midray) + 0.5) * dee;
 115
 116                 for (i = -midpix; i <= midpix; ++i)
 117                 {
 118                         // for each column of pixels, find the two edges of the strip
 119                         REAL x = i * GeoPar.pixsiz;
 120                         REAL y1 = (r1 + x * sin_1) / cos_1;
 121                         REAL y2 = (r2 + x * sin_1) / cos_1;
 122                         // swap if necessary
 123                         if (y1 > y2)
 124                         {
 125                                 REAL t = y1;
 126                                 y1 = y2;
 127                                 y2 = t;
 128                         }
 129                         // convert to pixel coordinates
 130                         INTEGER jlo = MAX0((INTEGER)ceil(-y2 / GeoPar.pixsiz), -midpix);
 131                         INTEGER jhi = MIN0((INTEGER)ceil(-y1 / GeoPar.pixsiz), midpix + 1);
 132                         // add to the list
 133                         for (j = jlo; j < jhi; ++j)
 134                         {
 135                                 list[*numb] = (j + midpix) * incj + (i + midpix) * inci;
 136                                 if (list[*numb] < 0 || list[*numb] >= GeoPar.area)
 137                                         printf("\nray: error pixel out of range = %d", list[*numb]);
 138                                 ++*numb;
 139                         }
 140                 }
 141                 // end of new algorithm - swr 2/20/05
 142         }
 143
 144         pix2 = GeoPar.pixsiz * GeoPar.pixsiz;
 145
 146         for (i = 0; i < *numb; i++)
 147         {
 148                 weight[i] = pix2;
 149         }
 150
 151         *snorm = pix2 * pix2 * (*numb);
 152         if (trace > 6)
 153         {
 154                 fprintf(output,
 155                                 "\n          ray     np = %5i  nr = %5i  numb = %5i  snorm = %10.4f",
 156                                 nproj, nray, *numb, *snorm);
 157         }
 158
 159         if ((*numb > 0) && (trace > 9))
 160         {
 161                 for (i = 0; i < *numb; i++)
 162                 {
 163                         if (i % 10 == 0)
 164                                 fprintf(output, "\n                  ");
 165                         fprintf(output, "%6i", list[i]);
 166                 }
 167         }
 168
 169         return;
 170 }