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) $
8 ***********************************************************
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
17 Algorithm rewritten and simplified as part of bug 82 - swr 2/20/05
32 void ray(INTEGER nproj, INTEGER nray, INTEGER* list, REAL* weight,
33 INTEGER* numb, REAL* snorm)
37 //INTEGER lo,hi; // bug82 - swr 2/20/05
39 static INTEGER lastp = -1; // made static to match the intent of the Fortran code. Lajos, Feb 10, 2005
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
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
53 static INTEGER midpix; // bug82 - swr 2/20/05
54 static INTEGER nrlo, nrhi;
56 // bug 168 - swr - 9/24/05
60 "\n ***** ray is not usable with divergent geometry\n");
65 posit(nproj, nray, &ax, &ay, &costh, &sinth);
66 //nr = (REAL) (nray - GeoPar.midray); // bug82 - swr 2/20/05
71 // INITIALIZATION FOR A NEW PROJECTION ANGLE
77 dee = GeoPar.pinc * ((REAL) MAX0(fabs(costh), fabs(sinth)));
78 pmid = (GeoPar.nelem + 1) / ((REAL) 2.0);
80 if (fabs(sinth) < fabs(costh))
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.
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;
104 if (nray >= nrlo && nray <= nrhi)
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
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;
116 for (i = -midpix; i <= midpix; ++i)
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;
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);
133 for (j = jlo; j < jhi; ++j)
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]);
141 // end of new algorithm - swr 2/20/05
144 pix2 = GeoPar.pixsiz * GeoPar.pixsiz;
146 for (i = 0; i < *numb; i++)
151 *snorm = pix2 * pix2 * (*numb);
155 "\n ray np = %5i nr = %5i numb = %5i snorm = %10.4f",
156 nproj, nray, *numb, *snorm);
159 if ((*numb > 0) && (trace > 9))
161 for (i = 0; i < *numb; i++)
164 fprintf(output, "\n ");
165 fprintf(output, "%6i", list[i]);