1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
|
/*
* $Id: track.h,v 1.1 2005-09-18 22:04:56 dhmunro Exp $
* Routines for tracking a straight ray in 3D through a cylindrical mesh.
*/
/* Copyright (c) 2005, The Regents of the University of California.
* All rights reserved.
* This file is part of yorick (http://yorick.sourceforge.net).
* Read the accompanying LICENSE file for details.
*/
#ifndef TRACK_H
#define TRACK_H
#include "bound.h"
/* ---------------------------------------------------------------------- */
/* The result of a ray tracking operation is a RayPath: */
typedef struct RayPath RayPath;
struct RayPath {
long maxcuts, ncuts; /* physical lengths, current number zones cut */
long *zone; /* [ncuts] array of zone entered, if last exit
* (into region 0) then zone=0 */
double *ds; /* [ncuts] path length in zone (0.0 if zone=0) */
long *pt1, *pt2; /* [ncuts] arrays, pt1->pt2 is edge of zone cut */
double *f; /* [ncuts] cut at pt1 + (f+0.5)*(pt2-pt1) */
double fi,ff; /* fraction of distance into initial (final)
* zone at which tracking started (stopped)
* ds in initial (final) zone is 1-fi (1-ff) of
* the full length of the ray within the zone */
};
/* Although the mesh lives in a 2-D cylindrical coordinate space, a Ray
* is a directed line in Cartesian 3-D space. The z-axis is common to
* the two coordinate systems. The x-axis coincides with the r-axis
* in a standard (z,r) mesh plot (although, of course, r>=0, while x
* can have either sign). The y-axis comes out of the page of the
* standard mesh plot, toward your face. Hence, (x,y,z) is a right
* handed coordinate frame.
*
* A ray is assumed to lie in a plane of constant y. The angle theta
* is the angle from the +z direction to the ray direction, with
* theta=+pi/2 corresponding to the +x-direction. (This is opposite
* the convention in TDG/DIRT, but corresponds to phi=0 azimuth in
* the usual spherical coordinate system, e.g.- Jackson's E&M book.)
*
* The redundant data in the Ray structure allows the ray tracker to
* step from edge to edge without having to compute a lot of extra
* square roots, e.g.- to get r=sqrt(x^2+y^2) or sin=sqrt(1-cos^2).
* (The alternative approach, used in TDG/DIRT, is to describe the ray
* in an invariant fashion, i.e.- independent of the current zone in
* the tracking process.)
*/
typedef struct Ray Ray;
struct Ray {
double cos, sin; /* cos(theta), sin(theta) */
double y, z, x; /* any point on the ray...
* y remains fixed,
* while z and x step from edge to edge */
double r; /* (z,r) are cylindrical coordinates corresponding
* to (x,y,z) -- r^2 = x^2 + y^2 */
};
/* Structure to store results of a ray-edge intersection calculation
* made by ExitEdge. There are potentially two solutions, the "exit"
* solution, representing a left-to-right crossing in (z,r) space, and
* the "entry" solution, representing a right-to-left crossing.
* Neither solution necessarily exists. Several intermediate results
* of the calculation are saved here, as well as the two roots.
*/
typedef struct RayEdgeInfo RayEdgeInfo;
struct RayEdgeInfo {
double dz, dr; /* edge vector in (z,r) plane */
double area; /* a useful cross product (see EdgeExit) */
double A, B, C; /* Af^2+Bf+C=0 is equation for fx, fn --
* neither B nor C valid unless D>0 */
double D; /* sqrt(discriminant) if >0, else discriminant */
double fx; /* the left-to-right crossing (edge fraction) */
int validx; /* fx not valid unless this is true */
double fn; /* the right-to-left crossing (edge fraction) */
int validn; /* fn not valid unless this is true */
};
/* A ray enters the mesh to begin tracking at one or more EntryPoint's: */
typedef struct EntryPoint EntryPoint;
struct EntryPoint {
EntryPoint *next;
Ray ray;
RayEdgeInfo info;
long zone;
int side;
double f;
double s0;
};
/* ---------------------------------------------------------------------- */
/* Given the mesh boundary (i.e.- list of boundary edges) and a ray,
* return a linked list of entry point(s), ordered by increasing time.
*/
extern EntryPoint *FindEntryPoints(Boundary *boundary, Ray *ray);
extern void FreeEntryPoints(EntryPoint *entry);
/* Starting at the given entry point(s), track the ray through the mesh.
* The path arrays will be lengthened if necessary, but never shortened.
* sLimits[0] <= s[1]<=s[ncuts-2] <= sLimits[1], if sLimits[0]<sLimits[1]
* on entry. (Note that s[0] and s[ncuts-1] may lie outside range.)
*/
extern void RayTrack(Mesh *mesh, EntryPoint *entry, RayPath *path,
double *sLimits);
/* Track the ray through the mesh, assuming it represents a sphere.
* The path arrays will be lengthened if necessary, but never shortened.
* sLimits[0] <= s[1]<=s[ncuts-2] <= sLimits[1], if sLimits[0]<sLimits[1]
* on entry. (Note that s[0] and s[ncuts-1] may lie outside range.)
*/
extern void RayTrackS(Mesh *mesh, Ray *ray, RayPath *path, double *sLimits);
/* Find the intersections of a ray with an edge (z,r); return true if
* and only if the ray cuts the edge segment from left to right.
*/
extern int
ExitEdge(Ray *ray, double z[2], double r[2], /* inputs */
int *after, /* int *notafter, */ /* updates */
RayEdgeInfo *info); /* output */
/* Compute the path length from the current point on the ray to the
* exit (left-to-right) intersection computed by ExitEdge.
*/
extern double RayPathLength(Ray *ray, RayEdgeInfo *info);
/* Compute the path length from the exit point to the entry point (>0
* if entry AFTER exit) for the intersections computed by ExitEdge.
*/
extern double RayPathDifference(RayEdgeInfo *info);
/* Given a ray with current point entering a particular zone and side of
* the mesh, find the corresponding exit. The path length ds is >=0
* unless the ray segment lies in the negative area part of a bowtied
* zone, in which case ds<0. Return the exit side index (0<=side<=3).
* Update the current point on the ray to the exit point.
*/
extern int
ExitZone(Mesh *mesh, long zone, /* input mesh and zone number */
int side, /* input side number on entry */
Ray *ray, /* updated to exit point */
RayEdgeInfo *info[4], /* info[3] is entry side on entry,
* updated to exit side on exit */
double *ds, /* output ray path length */
double *f); /* output edge fraction */
/* Although algebraically exact formulas are used for the ExitEdge
* calculation, roundoff errors may make the redundant information
* in the Ray structure (y,z,x, and r) inconsistent. PolishExit
* attempts to restore precise self-consistency of the ray, without
* moving the point off of the current edge (remembered in info).
*/
extern void PolishExit(Ray *ray, RayEdgeInfo *info, double *ds, double *f);
/* Through roundoff errors, it is possible for the exit point found
* in ExitZone (or FindEntryPoints) to lie slightly beyond the
* endpoints of the exit edge. AdjustRayXY moves the ray back to
* the nearby endpoint when this occurs.
*/
extern void AdjustRayXY(Ray *ray, double* z, double *r);
/* Last ditch effort by ExitZone to find the exit edge */
extern int FindLostRay(Ray *ray, RayEdgeInfo *info[4],
double z[4], double r[4], double dsx[4]);
/* Rearrange the linked list of entry points into increasing order */
extern EntryPoint *EntrySort(EntryPoint *entry);
/* Make ray path arrays longer */
extern void ExtendRayPath(RayPath *path, long pathinc);
/* Free all pointers in a RayPath (but NOT the input path pointer). */
extern void EraseRayPath(RayPath *path);
/* Adjust point on ray (x,y,z) and (z,r) to be point such that normal
* plane to ray passes through origin.
*/
void NormalizeRay(Ray *ray);
/* ---------------------------------------------------------------------- */
/********** DEFAULT VALUES SET IN track.c ********/
/* Global variables control the root polishing routine PolishExit */
extern int polishRoot;
extern double polishTol1;
extern double polishTol2;
/* Global variable controlling roundoff tolerance in FindLostRay */
/* Setting this too large could result in an infinite loop, but in
* any such situation, FindLostRay will be unable to find the ray
* no matter what the setting of findRayTol */
extern double findRayTol;
/* ---------------------------------------------------------------------- */
#endif
|
