FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_vis_ray_trace.c
Go to the documentation of this file.
1/*****************************************************************************
2 * Copyright (c) 2019 FrontISTR Commons
3 * This software is released under the MIT License, see LICENSE.txt
4 *****************************************************************************/
5
7
8#include <math.h>
11
12#if defined(old_version) || defined(old_intersection)
14#endif
15
16/*
17void transform_face_node(int ijkn[3], int current_ijk[3], int face_sect[3],
18double vv[3],
19 int next_ijk[3]);
20 */
21static int find_out_point(double orig_xyz[3], double r_dxyz[3], int ijk[3],
22 double in_point[3], double ray_direction[3],
23 int face_sect[3], double out_point[3]);
24static void find_next_cell(int ijkn[3], int current_ijk[3], int face_sect[3],
25 double vv[3], int next_ijk[3]);
26
27static int find_first_inter_point(double point_o[3], double view_p[3],
28 double ray_direction[3], double orig_xyz[3],
29 double dxyz[3], int r_level[3],
30 double first_p[3]) {
31 int i, j, mincomp, intersection;
32 double minmax[6], t[6], mint;
33
34 for (i = 0; i < 3; i++) {
35 minmax[i * 2] = orig_xyz[i];
36 minmax[i * 2 + 1] = orig_xyz[i] + dxyz[i];
37 }
38 for (i = 0; i < 3; i++) {
39 if (fabs(ray_direction[i]) < EPSILON) {
40 t[i * 2] = 1.0E+17;
41 t[i * 2 + 1] = 1.0E+17;
42 } else {
43 t[i * 2] = (minmax[i * 2] - view_p[i]) / ray_direction[i];
44 t[i * 2 + 1] = (minmax[i * 2 + 1] - view_p[i]) / ray_direction[i];
45 }
46 }
47#ifdef old_version
48 for (i = 0; i < 3; i++) {
49 t[i * 2] = (minmax[i * 2] - view_p[i]) / (ray_direction[i] + EPSILON);
50 t[i * 2 + 1] =
51 (minmax[i * 2 + 1] - view_p[i]) / (ray_direction[i] + EPSILON);
52 }
53#endif
54
55 mint = 1.0E+17;
56 mincomp = -1;
57 for (i = 0; i < 6; i++) {
58 if ((mint > t[i]) && (t[i] > EPSILON)) {
59 for (j = 0; j < 3; j++) first_p[j] = view_p[j] + t[i] * ray_direction[j];
60 if ((first_p[0] >= minmax[0] - EPSILON) &&
61 (first_p[0] <= minmax[1] + EPSILON) &&
62 (first_p[1] >= minmax[2] - EPSILON) &&
63 (first_p[1] <= minmax[3] + EPSILON) &&
64 (first_p[2] >= minmax[4] - EPSILON) &&
65 (first_p[2] <= minmax[5] + EPSILON)) {
66 mint = t[i];
67 mincomp = i;
68 }
69 }
70 }
71 if ((mincomp >= 0) && (mincomp <= 5)) {
72 intersection = 1;
73 for (i = 0; i < 3; i++) first_p[i] = view_p[i] + mint * ray_direction[i];
74 } else {
75 intersection = 0;
76 }
77 /* for(i=0;i<3;i++) {
78 if(fabs(out_point[i]-minmax[i*2])<EPSILON)
79 face_sect[i]=i*2;
80 else if(fabs(out_point[i]-minmax[i*2+1])<EPSILON)
81 face_sect[i]=i*2+1;
82 }
83 */
84
85 return (mincomp);
86}
87
88#ifdef old_version
89
90int find_first_inter_point(double point_o[3], double view_p[3], VR_data *vd,
91 double first_p[3]) {
92 int i, j, m, node[4], spm[6], minsp;
93 double vv[8 * 3], nn, a[6], b[6], c[6], d[6], abc, in_vv[3], in_point[3],
94 fp[4][3];
95 double t, st[6], mint, out_point[3], point[3];
96 int lc, con1;
97 int on_face_flag;
98
99 for (i = 0; i < 3; i++) in_vv[i] = point_o[i] - view_p[i];
100 nn = sqrt(SQR(in_vv[0]) + SQR(in_vv[1]) + SQR(in_vv[2]));
101 if (fabs(nn) > EPSILON) {
102 for (i = 0; i < 3; i++) in_vv[i] /= nn;
103 }
104 for (i = 0; i < 3; i++) in_point[i] = view_p[i];
105 vv[0 * 3] = vv[4 * 3] = vv[7 * 3] = vv[3 * 3] = vd->xyz0[0];
106 vv[1 * 3] = vv[5 * 3] = vv[6 * 3] = vv[2 * 3] =
107 vd->xyz0[0] + vd->nxyz[0] * vd->dxyz[0];
108 vv[0 * 3 + 2] = vv[1 * 3 + 2] = vv[2 * 3 + 2] = vv[3 * 3 + 2] = vd->xyz0[2];
109 vv[4 * 3 + 2] = vv[5 * 3 + 2] = vv[6 * 3 + 2] = vv[7 * 3 + 2] =
110 vd->xyz0[2] + vd->nxyz[2] * vd->dxyz[2];
111 vv[0 * 3 + 1] = vv[1 * 3 + 1] = vv[5 * 3 + 1] = vv[4 * 3 + 1] = vd->xyz0[1];
112 vv[3 * 3 + 1] = vv[2 * 3 + 1] = vv[6 * 3 + 1] = vv[7 * 3 + 1] =
113 vd->xyz0[1] + vd->nxyz[1] * vd->dxyz[1];
114
115 for (m = 0; m < 6; m++) {
116 transform_face_node(m, node);
117 for (i = 0; i < 4; i++)
118 for (j = 0; j < 3; j++) fp[i][j] = vv[node[i] * 3 + j];
119 a[m] = (fp[1][1] - fp[0][1]) * (fp[3][2] - fp[0][2]) -
120 (fp[3][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
121 b[m] = -(fp[1][0] - fp[0][0]) * (fp[3][2] - fp[0][2]) +
122 (fp[3][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
123 c[m] = (fp[1][0] - fp[0][0]) * (fp[3][1] - fp[0][1]) -
124 (fp[3][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
125 abc = sqrt(a[m] * a[m] + b[m] * b[m] + c[m] * c[m]);
126
127 a[m] /= abc;
128 b[m] /= abc;
129 c[m] /= abc;
130
131 d[m] = -(fp[0][0] * a[m] + fp[0][1] * b[m] + fp[0][2] * c[m]);
132 }
133 lc = -1;
134 con1 = 1;
135 for (m = 0; m < 6; m++) {
136 if (fabs(a[m] * in_vv[0] + b[m] * in_vv[1] + c[m] * in_vv[2]) > EPSILON) {
137 t = (-d[m] - a[m] * in_point[0] - b[m] * in_point[1] -
138 c[m] * in_point[2]) /
139 (a[m] * in_vv[0] + b[m] * in_vv[1] + c[m] * in_vv[2]);
140 /* printf("t=%f\n",t);*/
141 if (t > 0) {
142 for (j = 0; j < 3; j++) point[j] = in_point[j] + t * in_vv[j];
143 if ((point[0] > vd->xyz0[0] - EPSILON) &&
144 (point[0] < vd->xyz0[0] + vd->dxyz[0] * vd->nxyz[0] + EPSILON) &&
145 (point[1] > vd->xyz0[1] - EPSILON) &&
146 (point[1] < vd->xyz0[1] + vd->dxyz[1] * vd->nxyz[1] + EPSILON) &&
147 (point[2] > vd->xyz0[2] - EPSILON) &&
148 (point[2] < vd->xyz0[2] + vd->dxyz[2] * vd->nxyz[2] + EPSILON)) {
149 lc++;
150 st[lc] = t;
151 spm[lc] = m;
152 }
153
154 /* on_face_flag=0;
155if((m==0) || (m==1)) {
156if((point[1]>=vd->xyz0[1]-EPSILON)
157&& (point[1]<=vd->xyz0[1]+vd->dxyz[1]*vd->nxyz[1]+EPSILON) &&
158(point[2]>=vd->xyz0[2]-EPSILON)
159&& (point[2]<=vd->xyz0[2]+vd->dxyz[2]*vd->nxyz[2]+EPSILON))
160on_face_flag=1;
161}
162else if((m==2) || (m==3)) {
163if((point[0]>=vd->xyz0[0]-EPSILON) &&
164 (point[0]=<vd->xyz0[0]+vd->dxyz[0]*vd->nxyz[0]+EPSILON)
165&& (point[2]>=vd->xyz0[2]-EPSILON)
166&& (point[2]<=vd->xyz0[2]+vd->dxyz[2]*vd->nxyz[2]+EPSILON))
167 on_face_flag=1;
168}
169else if((m==4) || (m==5)) {
170if((point[0]>=vd->xyz0[0]-EPSILON) &&
171(point[0]=<vd->xyz0[0]+vd->dxyz[0]*vd->nxyz[0]+EPSILON)
172&& (point[1]>=vd->xyz0[1]-EPSILON)
173&& (point[1]<=vd->xyz0[1]+vd->dxyz[1]*vd->nxyz[1]+EPSILON))
174on_face_flag=1;
175}
176if(on_face_flag==1) {
177lc++;
178st[lc]=t;
179spm[lc]=m;
180}
181 */
182 }
183 }
184 } /*end for*/
185 if (lc == -1) {
186 minsp = -1;
187 } else {
188 mint = st[0];
189 minsp = spm[0];
190 for (m = 1; m <= lc; m++) {
191 if (st[m] < mint) {
192 mint = st[m];
193 minsp = spm[m];
194 }
195 }
196 for (j = 0; j < 3; j++) out_point[j] = in_point[j] + mint * in_vv[j];
197 /* transform_face_node(minsp,node);*/
198 /*map the out_point to ijk*/
199 for (i = 0; i < 3; i++) first_p[i] = out_point[i];
200 }
201 /* for(i=0;i<3;i++)
202 ijk[i]=(int)((first_p[i]-vd->xyz0[i])/vd->dxyz[i]);
203 */
204 return (minsp);
205}
206#endif
207
208#ifdef Octree
209void search2_leave(Tree_pointer *in_voxel, double out_point[3],
210 Tree_pointer_ptr *next_voxel, double ray_direction[3]) {
211 int i, index[3], child_no;
212 Tree_pointer *p1;
213 double mid_value;
214 p1 = in_voxel;
215 while (p1->child != NULL) {
216 for (i = 0; i < 3; i++) {
217 mid_value = (p1->bound_box[i * 2 + 1] - p1->bound_box[i * 2]) / 2.0 +
218 p1->bound_box[i * 2];
219 if (ray_direction[i] < 0) {
220 if (out_point[i] <= mid_value)
221 index[i] = 0;
222 else
223 index[i] = 1;
224 } else if (ray_direction[i] >= 0) {
225 if (out_point[i] < mid_value)
226 index[i] = 0;
227 else
228 index[i] = 1;
229 }
230 }
231 child_no = index[2] * 4 + index[1] * 2 + index[0];
232 p1 = &(p1->child[child_no]);
233 }
234 *next_voxel = p1;
235 return;
236}
237
238void search_leave(Tree_pointer *root_tree, int ini_index, double first_p[3],
239 Tree_pointer_ptr *voxel_p, double ray_direction[3]) {
240 int i;
241 int index[3], child_no;
242 Tree_pointer *p1;
243 double mid_value;
244
245 p1 = &(root_tree[ini_index]);
246 while (p1->child != NULL) {
247 for (i = 0; i < 3; i++) {
248 mid_value = (p1->bound_box[i * 2 + 1] - p1->bound_box[i * 2]) / 2.0 +
249 p1->bound_box[i * 2];
250 if (ray_direction[i] < 0) {
251 if (first_p[i] <= mid_value)
252 index[i] = 0;
253 else
254 index[i] = 1;
255 } else if (ray_direction[i] >= 0) {
256 if (first_p[i] < mid_value)
257 index[i] = 0;
258 else
259 index[i] = 1;
260 }
261 }
262 child_no = index[2] * 4 + index[1] * 2 + index[0];
263 p1 = &(p1->child[child_no]);
264 }
265 *voxel_p = p1;
266
267 return;
268}
269#endif
270
271int find_first_inter(double point_o[3], double view_point_d[3], int r_level[3],
272 double orig_xyz[3], double dxyz[3], double r_dxyz[3],
273 double ray_direction[3], double first_p[3], int ijk[3]) {
274 int i;
275 int ini_face, intersection, con1, next_ijk[3], face_sect[3];
276 double out_point[3];
277
278 ini_face = find_first_inter_point(point_o, view_point_d, ray_direction,
279 orig_xyz, dxyz, r_level, first_p);
280 if (ini_face == -1) { /*this ray has no intersection with dataset */
281 intersection = 0;
282 } else {
283 intersection = 1;
284 for (i = 0; i < 3; i++)
285 ijk[i] = (int)((first_p[i] - orig_xyz[i] - EPSILON) / r_dxyz[i]);
286 if ((ijk[0] < 0) || (ijk[0] > r_level[0] - 1) || (ijk[1] < 0) ||
287 (ijk[1] > r_level[1] - 1) || (ijk[2] < 0) || (ijk[2] > r_level[2] - 1))
289 "ERROR: HEC-MW-VIS-E2006:There is something wrong on finding the "
290 "first point");
291 for (i = 0; i < 3; i++) out_point[i] = 0.0;
292 con1 = find_out_point(orig_xyz, r_dxyz, ijk, first_p, ray_direction,
293 face_sect, out_point);
294 if (con1 == 0) {
295 /* start voxel maybe at the neighbour voxel*/
296 find_next_cell(r_level, ijk, face_sect, ray_direction, next_ijk);
297 if ((next_ijk[0] < 0) || (next_ijk[0] > r_level[0] - 1) ||
298 (next_ijk[1] < 0) || (next_ijk[1] > r_level[1] - 1) ||
299 (next_ijk[2] < 0) || (next_ijk[2] > r_level[2] - 1))
300 intersection = 0;
301
302 for (i = 0; i < 3; i++) ijk[i] = next_ijk[i];
303 }
304 }
305 return intersection;
306}
307
308#ifdef Octree
309int find_first_inter(double point_o[3], Parameter_vr *vr, VR_data *vd,
310 Tree_pointer *root_tree, Tree_pointer_ptr *voxel_p,
311 double ray_direction[3], double first_p[3])
312
313{
314 int i, ijk[3], ini_index, ijkn[3];
315 int ini_face, intersection, con1, face_sect[3], next_ijk[3];
316 double view_p[3], out_point[3];
318
319 for (i = 0; i < 3; i++) ijkn[i] = vd->nxyz[i];
320
321 for (i = 0; i < 3; i++) view_p[i] = vr->view_point_d[i];
322 ini_face = find_first_inter_point(point_o, view_p, vd, first_p);
323 if (ini_face == -1) { /*this ray has no intersection with dataset */
324 intersection = 0;
325 } else {
326 intersection = 1;
327 /* for(i=0;i<3;i++)
328ijk[i]=(int)((first_p[i]-vd->xyz0[i]-EPSILON)/vd->dxyz[i]);
329 */
330 if (ini_face == 0) {
331 ijk[0] = 0;
332 for (i = 1; i < 3; i++)
333 ijk[i] = (int)((first_p[i] - vd->xyz0[i]) / vd->dxyz[i]);
334 } else if (ini_face == 1) {
335 ijk[0] = vd->nxyz[0] - 1;
336 for (i = 1; i < 3; i++)
337 ijk[i] = (int)((first_p[i] - vd->xyz0[i]) / vd->dxyz[i]);
338 } else if (ini_face == 2) {
339 ijk[1] = 0;
340 ijk[0] = (int)((first_p[0] - vd->xyz0[0]) / vd->dxyz[0]);
341 ijk[2] = (int)((first_p[2] - vd->xyz0[2]) / vd->dxyz[2]);
342 } else if (ini_face == 3) {
343 ijk[1] = vd->nxyz[1] - 1;
344 ijk[0] = (int)((first_p[0] - vd->xyz0[0]) / vd->dxyz[0]);
345 ijk[2] = (int)((first_p[2] - vd->xyz0[2]) / vd->dxyz[2]);
346 } else if (ini_face == 4) {
347 ijk[2] = 0;
348 ijk[0] = (int)((first_p[0] - vd->xyz0[0]) / vd->dxyz[0]);
349 ijk[1] = (int)((first_p[1] - vd->xyz0[1]) / vd->dxyz[1]);
350 } else if (ini_face == 5) {
351 ijk[2] = vd->nxyz[2] - 1;
352 ijk[0] = (int)((first_p[0] - vd->xyz0[0]) / vd->dxyz[0]);
353 ijk[1] = (int)((first_p[1] - vd->xyz0[1]) / vd->dxyz[1]);
354 }
355
356 for (i = 0; i < 3; i++) out_point[i] = 0.0;
357 con1 = find_out_point(vd, first_p, ray_direction, face_sect, out_point);
358 if (con1 == 0) {
359 /* start voxel maybe at the neighbour voxel*/
360 find_next_cell(ijkn, ijk, face_sect, ray_direction, next_ijk);
361 for (i = 0; i < 3; i++) ijk[i] = next_ijk[i];
362 ini_index =
363 ijk[2] * vd->nxyz[1] * vd->nxyz[0] + ijk[1] * vd->nxyz[0] + ijk[0];
364 }
365 search_leave(root_tree, ini_index, first_p, &vvv, ray_direction);
366 vvv->local_face_in = ini_face;
367 *voxel_p = vvv;
368 }
369 return (intersection);
370}
371
372void build_adjacent_index(int connect[8][6], int local_face[8][6]) {
373 connect[0][0] = connect[0][2] = connect[0][4] = -1;
374 local_face[0][0] = 0;
375 local_face[0][2] = 2;
376 local_face[0][4] = 4;
377 connect[0][1] = 1;
378 local_face[0][1] = 0;
379 connect[0][3] = 2;
380 local_face[0][3] = 2;
381 connect[0][5] = 4;
382 local_face[0][5] = 4;
383 connect[1][1] = connect[1][2] = connect[1][4] = -1;
384 local_face[1][1] = 1;
385 local_face[1][2] = 2;
386 local_face[1][4] = 4;
387 connect[1][0] = 0;
388 local_face[1][0] = 1;
389 connect[1][3] = 3;
390 local_face[1][3] = 2;
391 connect[1][5] = 5;
392 local_face[1][5] = 4;
393 connect[2][0] = connect[2][3] = connect[2][4] = -1;
394 local_face[2][0] = 0;
395 local_face[2][3] = 3;
396 local_face[2][4] = 4;
397 connect[2][1] = 3;
398 local_face[2][1] = 0;
399 connect[2][2] = 0;
400 local_face[2][2] = 3;
401 connect[2][5] = 6;
402 local_face[2][5] = 4;
403 connect[3][1] = connect[3][3] = connect[3][4] = -1;
404 local_face[3][1] = 1;
405 local_face[3][3] = 3;
406 local_face[3][4] = 4;
407 connect[3][0] = 2;
408 local_face[3][0] = 1;
409 connect[3][2] = 1;
410 local_face[3][2] = 3;
411 connect[3][5] = 7;
412 local_face[3][5] = 4;
413 connect[4][0] = connect[4][2] = connect[4][5] = -1;
414 local_face[4][0] = 0;
415 local_face[4][2] = 2;
416 local_face[4][5] = 5;
417 connect[4][1] = 5;
418 local_face[4][1] = 0;
419 connect[4][3] = 6;
420 local_face[4][3] = 2;
421 connect[4][4] = 0;
422 local_face[4][4] = 5;
423 connect[5][1] = connect[5][2] = connect[5][5] = -1;
424 local_face[5][1] = 1;
425 local_face[5][2] = 2;
426 local_face[5][5] = 5;
427 connect[5][0] = 4;
428 local_face[5][0] = 1;
429 connect[5][3] = 7;
430 local_face[5][3] = 2;
431 connect[5][4] = 1;
432 local_face[5][4] = 5;
433 connect[6][0] = connect[6][3] = connect[6][5] = -1;
434 local_face[6][0] = 0;
435 local_face[6][3] = 3;
436 local_face[6][5] = 5;
437 connect[6][1] = 7;
438 local_face[6][1] = 0;
439 connect[6][2] = 4;
440 local_face[6][2] = 3;
441 connect[6][4] = 2;
442 local_face[6][4] = 5;
443 connect[7][1] = connect[7][3] = connect[7][5] = -1;
444 local_face[7][1] = 1;
445 local_face[7][3] = 3;
446 local_face[7][5] = 5;
447 connect[7][0] = 6;
448 local_face[7][0] = 1;
449 connect[7][2] = 5;
450 local_face[7][2] = 3;
451 connect[7][4] = 3;
452 local_face[7][4] = 5;
453 return;
454}
455#endif
456
457#ifdef old_intersection
458static void find_coordinates_of_cell(int r_level[3], double orig_xyz[3],
459 double r_dxyz[3], int ijk[3],
460 double vv[8 * 3]) {
461 vv[0 * 3] = vv[4 * 3] = vv[7 * 3] = vv[3 * 3] =
462 orig_xyz[0] + ijk[0] * r_dxyz[0];
463 vv[1 * 3] = vv[5 * 3] = vv[6 * 3] = vv[2 * 3] =
464 orig_xyz[0] + (ijk[0] + 1) * r_dxyz[0];
465 vv[0 * 3 + 2] = vv[1 * 3 + 2] = vv[2 * 3 + 2] = vv[3 * 3 + 2] =
466 orig_xyz[2] + ijk[2] * r_dxyz[2];
467 vv[4 * 3 + 2] = vv[5 * 3 + 2] = vv[6 * 3 + 2] = vv[7 * 3 + 2] =
468 orig_xyz[2] + (ijk[2] + 1) * r_dxyz[2];
469 vv[0 * 3 + 1] = vv[1 * 3 + 1] = vv[5 * 3 + 1] = vv[4 * 3 + 1] =
470 orig_xyz[1] + ijk[1] * r_dxyz[1];
471 vv[3 * 3 + 1] = vv[2 * 3 + 1] = vv[6 * 3 + 1] = vv[7 * 3 + 1] =
472 orig_xyz[1] + (ijk[1] + 1) * r_dxyz[1];
473 return;
474}
475
476int find_out_point(VR_data *vd, int ijk[3], double in_point[3],
477 double ray_direction[3], int face_sect[3],
478 double out_point[3])
479/*Calculate the coordiates of exit point*/
480{
481 double fp[4][3], t, a[6], b[6], c[6], d[6], abc, st[6], mint, vv[24];
482 int i, j, m, lc, spm[6], con1, minsp, node[4], face1_sect[3], lm;
483
484 find_coordinates_of_cell(vd, ijk, vv);
485 lc = -1;
486 con1 = 1;
487 for (m = 0; m < 6; m++) {
488 transform_face_node(m, node);
489 for (i = 0; i < 4; i++)
490 for (j = 0; j < 3; j++) {
491 fp[i][j] = vv[node[i] * 3 + j];
492 /* fv[i][j]=cell->v_data[node[i]*3+j];
493 */
494 }
495 a[m] = (fp[1][1] - fp[0][1]) * (fp[3][2] - fp[0][2]) -
496 (fp[3][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
497 b[m] = -(fp[1][0] - fp[0][0]) * (fp[3][2] - fp[0][2]) +
498 (fp[3][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
499 c[m] = (fp[1][0] - fp[0][0]) * (fp[3][1] - fp[0][1]) -
500 (fp[3][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
501 abc = sqrt(a[m] * a[m] + b[m] * b[m] + c[m] * c[m]);
502 a[m] /= abc;
503 b[m] /= abc;
504 c[m] /= abc;
505 d[m] = -(fp[0][0] * a[m] + fp[0][1] * b[m] + fp[0][2] * c[m]);
506 }
507 lm = -1;
508 for (m = 0; m < 3; m++) face_sect[m] = -1;
509 for (m = 0; m < 6; m++) {
510 if (fabs(a[m] * in_point[0] + b[m] * in_point[1] + c[m] * in_point[2] +
511 d[m]) < 0.0000001) {
512 lm++;
513 face_sect[lm] = m;
514 }
515 }
516 for (m = 0; m < 6; m++) {
517 if ((m != face_sect[0]) && (m != face_sect[1]) && (m != face_sect[2])) {
518 if (fabs(a[m] * ray_direction[0] + b[m] * ray_direction[1] +
519 c[m] * ray_direction[2]) > 0.000001) {
520 t = (-d[m] - a[m] * in_point[0] - b[m] * in_point[1] -
521 c[m] * in_point[2]) /
522 (a[m] * ray_direction[0] + b[m] * ray_direction[1] +
523 c[m] * ray_direction[2]);
524 /* printf("t=%f\n",t);*/
525 if (t > 0) {
526 lc++;
527 st[lc] = t;
528 spm[lc] = m;
529 }
530 }
531 }
532 } /*end for*/
533 if (lc == -1)
534 con1 = 0;
535 else {
536 mint = st[0];
537 minsp = spm[0];
538 for (m = 1; m <= lc; m++) {
539 if (st[m] < mint) {
540 mint = st[m];
541 minsp = spm[m];
542 }
543 }
544 /* *delta_t=mint;
545 */
546 for (j = 0; j < 3; j++)
547 out_point[j] = in_point[j] + mint * ray_direction[j];
548 /* transform_face_node(minsp,node);
549get_velocity(node,cell,out_point,out_vv);
550 */
551 /* con1=judge_on_face(out_point, node, cell);
552if(con1==0) {
553 printf("*****\n");
554
555}*/
556
557 lm = -1;
558 for (m = 0; m < 3; m++) face1_sect[m] = -1;
559 for (m = 0; m < 6; m++) {
560 if (fabs(a[m] * out_point[0] + b[m] * out_point[1] + c[m] * out_point[2] +
561 d[m]) < 0.000001) {
562 lm++;
563 face1_sect[lm] = m;
564 }
565 }
566
567 for (m = 0; m < 3; m++) face_sect[m] = face1_sect[m];
568 }
569 /* in_voxel.local_face_out=face_sect[0];
570 */
571
572 /* printf("con1=%d\n",con1);*/
573 return (con1);
574}
575
576#endif
577
578static int find_out_point(double orig_xyz[3], double r_dxyz[3], int ijk[3],
579 double in_point[3], double ray_direction[3],
580 int face_sect[3], double out_point[3]) {
581 int i, j, mincomp, intersection;
582 double minmax[6], t[3], mint;
583 int i1, j1;
584 for (i = 0; i < 3; i++) face_sect[i] = -1;
585 for (i = 0; i < 6; i++) {
586 i1 = i / 2;
587 j1 = i % 2;
588 minmax[i] = orig_xyz[i1] + (ijk[i1] + j1) * r_dxyz[i1];
589 }
590#ifdef old_version
591 for (i = 0; i < 3; i++) {
592 if (fabs(ray_direction[i]) < EPSILON)
593 t[i] = 1.0E+17;
594 else {
595 if (ray_direction[i] > 0)
596 t[i] = (minmax[i * 2 + 1] - in_point[i]) / ray_direction[i];
597 else if (ray_direction[i] < 0)
598 t[i] = (in_point[i] - minmax[i * 2]) / (-ray_direction[i]);
599 }
600 }
601#endif
602 for (i = 0; i < 3; i++) {
603 if (ray_direction[i] > 0)
604 t[i] = (minmax[i * 2 + 1] - in_point[i]) / (ray_direction[i] + EPSILON);
605 else if (ray_direction[i] < 0)
606 t[i] = (in_point[i] - minmax[i * 2]) / (-ray_direction[i] + EPSILON);
607 }
608
609 mint = 1.0E+17;
610 mincomp = -1;
611 for (i = 0; i < 3; i++) {
612 if ((mint > t[i]) && (t[i] > EPSILON)) {
613 for (j = 0; j < 3; j++)
614 out_point[j] = in_point[j] + t[i] * ray_direction[j];
615 if ((out_point[0] >= minmax[0] - EPSILON) &&
616 (out_point[0] <= minmax[1] + EPSILON) &&
617 (out_point[1] >= minmax[2] - EPSILON) &&
618 (out_point[1] <= minmax[3] + EPSILON) &&
619 (out_point[2] >= minmax[4] - EPSILON) &&
620 (out_point[2] <= minmax[5] + EPSILON)) {
621 mint = t[i];
622 mincomp = i;
623 }
624 }
625 }
626 if ((mincomp >= 0) && (mincomp <= 2)) {
627 intersection = 1;
628 for (i = 0; i < 3; i++)
629 out_point[i] = in_point[i] + mint * ray_direction[i];
630 } else {
631 intersection = 0;
632 for (i = 0; i < 3; i++) out_point[i] = in_point[i];
633 }
634 for (i = 0; i < 3; i++) {
635 if (fabs(out_point[i] - minmax[i * 2]) < EPSILON)
636 face_sect[i] = i * 2;
637 else if (fabs(out_point[i] - minmax[i * 2 + 1]) < EPSILON)
638 face_sect[i] = i * 2 + 1;
639 }
640
641 return (intersection);
642}
643
644static void find_next_cell(int ijkn[3], int current_ijk[3], int face_sect[3],
645 double vv[3], int next_ijk[3])
646/*get the code of next cell*/
647{
648 int j, m, flag;
649 int face1_sect[3];
650
651 flag = 1;
652 for (j = 0; j < 3; j++) face1_sect[j] = face_sect[j];
653 for (j = 0; j < 3; j++) next_ijk[j] = current_ijk[j];
654 for (m = 0; m < 3; m++) {
655 if (face_sect[m] != -1) {
656 if (face_sect[m] == 0) {
657 if (vv[0] < 0) {
658 next_ijk[0] = current_ijk[0] - 1;
659 face1_sect[m] = 1;
660 if (next_ijk[0] < 0) break;
661 }
662 } else if (face_sect[m] == 1) {
663 if (vv[0] > 0) {
664 next_ijk[0] = current_ijk[0] + 1;
665 face1_sect[m] = 0;
666 if (next_ijk[0] > ijkn[0] - 1) break;
667 }
668 } else if (face_sect[m] == 2) {
669 if (vv[1] < 0) {
670 next_ijk[1] = current_ijk[1] - 1;
671 face1_sect[m] = 3;
672 if (next_ijk[1] < 0) break;
673 }
674 } else if (face_sect[m] == 3) {
675 if (vv[1] > 0) {
676 next_ijk[1] = current_ijk[1] + 1;
677 face1_sect[m] = 2;
678 if (next_ijk[1] > ijkn[1] - 1) break;
679 }
680 } else if (face_sect[m] == 4) {
681 if (vv[2] < 0) {
682 next_ijk[2] = current_ijk[2] - 1;
683 face1_sect[m] = 5;
684 if (next_ijk[2] < 0) break;
685 }
686 } else if (face_sect[m] == 5) {
687 if (vv[2] > 0) {
688 next_ijk[2] = current_ijk[2] + 1;
689 face1_sect[m] = 4;
690 if (next_ijk[2] > ijkn[2] - 1) break;
691 }
692 }
693 }
694 }
695
696 for (m = 0; m < 3; m++) face_sect[m] = face1_sect[m];
697
698 return;
699}
700
701static int find_next_point(double orig_xyz[3], double r_dxyz[3], int r_level[3],
702 int current_ijk[3], double in_point[3],
703 double out_point[3], double ray_direction[3],
704 int next_ijk[3]) {
705 int flag, face_sect[3];
706
707 flag = 1;
708 flag = find_out_point(orig_xyz, r_dxyz, current_ijk, in_point, ray_direction,
709 face_sect, out_point);
710 /* con1=find_out_point(vd, current_ijk, in_point,ray_direction,face_sect,out_point);
711 */ if (flag == 0)
713 "ERROR: HEC-MW-VIS-E2007:There is some problem in finding intersection "
714 "point");
715
716 /* for(i=0;i<3;i++)
717 mid_point[i]=(out_point[i]+in_point[i])/2.0;
718 for(i=0;i<3;i++)
719 current_ijk[i]=(int)((mid_point[i]-vd->xyz0[i])/vd->dxyz[i]);
720 */
721 if (flag == 1) {
722 find_next_cell(r_level, current_ijk, face_sect, ray_direction, next_ijk);
723 if ((next_ijk[0] < 0) || (next_ijk[0] > r_level[0] - 1) ||
724 (next_ijk[1] < 0) || (next_ijk[1] > r_level[1] - 1) ||
725 (next_ijk[2] < 0) || (next_ijk[2] > r_level[2] - 1))
726 flag = -1;
727 }
728 return (flag);
729}
730
731void ray_trace(int remove_0_display_on, int color_mapping_style,
732 double *interval_point, int transfer_function_style,
733 double opa_value, int num_of_features, double *fea_point,
734 double view_point_d[3], int interval_mapping_num,
735 int color_system_type, int num_of_lights, double *light_point,
736 double k_ads[3], double orig_xyz[3], double dxyz[3],
737 double r_dxyz[3], int r_level[3], int *empty_flag, double *var,
738 double *grad_var, double first_p[3], int first_ijk[3],
739 double ray_direction[3], double mincolor, double maxcolor,
740 double accum_rgba[4], double grad_minmax[2],
741 double feap_minmax[2], double feai_minmax[2],
742 double dis_minmax[2], double *opa_table, double tav_length,
743 int time_step, int test_i, int test_j)
744
745/*void ray_trace(Parameter_vr *vr, VR_data *vd, Tree_pointer *root_tree,
746 Tree_pointer *voxel_p,
747 double point_o[3], double ray_direction[3], double
748 mincolor, double maxcolor,
749 int connect[8][6], int local_face[8][6], double
750 accum_rgba[4], double grad_minmax[2],
751 double feap_minmax[2], double feai_minmax[2],
752 double dis_minmax[2], double *opa_table, double
753 tav_length, int time_step)
754 */
755{
756 int i, j;
757 double in_point[3], out_point[3];
758 int flag, single_empty_flag, value_flag, index, current_ijk[3], next_ijk[3];
759 int print_flag;
760
761 for (i = 0; i < 4; i++) accum_rgba[i] = 0.0;
762 flag = 1;
763 for (i = 0; i < 3; i++) in_point[i] = first_p[i];
764 for (i = 0; i < 3; i++) current_ijk[i] = first_ijk[i];
765 while ((accum_rgba[3] < 0.99) && (flag == 1)) {
766 /*
767while((accum_rgba[3]<0.99) && (flag==1) && (accum_rgba[0]<0.99) &&
768(accum_rgba[1]<0.99)
769&& (accum_rgba[2]<0.99)) {
770 */
771 flag = find_next_point(orig_xyz, r_dxyz, r_level, current_ijk, in_point,
772 out_point, ray_direction, next_ijk);
773
774 single_empty_flag = 1;
775 value_flag = 1;
776 /* for(k=0;k<2;k++)
777 for(j=0;j<2;j++)
778 for(i=0;i<2;i++) {
779
780 index=(current_ijk[2]+k)*(r_level[0]+1)*(r_level[1]+1)+(current_ijk[1]+j)*(r_level[0]+1)+
781 current_ijk[0]+i;
782 */
783 for (i = 0; i < 8; i++) {
784 j = (int)i % 4;
785 index = (current_ijk[2] + (int)(i / 4)) * (r_level[0] + 1) *
786 (r_level[1] + 1) +
787 (current_ijk[1] + (int)(j / 2)) * (r_level[0] + 1) +
788 current_ijk[0] + j % 2;
789 if (empty_flag[index] == 0) single_empty_flag = 0;
790 /* if(fabs(var[index])>EPSILON)
791 value_flag=1;
792 */
793 }
794 /* value_flag=1;
795if(remove_0_display_on==1) {
796 value_flag=0;
797for(i=0;i<8;i++) {
798j
799 if(fabs(var[index])>EPSILON) {
800 value_flag=1;
801 break;
802 }
803 }
804}
805
806overlap_flag=1;
807for(k=0;k<3;k++) {
808 if(fabs(in_point[k]-(vd->xyz0[k]+vd->dxyz[k]))<EPSILON) {
809 overlap_flag=0;
810 break;
811 }
812}
813 */
814
815 /*
816if((empty_flag==1) && (value_flag==1))
817 */
818 if ((single_empty_flag == 1) && (flag != 0) && (value_flag == 1))
819 compute_color_vr(current_ijk, color_mapping_style, interval_point,
820 transfer_function_style, opa_value, num_of_features,
821 fea_point, view_point_d, interval_mapping_num,
822 color_system_type, num_of_lights, light_point, k_ads,
823 r_level, orig_xyz, r_dxyz, var, grad_var, accum_rgba,
824 mincolor, maxcolor, grad_minmax, feap_minmax,
825 feai_minmax, dis_minmax, opa_table, in_point, out_point,
826 tav_length, time_step, print_flag);
827 if (flag == 1) {
828 for (i = 0; i < 3; i++) in_point[i] = out_point[i];
829 for (i = 0; i < 3; i++) current_ijk[i] = next_ijk[i];
830 }
831 /* if((vd->empty_flag[current_voxel->cell_id-vd->nxyz[0]*vd->nxyz[1]*vd->nxyz[2]]==1)
832&& ((fabs(vd->var[current_voxel->cell_id])>EPSILON) ||
833 (fabs(vd->grad_var[current_voxel->cell_id*3])>EPSILON) ||
834 (fabs(vd->grad_var[current_voxel->cell_id*3+1])>EPSILON) ||
835 (fabs(vd->grad_var[current_voxel->cell_id*3+2])>EPSILON)))
836 */
837 /* if((vd->empty_flag[current_voxel->cell_id-vd->nxyz[0]*vd->nxyz[1]*vd->nxyz[2]]==1)
838 && (fabs(vd->var[current_voxel->cell_id])>EPSILON))
839 */
840 }
841 return;
842}
if(!(yy_init))
Definition: hecmw_ablex.c:1305
#define NULL
void compute_color_vr(int current_ijk[3], int color_mapping_style, double *interval_point, int transfer_function_style, double opa_value, int num_of_features, double *fea_point, double view_point_d[3], int interval_mapping_num, int color_system_type, int num_of_lights, double *light_point, double k_ads[3], int r_level[3], double orig_xyz[3], double r_dxyz[3], double *var, double *grad_var, double accum_rgba[4], double mincolor, double maxcolor, double grad_minmax[2], double feap_minmax[2], double feai_minmax[2], double dis_minmax[2], double *opa_table, double in_point[3], double out_point[3], double tav_length, int time_step, int print_flag)
#define EPSILON
void HECMW_vis_print_exit(char *var)
void transform_face_node(int face, int node[4])
#define SQR(x)
int find_first_inter(double point_o[3], double view_point_d[3], int r_level[3], double orig_xyz[3], double dxyz[3], double r_dxyz[3], double ray_direction[3], double first_p[3], int ijk[3])
void ray_trace(int remove_0_display_on, int color_mapping_style, double *interval_point, int transfer_function_style, double opa_value, int num_of_features, double *fea_point, double view_point_d[3], int interval_mapping_num, int color_system_type, int num_of_lights, double *light_point, double k_ads[3], double orig_xyz[3], double dxyz[3], double r_dxyz[3], int r_level[3], int *empty_flag, double *var, double *grad_var, double first_p[3], int first_ijk[3], double ray_direction[3], double mincolor, double maxcolor, double accum_rgba[4], double grad_minmax[2], double feap_minmax[2], double feai_minmax[2], double dis_minmax[2], double *opa_table, double tav_length, int time_step, int test_i, int test_j)
struct _tree_pointer_struct * child