FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 
6 #include "hecmw_vis_ray_trace.h"
7 
8 #include <math.h>
9 #include "hecmw_vis_mem_util.h"
10 #include "hecmw_vis_color_composite_vr.h"
11 
12 #if defined(old_version) || defined(old_intersection)
13 #include "hecmw_vis_new_refine.h"
14 #endif
15 
16 /*
17 void transform_face_node(int ijkn[3], int current_ijk[3], int face_sect[3],
18 double vv[3],
19  int next_ijk[3]);
20  */
21 static 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]);
24 static void find_next_cell(int ijkn[3], int current_ijk[3], int face_sect[3],
25  double vv[3], int next_ijk[3]);
26 
27 static 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 
90 int 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;
155 if((m==0) || (m==1)) {
156 if((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))
160 on_face_flag=1;
161 }
162 else if((m==2) || (m==3)) {
163 if((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 }
169 else if((m==4) || (m==5)) {
170 if((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))
174 on_face_flag=1;
175 }
176 if(on_face_flag==1) {
177 lc++;
178 st[lc]=t;
179 spm[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
209 void 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 
238 void 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 
271 int 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))
288  HECMW_vis_print_exit(
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
309 int 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];
317  Tree_pointer_ptr vvv;
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++)
328 ijk[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 
372 void 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
458 static 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 
476 int 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 coordinates 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);
549 get_velocity(node,cell,out_point,out_vv);
550  */
551  /* con1=judge_on_face(out_point, node, cell);
552 if(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 
578 static 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 
644 static 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 
701 static 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)
712  HECMW_vis_print_exit(
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 
731 void 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  /*
767 while((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;
795 if(remove_0_display_on==1) {
796  value_flag=0;
797 for(i=0;i<8;i++) {
798 j
799  if(fabs(var[index])>EPSILON) {
800  value_flag=1;
801  break;
802  }
803  }
804 }
805 
806 overlap_flag=1;
807 for(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  /*
816 if((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 }
_data_structured_vr_struct::xyz0
double xyz0[3]
Definition: hecmw_vis_ray_trace.h:144
_data_structured_vr_struct::nxyz
int nxyz[3]
Definition: hecmw_vis_ray_trace.h:143
_tree_pointer_struct
Definition: hecmw_vis_ray_trace.h:158
if
if(!(yy_init))
Definition: hecmw_ablex.c:1823
ray_trace
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)
Definition: hecmw_vis_ray_trace.c:731
hecmw_vis_mem_util.h
hecmw_vis_ray_trace.h
_tree_pointer_struct::bound_box
double bound_box[6]
Definition: hecmw_vis_ray_trace.h:162
_tree_pointer_struct::local_face_in
int local_face_in
Definition: hecmw_vis_ray_trace.h:164
_vr_parameter_struct
Definition: hecmw_vis_ray_trace.h:56
HECMW_vis_print_exit
void HECMW_vis_print_exit(char *var)
Definition: hecmw_vis_mem_util.c:21
hecmw_vis_color_composite_vr.h
hecmw_vis_new_refine.h
compute_color_vr
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)
Definition: hecmw_vis_color_composite_vr.c:405
find_first_inter
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])
Definition: hecmw_vis_ray_trace.c:271
_data_structured_vr_struct
Definition: hecmw_vis_ray_trace.h:138
_data_structured_vr_struct::dxyz
double dxyz[3]
Definition: hecmw_vis_ray_trace.h:139
_vr_parameter_struct::view_point_d
double view_point_d[3]
Definition: hecmw_vis_ray_trace.h:64
EPSILON
#define EPSILON
Definition: hecmw_vis_color_mapping.c:10
NULL
#define NULL
Definition: hecmw_io_nastran.c:30
SQR
#define SQR(x)
Definition: hecmw_vis_psf_rendering.h:15
transform_face_node
void transform_face_node(int face, int node[4])
Definition: hecmw_vis_new_refine.c:19
_tree_pointer_struct::child
struct _tree_pointer_struct * child
Definition: hecmw_vis_ray_trace.h:166