81 |
81 |
82 root->subdivide(bbox, max_depth); |
82 root->subdivide(bbox, max_depth); |
83 built = true; |
83 built = true; |
84 } |
84 } |
85 |
85 |
86 static inline int first_node(const Float tx0, const Float ty0, const Float tz0, |
86 |
87 const Float txm, const Float tym, const Float tzm) |
87 /******************************************************* |
88 { |
88 octree traversal algorithm as described in paper |
89 int res = 0; |
89 "An Efficient Parametric Algorithm for Octree Traversal" |
90 if (tx0 > ty0) |
90 by J. Revelles, C. Urena and M. Lastra. |
91 { |
91 |
92 if (tx0 > tz0) |
92 see revision 37 for original recursive version |
93 { // YZ |
93 *******************************************************/ |
94 if (tym < tx0) |
94 |
95 res |= 2; |
95 struct OctreeTravState |
96 if (tzm < tx0) |
96 { |
97 res |= 1; |
97 Float tx0,ty0,tz0,tx1,ty1,tz1,txm,tym,tzm; |
98 } |
98 OctreeNode *node; |
99 else |
99 int next; |
100 { // XY |
100 OctreeTravState() {}; |
101 if (txm < tz0) |
101 OctreeTravState( |
102 res |= 4; |
102 const Float atx0, const Float aty0, const Float atz0, |
103 if (tym < tz0) |
103 const Float atx1, const Float aty1, const Float atz1, |
104 res |= 2; |
104 const Float atxm, const Float atym, const Float atzm, |
105 } |
105 OctreeNode *const anode, const int anext): |
106 } |
106 tx0(atx0), ty0(aty0), tz0(atz0), tx1(atx1), ty1(aty1), tz1(atz1), |
107 else |
107 txm(atxm), tym(atym), tzm(atzm), node(anode), next(anext) {}; |
108 { |
108 }; |
109 if (ty0 > tz0) |
109 |
110 { // XZ |
110 inline const int &next_node(const Float &txm, const int &xnode, |
111 if (txm < ty0) |
111 const Float &tym, const int &ynode, const Float &tzm, const int &znode) |
112 res |= 4; |
|
113 if (tzm < ty0) |
|
114 res |= 1; |
|
115 return res; |
|
116 } |
|
117 else |
|
118 { // XY |
|
119 if (txm < tz0) |
|
120 res |= 4; |
|
121 if (tym < tz0) |
|
122 res |= 2; |
|
123 } |
|
124 } |
|
125 return res; |
|
126 } |
|
127 |
|
128 static inline int next_node(const Float txm, const int xnode, |
|
129 const Float tym, const int ynode, const Float tzm, const int znode) |
|
130 { |
112 { |
131 if (txm < tym) |
113 if (txm < tym) |
132 { |
114 { |
133 if (txm < tzm) |
115 if (txm < tzm) |
134 return xnode; |
116 return xnode; |
142 else |
124 else |
143 return znode; |
125 return znode; |
144 } |
126 } |
145 } |
127 } |
146 |
128 |
147 static Shape *proc_subtree(const int a, const Float tx0, const Float ty0, const Float tz0, |
|
148 const Float tx1, const Float ty1, const Float tz1, OctreeNode *node, |
|
149 const Shape *origin_shape, const Ray &ray, Float &nearest_distance) |
|
150 { |
|
151 Float txm, tym, tzm; |
|
152 int curr_node; |
|
153 |
|
154 // if ray does not intersect this node |
|
155 if (tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0) |
|
156 return NULL; |
|
157 |
|
158 if (node->isLeaf()) |
|
159 { |
|
160 Shape *nearest_shape = NULL; |
|
161 ShapeList::iterator shape; |
|
162 Float mindist = max(max(tx0,ty0),tz0); |
|
163 Float dist = min(min(min(tx1,ty1),tz1),nearest_distance); |
|
164 for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++) |
|
165 if (*shape != origin_shape && (*shape)->intersect(ray, dist) && dist >= mindist) |
|
166 { |
|
167 nearest_shape = *shape; |
|
168 nearest_distance = dist; |
|
169 } |
|
170 return nearest_shape; |
|
171 } |
|
172 |
|
173 txm = 0.5 * (tx0+tx1); |
|
174 tym = 0.5 * (ty0+ty1); |
|
175 tzm = 0.5 * (tz0+tz1); |
|
176 |
|
177 curr_node = first_node(tx0,ty0,tz0,txm,tym,tzm); |
|
178 Shape *shape = NULL; |
|
179 while (curr_node < 8) |
|
180 { |
|
181 switch (curr_node) |
|
182 { |
|
183 case 0: |
|
184 shape =proc_subtree (a,tx0,ty0,tz0,txm,tym,tzm,node->getChild(a), origin_shape, ray, nearest_distance); |
|
185 curr_node = next_node(txm, 4, tym, 2, tzm, 1); |
|
186 break; |
|
187 case 1: |
|
188 shape =proc_subtree (a,tx0,ty0,tzm,txm,tym,tz1,node->getChild(1^a), origin_shape, ray, nearest_distance); |
|
189 curr_node = next_node(txm, 5, tym, 3, tz1, 8); |
|
190 break; |
|
191 case 2: |
|
192 shape =proc_subtree (a,tx0,tym,tz0,txm,ty1,tzm,node->getChild(2^a), origin_shape, ray, nearest_distance); |
|
193 curr_node = next_node(txm, 6, ty1, 8, tzm, 3); |
|
194 break; |
|
195 case 3: |
|
196 shape =proc_subtree (a,tx0,tym,tzm,txm,ty1,tz1,node->getChild(3^a), origin_shape, ray, nearest_distance); |
|
197 curr_node = next_node(txm, 7, ty1, 8, tz1, 8); |
|
198 break; |
|
199 case 4: |
|
200 shape =proc_subtree (a,txm,ty0,tz0,tx1,tym,tzm,node->getChild(4^a), origin_shape, ray, nearest_distance); |
|
201 curr_node = next_node(tx1, 8, tym, 6, tzm, 5); |
|
202 break; |
|
203 case 5: |
|
204 shape =proc_subtree (a,txm,ty0,tzm,tx1,tym,tz1,node->getChild(5^a), origin_shape, ray, nearest_distance); |
|
205 curr_node = next_node(tx1, 8, tym, 7, tz1, 8); |
|
206 break; |
|
207 case 6: |
|
208 shape =proc_subtree (a,txm,tym,tz0,tx1,ty1,tzm,node->getChild(6^a), origin_shape, ray, nearest_distance); |
|
209 curr_node = next_node(tx1, 8, ty1, 8, tzm, 7); |
|
210 break; |
|
211 case 7: |
|
212 shape =proc_subtree (a,txm,tym,tzm,tx1,ty1,tz1,node->getChild(7^a), origin_shape, ray, nearest_distance); |
|
213 curr_node = 8; |
|
214 break; |
|
215 } |
|
216 if (shape != NULL) |
|
217 return shape; |
|
218 } |
|
219 return NULL; |
|
220 } |
|
221 |
|
222 /* |
|
223 traversal algorithm paper as described in paper |
|
224 "An Efficient Parametric Algorithm for Octree Traversal" |
|
225 by J. Revelles, C. Urena and M. Lastra. |
|
226 */ |
|
227 Shape * Octree::nearest_intersection(const Shape *origin_shape, const Ray &ray, |
129 Shape * Octree::nearest_intersection(const Shape *origin_shape, const Ray &ray, |
228 Float &nearest_distance) |
130 Float &nearest_distance) |
229 { |
131 { |
230 /* if we have no tree, fall back to naive test */ |
132 /* if we have no tree, fall back to naive test */ |
231 if (!built) |
133 if (!built) |
232 return Container::nearest_intersection(origin_shape, ray, nearest_distance); |
134 return Container::nearest_intersection(origin_shape, ray, nearest_distance); |
233 |
135 |
|
136 OctreeTravState st[max_depth+1]; |
|
137 register OctreeTravState *st_cur = st; |
|
138 |
|
139 # define node st_cur->node |
|
140 # define tx0 st_cur->tx0 |
|
141 # define ty0 st_cur->ty0 |
|
142 # define tz0 st_cur->tz0 |
|
143 # define tx1 st_cur->tx1 |
|
144 # define ty1 st_cur->ty1 |
|
145 # define tz1 st_cur->tz1 |
|
146 # define txm st_cur->txm |
|
147 # define tym st_cur->tym |
|
148 # define tzm st_cur->tzm |
|
149 |
234 int a = 0; |
150 int a = 0; |
235 Vector3 ro = ray.o; |
151 Vector3 ro(ray.o); |
236 Vector3 rdir = ray.dir; |
152 Vector3 rdir(1.0/ray.dir.x, 1.0/ray.dir.y, 1.0/ray.dir.z); |
237 |
153 |
238 if (rdir.x < 0.0) |
154 if (rdir.x < 0.0) |
239 { |
155 { |
240 ro.x = (bbox.L.x+bbox.H.x) - ro.x; |
156 ro.x = (bbox.L.x+bbox.H.x) - ro.x; |
241 rdir.x = -rdir.x; |
157 rdir.x = -rdir.x; |
251 { |
167 { |
252 ro.z = (bbox.L.z+bbox.H.z) - ro.z; |
168 ro.z = (bbox.L.z+bbox.H.z) - ro.z; |
253 rdir.z = -rdir.z; |
169 rdir.z = -rdir.z; |
254 a |= 1; |
170 a |= 1; |
255 } |
171 } |
256 Float tx0 = (bbox.L.x - ro.x) / rdir.x; |
172 |
257 Float tx1 = (bbox.H.x - ro.x) / rdir.x; |
173 tx0 = (bbox.L.x - ro.x) * rdir.x; |
258 Float ty0 = (bbox.L.y - ro.y) / rdir.y; |
174 tx1 = (bbox.H.x - ro.x) * rdir.x; |
259 Float ty1 = (bbox.H.y - ro.y) / rdir.y; |
175 ty0 = (bbox.L.y - ro.y) * rdir.y; |
260 Float tz0 = (bbox.L.z - ro.z) / rdir.z; |
176 ty1 = (bbox.H.y - ro.y) * rdir.y; |
261 Float tz1 = (bbox.H.z - ro.z) / rdir.z; |
177 tz0 = (bbox.L.z - ro.z) * rdir.z; |
262 |
178 tz1 = (bbox.H.z - ro.z) * rdir.z; |
263 if (max(max(tx0,ty0),tz0) < min (min(tx1,ty1),tz1)) |
179 |
264 return proc_subtree(a,tx0,ty0,tz0,tx1,ty1,tz1,root, |
180 if (max3(tx0,ty0,tz0) > min3(tx1,ty1,tz1)) |
265 origin_shape, ray, nearest_distance); |
|
266 else |
|
267 return NULL; |
181 return NULL; |
268 } |
182 |
|
183 node = root; |
|
184 st_cur->next = -1; |
|
185 |
|
186 Shape *nearest_shape = NULL; |
|
187 while (nearest_shape == NULL) |
|
188 { |
|
189 if (st_cur->next == -1) |
|
190 { |
|
191 st_cur->next = 8; |
|
192 |
|
193 // if ray does intersect this node |
|
194 if (!(tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0)) |
|
195 { |
|
196 if (node->isLeaf()) |
|
197 { |
|
198 ShapeList::iterator shape; |
|
199 register Float mindist = max3(tx0,ty0,tz0); |
|
200 /* correct & slow */ |
|
201 //Float dist = min(min3(tx1,ty1,tz1),nearest_distance); |
|
202 /* faster */ |
|
203 register Float dist = nearest_distance; |
|
204 for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++) |
|
205 if (*shape != origin_shape && (*shape)->intersect(ray, dist) && dist >= mindist) |
|
206 { |
|
207 nearest_shape = *shape; |
|
208 nearest_distance = dist; |
|
209 } |
|
210 } |
|
211 else |
|
212 { |
|
213 txm = 0.5 * (tx0+tx1); |
|
214 tym = 0.5 * (ty0+ty1); |
|
215 tzm = 0.5 * (tz0+tz1); |
|
216 |
|
217 // first node |
|
218 st_cur->next = 0; |
|
219 if (tx0 > ty0) |
|
220 { |
|
221 if (tx0 > tz0) |
|
222 { // YZ |
|
223 if (tym < tx0) |
|
224 st_cur->next |= 2; |
|
225 if (tzm < tx0) |
|
226 st_cur->next |= 1; |
|
227 } |
|
228 else |
|
229 { // XY |
|
230 if (txm < tz0) |
|
231 st_cur->next |= 4; |
|
232 if (tym < tz0) |
|
233 st_cur->next |= 2; |
|
234 } |
|
235 } |
|
236 else |
|
237 { |
|
238 if (ty0 > tz0) |
|
239 { // XZ |
|
240 if (txm < ty0) |
|
241 st_cur->next |= 4; |
|
242 if (tzm < ty0) |
|
243 st_cur->next |= 1; |
|
244 } |
|
245 else |
|
246 { // XY |
|
247 if (txm < tz0) |
|
248 st_cur->next |= 4; |
|
249 if (tym < tz0) |
|
250 st_cur->next |= 2; |
|
251 } |
|
252 } |
|
253 } |
|
254 } |
|
255 } |
|
256 |
|
257 while (st_cur->next == 8) |
|
258 { |
|
259 // pop state from stack |
|
260 if (st_cur == st) |
|
261 return NULL; // nothing to pop, finish |
|
262 --st_cur; |
|
263 } |
|
264 |
|
265 // push current state |
|
266 *(st_cur+1) = *st_cur; |
|
267 ++st_cur; |
|
268 |
|
269 switch (st_cur->next) |
|
270 { |
|
271 case 0: |
|
272 tx1 = txm; |
|
273 ty1 = tym; |
|
274 tz1 = tzm; |
|
275 node = node->getChild(a); |
|
276 (st_cur-1)->next = next_node(txm, 4, tym, 2, tzm, 1); |
|
277 break; |
|
278 case 1: |
|
279 tz0 = tzm; |
|
280 tx1 = txm; |
|
281 ty1 = tym; |
|
282 node = node->getChild(1^a); |
|
283 (st_cur-1)->next = next_node(txm, 5, tym, 3, tz1, 8); |
|
284 break; |
|
285 case 2: |
|
286 ty0 = tym; |
|
287 tx1 = txm; |
|
288 tz1 = tzm; |
|
289 node = node->getChild(2^a); |
|
290 (st_cur-1)->next = next_node(txm, 6, ty1, 8, tzm, 3); |
|
291 break; |
|
292 case 3: |
|
293 ty0 = tym; |
|
294 tz0 = tzm; |
|
295 tx1 = txm; |
|
296 node = node->getChild(3^a); |
|
297 (st_cur-1)->next = next_node(txm, 7, ty1, 8, tz1, 8); |
|
298 break; |
|
299 case 4: |
|
300 tx0 = txm; |
|
301 ty1 = tym; |
|
302 tz1 = tzm; |
|
303 node = node->getChild(4^a); |
|
304 (st_cur-1)->next = next_node(tx1, 8, tym, 6, tzm, 5); |
|
305 break; |
|
306 case 5: |
|
307 tx0 = txm; |
|
308 tz0 = tzm; |
|
309 ty1 = tym; |
|
310 node = node->getChild(5^a); |
|
311 (st_cur-1)->next = next_node(tx1, 8, tym, 7, tz1, 8); |
|
312 break; |
|
313 case 6: |
|
314 tx0 = txm; |
|
315 ty0 = tym; |
|
316 tz1 = tzm; |
|
317 node = node->getChild(6^a); |
|
318 (st_cur-1)->next = next_node(tx1, 8, ty1, 8, tzm, 7); |
|
319 break; |
|
320 case 7: |
|
321 tx0 = txm; |
|
322 ty0 = tym; |
|
323 tz0 = tzm; |
|
324 node = node->getChild(7^a); |
|
325 (st_cur-1)->next = 8; |
|
326 break; |
|
327 } |
|
328 st_cur->next = -1; |
|
329 } |
|
330 return nearest_shape; |
|
331 } |