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/raylen.cpp $
5 $LastChangedRevision: 91 $
6 $Date: 2014-07-02 17:34:27 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 RAYLEN RETURNS THE LENGTH OF THE RAY IN AN ELEMENTAL OBJECT
12 SEE SNARK05 MANUAL FOR A FULL DESCRIPTION OF THE FUNCTIONING OF
15 TYPE = 1 FOR ELLIPSE/CIRCLE
16 TYPE = 2 FOR RECTANGLE/SQUARE
17 TYPE = 3 FOR ISOS. TRIANGLE
18 TYPE = 4 FOR SEGMENT OF CIRCLE
19 TYPE = 5 FOR SECTOR OF CIRCLE
20 AX,AY AND MX,MY DEFINE THE RAY:
21 (AX,AY) COORDINATES OF ANY POINT ON THE RAY;
22 (MX,MY) ARE THE X- AND Y- DIRECTION COSINES OF THE RAY
23 IF THE RAY MAKES AN ANGLE THETA WITH X-AXIS THEN MX = COS(THETA)
26 (CX,CY) ARE THE COORDINATES OF THE CENTER OF THE OBJECT
28 (U,V) ARE THE HALF- AXIS LENGTHS OF THE OBJECT
30 (PX,PY) ARE THE COSINE AND SINE OF THE ANGLE OF ROTATION OF THE
31 OBJECT ( THE ANGLE U MAKES WITH THE X-AXIS)
42 REAL raylen(SNARK_Object_type type, REAL ax, REAL ay, REAL mx, REAL my, REAL cx,
43 REAL cy, REAL u, REAL v, REAL px, REAL py)
48 REAL xnew, ynew, unew, vnew, xdc, ydc;
66 // SHIFT THE ORIGIN TO THE CENTER OF THE OBJECT
71 // ROTATE ABOUT ORIGIN BY AN ANGLE = (U,X-AXIS) SO THAT U AND V ARE
72 // ALONG THE NEW X AND Y AXES RESPECTIVELY
74 xnew = xo * px + yo * py;
75 ynew = -xo * py + yo * px;
78 xdc = mx * px + my * py;
79 ydc = -mx * py + my * px;
82 // FIND THE CIRCLE COMPLETELY ENCLOSING OBJECT
84 diam = (REAL) 2.1 * MAX0(unew, vnew);
86 // CHECK IF RAY HAS SLOPE .LE. 1.0 FOR COMPUTATIONAL PURPOSES
87 // REMEMBER THIS FOR TRIA, SEGM, AND SECT SINCE THEY DO NOT
88 // HAVE TWO-FOLD SYMMETRY
91 if (!((REAL) fabs(xdc) >= (REAL) fabs(ydc)))
94 // ROTATE ABOUT ORIGIN BY 90 DEG TO INSURE THAT SLOPE IS LE 1.0 IN
111 // SEE IF RAY IS ANYWHERE NEAR OBJECT CIRCLE
114 if ((REAL) fabs(slope) <= (REAL) 0.000001)
119 ycept = ynew - slope * xnew;
121 if ((REAL) fabs(ycept) < (REAL) 0.000001)
128 if ((REAL) fabs(slope) > Consts.zero)
130 xcept = xnew - ynew / slope;
133 if (((REAL) fabs(xcept) > diam) && ((REAL) fabs(ycept) > diam))
138 // if RECTANGLE or CIRCLE/ELLIPSE
140 if (type == SOT_rect)
142 // CHECK IF RAY THROUGH RECTANGLE/SQUARE IS NEARLY PARALLEL TO X AXIS
143 if (!((REAL) fabs(ydc) > (REAL) 0.0001))
145 if ((REAL) fabs(ynew) <= vnew)
146 raylen_temp = (REAL) 2.0 * unew;
150 // GENERAL LINE INSIDE A RECTANGLE
154 // there is still difference here between fortran
155 raylen_temp = (REAL) fabs(
157 MAX0(-unew, MIN0(unew, ((REAL) (vnew - ynew) / slope) + xnew))
158 - MAX0(-unew, MIN0(unew, ((REAL) -(vnew + ynew) / slope) + xnew))
167 MAX0(-unew, MIN0(unew, (vnew - ynew) / slope + xnew))
169 MIN0(unew, -(vnew + ynew) / slope + xnew)))
176 // GENERAL LINE INSIDE ELLIPSE
177 if (type == SOT_elip)
179 vsu = vnew * vnew + (unew * slope) * (unew * slope);
180 d = vsu - ycept * ycept;
182 // IF D POSITIVE, RAY INTERSECTS THE CIRCLE/ELLIPSE
185 raylen_temp = (REAL) 2.0 * unew * vnew * (REAL) sqrt(d)
186 / (REAL) fabs(vsu * xdc);
188 REAL zzz = (REAL) sqrt(d);
189 REAL xxx = (REAL) 2.0 * unew * vnew * zzz;
190 REAL yyy = (REAL) fabs(vsu * xdc);
191 REAL aaa = (REAL) xxx / yyy;
192 raylen_temp = (REAL) xxx / yyy;
198 // TRIANGLE, SEGMENT, SECTOR
200 if (!(type == SOT_segm))
203 // TRIANGLE AND TRIANGULAR PORTION OF SECTOR
205 // SPECIAL CASE FOR TRIANGLE AND SECTOR
206 // CHECK IF RAY IS NEARLY PARALLEL TO ONE OF THE SIDES
209 if (!((REAL) fabs(vtou - (REAL) fabs(slope)) >= (REAL) 0.0001))
213 if (!((REAL) fabs(xcept) > unew))
216 raylen_temp = (REAL) 0.5
217 * (REAL) fabs((unew - xcept) / xdc);
219 raylen_temp = (REAL) 0.5
220 * (REAL) fabs((unew + xcept) / xdc);
225 if (!((REAL) fabs(ycept) > vnew))
228 raylen_temp = (REAL) 0.5
229 * (REAL) fabs((vnew - ycept) / ydc);
231 raylen_temp = (REAL) 0.5
232 * (REAL) fabs((vnew + ycept) / ydc);
238 // TRIANGLE GENERAL CASE
242 x1 = (vnew - ycept) / (slope + vtou);
243 if ((x1 <= -Consts.zero) || (x1 > unew)) //bug 273, jklukowska
245 x2 = (vnew - ycept) / (slope - vtou);
246 if ((x2 < -unew) || (x2 >= Consts.zero)) //bug 273, jklukowska
249 // RAY THROUGH VERTEX OF ISOSCELES TRIANGLE
250 if (((REAL) fabs(x1) < (REAL) 0.00001)
251 && ((REAL) fabs(x2) < (REAL) 0.00001)
252 && ((REAL) fabs(xcept) <= unew))
254 raylen_temp = (REAL) fabs((x1 - x2) / xdc);
259 x1 = (vnew - ycept) / (slope + vtou);
260 if ((x1 < 0.0) || (x1 >= unew + 0.00001))
262 x2 = (-vnew - ycept) / (slope - vtou);
263 if ((x2 < 0.0) || (x2 >= unew + 0.00001))
265 // RAY THROUGH VERTEX OF ISOSCELES TRIANGLE
267 if (((REAL) fabs(x1 - unew) < (REAL) 0.000001)
268 && ((REAL) fabs(x2 - unew) < (REAL) 0.000001)
269 && ((REAL) fabs(ycept) <= vnew))
271 raylen_temp = (REAL) fabs((x1 - x2) / xdc);
276 if (type == SOT_tria)
282 // SEGMENT OR THE SEGMENT PORTION OF THE SECTOR OF A CIRCLE
286 if (((REAL) fabs(xcept) > unew) && (ycept > 0.0))
291 a = (REAL) 1.0 + slope * slope;
294 c = bb * bb - (unew * unew + vnew * vnew);
301 // THE RAY DOES INTERSECT THE CIRCLE
308 if (x1 * slope + ycept > 0.0)
310 if (x2 * slope + ycept > 0.0)
313 // X COMPONENT OF INTERSECTION LENGTH IS DIFFERENCE OF ROOTS OF QUADR
314 if (raylen_temp > 0.0)
318 raylen_temp += (REAL) 2.0 * d / (REAL) fabs(a * xdc);
321 // RAY INTERSECTS LINEAR BOUNDARY OF SEGMENT
323 if (x2 * slope + ycept <= 0.0)
328 L113: raylen_temp += (REAL) fabs((x1 - x2) / xdc);
332 if (((REAL) fabs(ycept) > vnew) && (xcept > 0.0))
337 a = (REAL) 1.0 + slope * slope;
338 b = slope * ycept - unew;
339 c = ycept * ycept - vnew * vnew;
358 // X COMPONENT OF INTERSECTION LENGTH IS DIFFERENCE OF ROOTS OF QUADR
359 if (raylen_temp > 0.0)
363 raylen_temp += (REAL) 2.0 * d / (REAL) fabs(a * xdc);
366 // RAY INTERSECTS LINEAR BOUNDARY OF SEGMENT
373 L123: raylen_temp += (REAL) fabs((x1 - x2) / xdc);