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