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) $
8 ***********************************************************
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
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
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()
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
72 scale = GeoPar.pixsiz;
73 elem = (REAL) GeoPar.nelem;
78 //if the values for _sorx and _sory are provided, do not call posit()
79 if (_sorx != 0 && _sory != 0 && _cf != 0 && _sf != 0)
87 posit(np, nr, &sorx, &sory, &cf, &sf);
89 sorx /= GeoPar.pixsiz;
90 sory /= GeoPar.pixsiz;
106 IF SIN(PHI) IS TOO CLOSE TO ZERO SET TO SOME SMALL VALUE
108 LIKEWISE WITH COS(PHI) THIS AVOIDS TROUBLES WITH DIVISION
112 if (fabs(sf) <= 0.0000001)
116 sf = (REAL) -0.0000001;
120 sf = (REAL) 0.0000001;
124 hx = ((REAL) 1.0) / sf;
126 if (((REAL) fabs(cf)) <= (REAL) 0.0000001)
130 cf = (REAL) -0.0000001;
134 cf = (REAL) 0.0000001;
139 hy = ((REAL) 1.0) / cf;
147 // CALCULATE THE INTERSECTION OF REFERANCE LINE AND THE RAY IN QUESTI
149 sory += elem / ((REAL) 2.0);
153 // REFLECT IN Y=0 IF NEEDED
157 sorx += elem / ((REAL) 2.0);
158 c = sory - sorx * tanphi;
160 // FLIP LINE 180 DEGREES
162 c = (((REAL) 1.0) - tanphi) * elem - c;
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
173 // IF Y INTERCEPT IS NEGATIVE THEN TAKE X INTERCEPT (WHICH MUST
174 // BE POSITIVE SINCE SLOPE IS ALWAYS POSITIVE
176 // RAY ENTERS FROM BOTTOM
180 // RECALCULATE C SO THAT EQUATION OF LINE IS NOW
182 // IF X INTERCEPT IS GREATER TNAN N THE RAY HAS MISSED THE RECONSTRUC
188 // INITIAL GRID (IX,IY) COORDINATE IS NOW CALCULATED
193 // FP IS ACTUAL DISPLACEMENT OF INTERCEPT FROM THE NEXT GREATER
196 fp = ((REAL) (ix + 1)) - c;
198 // START AT LOW X VALUES AND INCREMENT ACROSS PICTURE UNLESS
203 // FLIP INDICATOR IS EXACTLY WRONG FOR HERMANS MATRIX STORAGE
208 // START AT HIGH X VALUES AND DECREMENT (FROM RIGHT TO LEFT)
210 ix = GeoPar.nelem - ix - 1;
214 // SET INITIAL VALUES OF THE LENGTH OF THE RAY IN CURRENT CELL
220 // INDEX IS THE LINEAR ARRAY SUBSCRIPT OF THE PICTURE ARRAY
221 // INDEX INTO FORTRAN TYPE ARRAY
223 index = iy * GeoPar.nelem + ix + 1;
225 // INITIAL CONDITIONS SET FOR ENTRANCE TO LOOPS WE HAVE ENTERED
226 // THE PICTURE AREA FROM THE BOTTOM
233 // ENTER PICTURE FROM SIDE (LEFT)
234 // CALCULATE INITIAL GRID COORDINATES (IX,IY)
239 // START WITH LOW X VALUES AND GO UP UNLESS THE FLIP INDICATOR
241 // FLIP INDICATOR IS EXACTLY WRONG FOR HERMANS MATRIX STORAGE
245 // START WITH HIGH X VALUES AND DECREMENT
247 ix = GeoPar.nelem - 1;
253 // FP IS DISTANCE FROM Y INTERCEPT TO NEXT GREATER GRID LINE
255 fp = ((REAL) (iy + 1)) - c;
257 // CALCULATE ARRAY SUBSCRIPT IN PICTURE AREA OF INITIAL CELL
258 // INDEX INTO FORTRAN TYPE ARRAY
260 index = iy * GeoPar.nelem + ix + 1;
262 // CALCULATE INITIAL DESCRIPTION OF RAY IN FIRST CELL FOR LOOP
268 // RAY ENTERS FROM LEFT SIDE
278 // MAKE ROOM FOR NEW PICTURE ELEMENT
280 list[(*numb) - 1] = index - 1;
282 // STEP ONE IN X DIRECTION AND CORRESPONDING STEP IN Y DIRECTION
287 // IF XC IS NEGATIVE THE Y GRID LINE HAS BEEN CROSSED THE RAY
288 // HAS LEFT THROUGHTHE TOP
292 // ENTERED SIDE,LEFT THROUGH SIDE --ASSIGN MAXIMUM Y-TYPE WEIGHT
293 weight[(*numb) - 1] = hy;
300 weight[(*numb) - 1] = w;
303 if ((w - (REAL) 0.0000001) < 0)
316 //////////////////////////////////////////////////////////////
320 // RAY ENTERS FROM BOTTOM
321 // STEP ONE CELL IN Y DIRECTION AND TEST FOR COMPLETION
323 if ((iy - GeoPar.nelem) >= 0)
326 index += GeoPar.nelem;
331 //////////////////////////////////////////////////////////////
333 // NEW LIST ELEMENT MEANS ANOTHER PIXEL SELECTED
336 list[(*numb) - 1] = index - 1;
338 // STEP ONE IN Y DIRECTION AND CORRESPONDING X INCREMENT
343 // IF YC IS NEGATIVE WH HAVE CROSSED A X-GRID LINE AND THAT
344 // WE LEAVE THROUGE THE SIDE
350 // ENTERED BOTTOM LEFT TOP ASSIGN MAXIMUM X-TYPE WEIGHT
352 weight[(*numb) - 1] = hx;
358 weight[(*numb) - 1] = w;
361 // IF THE WEIGHT SO CALCULATED IS TOO SMALL DELETE THIS PIXEL
362 // THIS CAUSED BY CROSSING GRID LINES AT A CORNER
364 if ((w - (REAL) 0.000001) <= 0.0)
377 /////////////////////////////////////////////////////////////
379 // ENTER PIXEL FROM SIDE
382 // STEP ONE IN X DIRECTION
383 // AND TEST FOR TERMINATION
388 if ((ix - GeoPar.nelem) >= 0)
402 fprintf(output, "\n wray np = %5i nr = %5i numb = %5i snorm = %9.4f", np, nr, *numb, *snorm);
404 if ((*numb > 0) && (trace > 9))
406 for (int i = 0; i < *numb; i++)
408 fprintf(output, "\n %5i %9.4f", list[i], weight[i]);