FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_vis_color_composite_vr.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>
9 #include "hecmw_vis_ray_trace.h"
10 #include "hecmw_malloc.h"
11 
12 #if 0
13 void find_sample_point(double in_point[3], double out_point[3], int num_sample, double *sample_point)
14 {
15  int i,j;
16  double t;
17 
18  t=1.0/(double)(num_sample+1);
19  for(i=0;i<num_sample;i++) {
20  for(j=0;j<3;j++)
21  sample_point[i*3+j]=in_point[j]+t*(i+1)*(out_point[j]-in_point[j]);
22  }
23  return;
24 }
25 
26 double value_compute(double point[3], int current_ijk[3], double orig_xyz[3], double r_dxyz[3], int r_level[3], double *var)
27 {
28  int i, j,k, index, m;
29  double value[8], vv[8*3], v;
30  double dis[8], d;
31  int return_flag;
32  return_flag=0;
33  vv[0*3]=vv[4*3]=vv[7*3]=vv[3*3]=orig_xyz[0]+current_ijk[0]*r_dxyz[0];
34  vv[1*3]=vv[5*3]=vv[6*3]=vv[2*3]=orig_xyz[0]+(current_ijk[0]+1)*r_dxyz[0];
35  vv[0*3+2]=vv[1*3+2]=vv[2*3+2]=vv[3*3+2]=orig_xyz[2]+current_ijk[2]*r_dxyz[2];
36  vv[4*3+2]=vv[5*3+2]=vv[6*3+2]=vv[7*3+2]=orig_xyz[2]+(current_ijk[2]+1)*r_dxyz[2];
37  vv[0*3+1]=vv[1*3+1]=vv[5*3+1]=vv[4*3+1]=orig_xyz[1]+current_ijk[1]*r_dxyz[1];
38  vv[3*3+1]=vv[2*3+1]=vv[6*3+1]=vv[7*3+1]=orig_xyz[1]+(current_ijk[1]+1)*r_dxyz[1];
39 
40 
41  for(m=0;m<8;m++) {
42  i=m % 2;
43  j=(m/2) % 2;
44  k=(m/2) /2;
45  index=(current_ijk[2]+k)*(r_level[0]+1)*(r_level[1]+1)+(current_ijk[1]+j)*(r_level[0]+1)+
46  current_ijk[0]+i;
47  value[k*4+j*2+i]=var[index];
48  }
49 
50  for(i=0;i<8;i++) {
51  dis[i]=sqrt(SQR(point[0]-vv[i*3])+SQR(point[1]-vv[i*3+1])+SQR(point[2]-vv[i*3+2]));
52  if(fabs(dis[i])<EPSILON) {
53  v=value[i];
54  return_flag=1;
55  }
56  }
57  if(return_flag==0) {
58  d=0.0;
59  for(i=0;i<8;i++)
60  d+=1.0/dis[i];
61  v=0.0;
62  for(i=0;i<8;i++) {
63  v+=value[i]/(dis[i]*d);
64 
65  }
66  }
67  return v;
68 }
69 
70 
71 /*
72 void find_grad_minmax(int r_nxyz[3], grad_var, double grad_minmax[2])
73 {
74  int i;
75  double grad;
76  grad_minmax[0]=grad_minmax[1]=sqrt(SQR(grad_var[0])+SQR(grad_var[1])
77  +SQR(grad_var[2]));
78  for(i=0;i<(r_nxyz[0]+1)*(r_nxyz[1]+1)*(r_nxyz[2]+1);i++) {
79  grad=sqrt(SQR(grad_var[i*3])+SQR(grad_var[i*3+1])
80  +SQR(grad_var[i*3+2]));
81  if(grad<grad_minmax[0]) grad_minmax[0]=grad;
82  if(grad>grad_minmax[1]) grad_minmax[1]=grad;
83  }
84  return;
85 }
86  */
87 
88 
89 
90 
91 
92 double opacity_decision(int current_ijk[3], int transfer_function_style, double opa_value, int num_of_features, double *fea_point,
93  double view_point_d[3], int color_mapping_style, double *interval_point, int interval_mapping_num,
94  double orig_xyz[3], double r_dxyz[3],
95  double value, double n[3], double grad_minmax[2],
96  double feap_minmax[2], double feai_minmax[2],
97  double dis_minmax[2], double *opa_table, double mincolor, double maxcolor, int time_step)
98 {
99  double opacity, grad;
100  double t, mint;
101  int i, /* cell_id, */ level, min_type;
102  double cp[3], del_l, dist;
103  /* cell_id=in_voxel->cell_id;
104  */
105  if(transfer_function_style==1)
106  opacity=opa_value;
107  else if(transfer_function_style==2) {
108  grad=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2]));
109  opacity=(grad-grad_minmax[0])/(grad_minmax[1]-grad_minmax[0])/200.0+0.0002;
110  }
111  else if(transfer_function_style==3) {
112  /* value=vd->var[cell_id*vd->varnumtot+vr->color_comp];
113  */
114  mint=1.0E17;
115  for(i=0;i<num_of_features;i++) {
116  t=fabs(value- fea_point[i*3]);
117  if(t<mint) {
118  mint=t;
119  min_type=i;
120  }
121  }
122  /* if(mint<0.00001) {
123  t=feap_minmax[1]-mint+feap_minmax[0];
124  opacity=(t-feap_minmax[0])/(feap_minmax[1]-feap_minmax[0])*3.05+0.15;
125  }
126  */
127  if(mint<fea_point[min_type*3+1]) {
128  opacity=fea_point[min_type*3+2]*(fea_point[min_type*3+1]-mint)/fea_point[min_type*3+1]+opa_value;
129  }
130  else opacity=opa_value;
131  }
132  else if(transfer_function_style==4) {
133  /* value=vd->var[cell_id*vd->varnumtot+vr->color_comp];
134  */
135  /* mint=1.0E17;
136  for(i=0;i<vr->num_of_features;i++) {
137  if((value>=vr->fea_point[i*3]) && (value<=vr->fea_point[i*3+1]))
138  t=0.0;
139  else {
140  t1=fabs(value-vr->fea_point[i*2]);
141  t2=fabs(value-vr->fea_point[i*2+1]);
142  if(t1<t2) t=t1;
143  else t=t2;
144  }
145  if(t<mint) mint=t;
146  }
147  t=feai_minmax[1]-mint+feai_minmax[0];
148  opacity=(t-feai_minmax[0])/(feai_minmax[1]-feai_minmax[0])*0.4+0.1;
149  */
150  opacity=opa_value;
151  for(i=0;i<num_of_features;i++) {
152  if((value>=fea_point[i*3]) && (value<=fea_point[i*3+1])) {
153  opacity=fea_point[i*3+2];
154  break;
155  }
156  }
157 
158  }
159  else if(transfer_function_style==5) {
160  /* for(i=0;i<3;i++)
161  cp[i]=(in_voxel->bound_box[i*2+1]+in_voxel->bound_box[i*2])/2.0;
162  */
163  for(i=0;i<3;i++)
164  cp[i]=(current_ijk[i]+0.5)*r_dxyz[i]+orig_xyz[i];
165  dist=sqrt(SQR(cp[0]-view_point_d[0])+SQR(cp[1]-view_point_d[1])
166  +SQR(cp[2]-view_point_d[2]));
167  dist=dis_minmax[1]-dist+dis_minmax[0];
168  opacity=(dist-dis_minmax[0])/(dis_minmax[1]-dis_minmax[0])/200.0+0.0002;
169  }
170  else if(transfer_function_style==6) {
171  /* for(i=0;i<3;i++)
172  cp[i]=(in_voxel->bound_box[i*2+1]+in_voxel->bound_box[i*2])/2.0;
173  */
174  for(i=0;i<3;i++)
175 
176  cp[i]=(current_ijk[i]+0.5)*r_dxyz[i]+orig_xyz[i];
177 
178  dist=sqrt(SQR(cp[0]-view_point_d[0])+SQR(cp[1]-view_point_d[1])
179  +SQR(cp[2]-view_point_d[2]));
180  opacity=(dist-dis_minmax[0])/(dis_minmax[1]-dis_minmax[0])/200.0+0.0002;
181  }
182  else if(transfer_function_style==7) {
183 
184 
185  if(color_mapping_style==1) {
186  if(fabs(maxcolor-mincolor)>EPSILON)
187  value=(value-mincolor)/(maxcolor-mincolor);
188  }
189 
190  if(color_mapping_style==2) {
191  mincolor=interval_point[0];
192  maxcolor=interval_point[1];
193  if(fabs(maxcolor-mincolor)>EPSILON)
194  value=(value-mincolor)/(maxcolor-mincolor);
195  }
196  if((color_mapping_style==3) || (color_mapping_style==4)) {
197  for(i=1;i<interval_mapping_num+1;i++) {
198  if((value<=interval_point[i*2]) && (value>interval_point[(i-1)*2])) {
199  value=(value-interval_point[(i-1)*2])/(interval_point[i*2]-interval_point[(i-1)*2])*
200  (interval_point[i*2+1]-interval_point[(i-1)*2+1])+interval_point[(i-1)*2+1];
201  break;
202  }
203  }
204  }
205 
206  if(value>1.0) value=1.0;
207  if(value<0.0) value=0.0;
208  time_step+=0;
209  /* ss=(sqrt((double)time_step)-6.0)/60.0;
210  opacity=((value-0.6)/6.0*((double)time_step-40.0)/7.0)+0.23-ss;
211 
212  ss=((double)time_step)/((time_step+280));
213  */
214  opacity=value/200.0+0.0002;;
215  if(opacity<0) opacity=0.0;
216  }
217  else if(transfer_function_style==8) {
218  del_l=(maxcolor-mincolor)/255.0;
219  level=(int)(value-mincolor)/del_l;
220  if(level<0) level=0;
221  if(level>255) level=255;
222  opacity=opa_table[level];
223  }
224 
225  return(opacity);
226 }
227 #endif
228 
229 #ifdef surface_vr
230 int find_surface_point(int index, VR_data *vd, double in[3], double out[3],
231  double *surf_p) {
232  int i, j, k, m, n1, n2, flag, flag_inside;
233  double fp[3][3], n_f[4], v_f[3], d, point[3], s[3], sina, cosa, ss, t, color,
234  dis[3], dist, n_norm, v_norm, nv_sign;
235  int num_surp, num_surp1;
236  double *dd, dd_norm;
237  int *o_flag;
238  double *surf_p1;
239  double tmp_dd, tmp_p[7];
240 
241  num_surp = 0;
242  surf_p1 = (double *)HECMW_calloc(7 * vd->surface[index].num, sizeof(double));
243  if (surf_p1 == NULL) {
244  fprintf(stderr, "There is no enough memory for surf_p\n"), exit(0);
245  }
246 
247  for (m = 0; m < vd->surface[index].num; m++) {
248  /*find the equation of the patch first */
249 
250  for (i = 0; i < 3; i++)
251  for (j = 0; j < 3; j++)
252  fp[i][j] = vd->surface[index].surf_data[m * 12 + i * 3 + j];
253  n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
254  (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
255  n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
256  (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
257  n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
258  (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
259  n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
260  if (fabs(n_norm) > EPSILON) {
261  for (j = 0; j < 3; j++) n_f[j] /= n_norm;
262  }
263  for (i = 0; i < 3; i++) v_f[i] = in[i] - out[i];
264  v_norm = sqrt(v_f[0] * v_f[0] + v_f[1] * v_f[1] + v_f[2] * v_f[2]);
265  if (fabs(v_norm) > EPSILON) {
266  for (i = 0; i < 3; i++) v_f[i] /= v_norm;
267  }
268  nv_sign = n_f[0] * v_f[0] + n_f[1] * v_f[1] + n_f[2] * v_f[2];
269  if (nv_sign < 0.0) {
270  for (i = 0; i < 3; i++) n_f[i] = -n_f[i];
271  }
272  n_f[3] = -n_f[0] * fp[0][0] - n_f[1] * fp[0][1] - n_f[2] * fp[0][2];
273  /*find intersection point*/
274  if (fabs(n_f[0] * (out[0] - in[0]) + n_f[1] * (out[1] - in[1]) +
275  n_f[2] * (out[2] - in[2])) > EPSILON) {
276  t = (-n_f[3] - n_f[0] * in[0] - n_f[1] * in[1] - n_f[2] * in[2]) /
277  (n_f[0] * (out[0] - in[0]) + n_f[1] * (out[1] - in[1]) +
278  n_f[2] * (out[2] - in[2]));
279  if ((t > -EPSILON) && (t < 1.0)) {
280  for (j = 0; j < 3; j++) point[j] = in[j] + t * (out[j] - in[j]);
281  /*judge whether it is inside the range of patch by equal area method*/
282  flag_inside = 0;
283  for (j = 0; j < 3; j++) {
284  n1 = j;
285  n2 = j + 1;
286  if (n2 == 3) n2 = 0;
287  if ((fabs(point[0] - fp[n1][0]) < EPSILON) &&
288  (fabs(point[1] - fp[n1][1]) < EPSILON) &&
289  (fabs(point[2] - fp[n1][2]) < EPSILON)) {
290  flag_inside = 1;
291  break;
292  } else
293  cosa =
294  ((point[0] - fp[n1][0]) * (fp[n2][0] - fp[n1][0]) +
295  (point[1] - fp[n1][1]) * (fp[n2][1] - fp[n1][1]) +
296  (point[2] - fp[n1][2]) * (fp[n2][2] - fp[n1][2])) /
297  (sqrt(SQR(point[0] - fp[n1][0]) + SQR(point[1] - fp[n1][1]) +
298  SQR(point[2] - fp[n1][2])) *
299  sqrt(SQR(fp[n2][0] - fp[n1][0]) + SQR(fp[n2][1] - fp[n1][1]) +
300  SQR(fp[n2][2] - fp[n1][2])));
301  sina = sqrt(1 - cosa * cosa);
302  s[j] = sqrt(SQR(fp[n2][0] - fp[n1][0]) + SQR(fp[n2][1] - fp[n1][1]) +
303  SQR(fp[n2][2] - fp[n1][2])) *
304  sqrt(SQR(point[0] - fp[n1][0]) + SQR(point[1] - fp[n1][1]) +
305  SQR(point[2] - fp[n1][2])) *
306  sina / 2.0;
307  }
308  cosa = ((fp[2][0] - fp[0][0]) * (fp[1][0] - fp[0][0]) +
309  (fp[2][1] - fp[0][1]) * (fp[1][1] - fp[0][1]) +
310  (fp[2][2] - fp[0][2]) * (fp[1][2] - fp[0][2])) /
311  (sqrt(SQR(fp[2][0] - fp[0][0]) + SQR(fp[2][1] - fp[0][1]) +
312  SQR(fp[2][2] - fp[0][2])) *
313  sqrt(SQR(fp[1][0] - fp[0][0]) + SQR(fp[1][1] - fp[0][1]) +
314  SQR(fp[1][2] - fp[0][2])));
315  sina = sqrt(1 - cosa * cosa);
316  ss = sqrt(SQR(fp[1][0] - fp[0][0]) + SQR(fp[1][1] - fp[0][1]) +
317  SQR(fp[1][2] - fp[0][2])) *
318  sqrt(SQR(fp[2][0] - fp[0][0]) + SQR(fp[2][1] - fp[0][1]) +
319  SQR(fp[2][2] - fp[0][2])) *
320  sina / 2.0;
321  if (flag_inside == 0) {
322  /* if(fabs(ss-s[0]-s[1]-s[2])<EPSILON*ss/100.0)
323  */
324  if (fabs(ss - s[0] - s[1] - s[2]) < EPSILON * ss / 100.0)
325  flag_inside = 1;
326  }
327  /*this intersection point on the patch*/
328  if (flag_inside == 1) {
329  flag = 1;
330  for (j = 0; j < 3; j++) {
331  dis[j] = sqrt(SQR(point[0] - fp[j][0]) + SQR(point[1] - fp[j][1]) +
332  SQR(point[2] - fp[j][2]));
333  if (dis[j] < EPSILON) {
334  flag = 0;
335  color = vd->surface[index].surf_data[m * 12 + 9 + j];
336  break;
337  }
338  }
339  if (flag == 1) {
340  dist = 1.0 / (1.0 / dis[0] + 1.0 / dis[1] + 1.0 / dis[2]);
341  color = vd->surface[index].surf_data[m * 12 + 9] * dist / dis[0] +
342  vd->surface[index].surf_data[m * 12 + 10] * dist / dis[1] +
343  vd->surface[index].surf_data[m * 12 + 11] * dist / dis[2];
344  }
345 
346  for (j = 0; j < 3; j++) {
347  surf_p1[num_surp * 7 + j] = point[j];
348  surf_p1[num_surp * 7 + 3 + j] = n_f[j];
349  surf_p1[num_surp * 7 + 6] = color;
350  }
351  num_surp++;
352  }
353  }
354  } /* end of if fabs >epsilon */
355  } /*end of m*/
356  /*delete overlapped point*/
357  num_surp1 = 0;
358  if (num_surp > 0) {
359  dd = (double *)HECMW_calloc(num_surp, sizeof(double));
360  o_flag = (int *)HECMW_calloc(num_surp, sizeof(int));
361  if ((dd == NULL) || (o_flag == NULL)) {
362  fprintf(stderr, "There is no enough memory for dd and overlap_flag\n");
363  exit(0);
364  }
365  for (i = 0; i < num_surp; i++)
366  dd[i] = sqrt(SQR(surf_p1[i * 7 + 0] - in[0]) +
367  SQR(surf_p1[i * 7 + 1] - in[1]) +
368  SQR(surf_p1[i * 7 + 2] - in[2]));
369  dd_norm =
370  sqrt(SQR(out[0] - in[0]) + SQR(out[1] - in[1]) + SQR(out[2] - in[2]));
371  for (i = 0; i < num_surp - 1; i++)
372  for (j = i + 1; j < num_surp; j++) {
373  if (dd[i] > dd[j]) {
374  tmp_dd = dd[i];
375  dd[i] = dd[j];
376  dd[j] = tmp_dd;
377  for (k = 0; k < 7; k++) tmp_p[k] = surf_p1[i * 7 + k];
378  for (k = 0; k < 7; k++) surf_p1[i * 7 + k] = surf_p1[j * 7 + k];
379  for (k = 0; k < 7; k++) surf_p1[j * 7 + k] = tmp_p[k];
380  }
381  }
382  for (i = 0; i < num_surp; i++) o_flag[i] = 1;
383  for (i = 1; i < num_surp; i++) {
384  if (fabs(dd[i] - dd[i - 1]) < EPSILON * dd_norm * 100.0) o_flag[i] = 0;
385  }
386  /* for(i=1;i<num_surp;i++)
387 o_flag[i]=0;
388  */
389  num_surp1 = 0;
390  for (i = 0; i < num_surp; i++) {
391  if (o_flag[i] == 1) {
392  for (j = 0; j < 7; j++) surf_p[num_surp1 * 7 + j] = surf_p1[i * 7 + j];
393  num_surp1++;
394  }
395  }
396  HECMW_free(dd);
397  HECMW_free(o_flag);
398  HECMW_free(surf_p1);
399  }
400  return (num_surp1);
401 }
402 
403 #endif
404 
405 void compute_color_vr(int current_ijk[3], int color_mapping_style,
406  double *interval_point, int transfer_function_style,
407  double opa_value, int num_of_features, double *fea_point,
408  double view_point_d[3], int interval_mapping_num,
409  int color_system_type, int num_of_lights,
410  double *light_point, double k_ads[3], int r_level[3],
411  double orig_xyz[3], double r_dxyz[3], double *var,
412  double *grad_var, double accum_rgba[4], double mincolor,
413  double maxcolor, double grad_minmax[2],
414  double feap_minmax[2], double feai_minmax[2],
415  double dis_minmax[2], double *opa_table,
416  double in_point[3], double out_point[3],
417  double tav_length, int time_step, int print_flag) {
418  int i, j, k;
419  double a_current;
420  double color[3];
421  double value;
422  double length, coff_i;
423  int index;
424  int m;
425  double value2[8], vv[8 * 3], v;
426  double dis[8], d;
427  int return_flag;
428  double t;
429  double r, g, b;
430 
431  length =
432  sqrt(SQR(out_point[0] - in_point[0]) + SQR(out_point[1] - in_point[1]) +
433  SQR(out_point[2] - in_point[2]));
434  coff_i = length / tav_length;
435  index = current_ijk[2] * r_level[0] * r_level[1] +
436  current_ijk[1] * r_level[0] + current_ijk[0];
437 
438  /*---------------------start computing value of in_point
439  * --------------------*/
440  return_flag = 0;
441  vv[0 * 3] = vv[4 * 3] = vv[7 * 3] = vv[3 * 3] =
442  orig_xyz[0] + current_ijk[0] * r_dxyz[0];
443  vv[1 * 3] = vv[5 * 3] = vv[6 * 3] = vv[2 * 3] =
444  orig_xyz[0] + (current_ijk[0] + 1) * r_dxyz[0];
445  vv[0 * 3 + 2] = vv[1 * 3 + 2] = vv[2 * 3 + 2] = vv[3 * 3 + 2] =
446  orig_xyz[2] + current_ijk[2] * r_dxyz[2];
447  vv[4 * 3 + 2] = vv[5 * 3 + 2] = vv[6 * 3 + 2] = vv[7 * 3 + 2] =
448  orig_xyz[2] + (current_ijk[2] + 1) * r_dxyz[2];
449  vv[0 * 3 + 1] = vv[1 * 3 + 1] = vv[5 * 3 + 1] = vv[4 * 3 + 1] =
450  orig_xyz[1] + current_ijk[1] * r_dxyz[1];
451  vv[3 * 3 + 1] = vv[2 * 3 + 1] = vv[6 * 3 + 1] = vv[7 * 3 + 1] =
452  orig_xyz[1] + (current_ijk[1] + 1) * r_dxyz[1];
453 
454  for (m = 0; m < 8; m++) {
455  i = m % 2;
456  j = (m / 2) % 2;
457  k = (m / 2) / 2;
458  index = (current_ijk[2] + k) * (r_level[0] + 1) * (r_level[1] + 1) +
459  (current_ijk[1] + j) * (r_level[0] + 1) + current_ijk[0] + i;
460  value2[k * 4 + j * 2 + i] = var[index];
461  }
462 
463  for (i = 0; i < 8; i++) {
464  dis[i] =
465  sqrt(SQR(in_point[0] - vv[i * 3]) + SQR(in_point[1] - vv[i * 3 + 1]) +
466  SQR(in_point[2] - vv[i * 3 + 2]));
467  }
468  /* if(return_flag==0) {
469  */
470  d = 0.0;
471  for (i = 0; i < 8; i++) d += 1.0 / (dis[i] + EPSILON);
472  v = 0.0;
473  for (i = 0; i < 8; i++) {
474  v += value2[i] / ((dis[i] + EPSILON) * d);
475  }
476  /* }*/
477  /*---------------------end computing value of in_point --------------------*/
478 
479  value = v;
480 
481 #ifdef transfer_change
482 
483 /* if(transfer_function_style==1)
484  opacity=opa_value;
485 if(transfer_function_style==2) {
486  grad=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2]));
487  opacity=(grad-grad_minmax[0])/(grad_minmax[1]-grad_minmax[0])/200.0+0.0002;
488 }
489 
490  */
491 #endif
492  a_current = opa_value;
493  if (transfer_function_style == 3) {
494  a_current = opa_value;
495 
496  for (i = 0; i < num_of_features; i++) {
497  t = fabs(value - fea_point[i * 3]);
498  if (t < fea_point[i * 3 + 1])
499  a_current += fea_point[i * 3 + 2] * (fea_point[i * 3 + 1] - t) /
500  fea_point[i * 3 + 1];
501  }
502  }
503 
504  if (transfer_function_style == 4) {
505  a_current = opa_value;
506  for (i = 0; i < num_of_features; i++) {
507  if ((value >= fea_point[i * 3]) && (value <= fea_point[i * 3 + 1])) {
508  a_current = fea_point[i * 3 + 2];
509  }
510  }
511  }
512 
513  /*--------------map value to rgb ------------------- */
514  if (color_mapping_style == 1) {
515  if (fabs(maxcolor - mincolor) > EPSILON)
516  value = (value - mincolor) / (maxcolor - mincolor);
517  }
518 
519  if (color_mapping_style == 2) {
520  mincolor = interval_point[0];
521  maxcolor = interval_point[1];
522  if (fabs(maxcolor - mincolor) > EPSILON)
523  value = (value - mincolor) / (maxcolor - mincolor);
524  }
525 
526  if (color_mapping_style == 3) {
527  if (value < interval_point[0])
528  value = 0.0;
529  else if (value > interval_point[interval_mapping_num * 2])
530  value = 1.0;
531  else {
532  for (i = 1; i < interval_mapping_num + 1; i++) {
533  if ((value <= interval_point[i * 2]) &&
534  (value > interval_point[(i - 1) * 2])) {
535  value = (value - interval_point[(i - 1) * 2]) /
536  (interval_point[i * 2] - interval_point[(i - 1) * 2]) *
537  (interval_point[i * 2 + 1] -
538  interval_point[(i - 1) * 2 + 1]) +
539  interval_point[(i - 1) * 2 + 1];
540  break;
541  }
542  }
543  }
544  }
545 
546  if (color_system_type == 1) {
547  if (value < 0.0) value = 0.0;
548  if (value > 1.0) value = 1.0;
549  if (value <= 0.25) {
550  r = 0.0;
551  g = value * 4.0;
552  b = 1.0;
553  } else if ((value > 0.25) && (value <= 0.5)) {
554  r = 0.0;
555  g = 1.0;
556  b = (0.5 - value) * 4.0;
557  } else if ((value > 0.5) && (value <= 0.75)) {
558  r = (value - 0.5) * 4.0;
559  g = 1.0;
560  b = 0.0;
561  } else if (value > 0.75) {
562  r = 1.0;
563  g = (1.0 - value) * 4.0;
564  b = 0.0;
565  }
566 
567  }
568 
569  else if (color_system_type == 2) {
570  if (value < 0.0) value = 0.0;
571  if (value > 1.0) value = 1.0;
572  if (value <= 0.2) {
573  g = 0.0;
574  b = 1.0;
575  r = (0.2 - value) * 5.0;
576  } else if ((value > 0.2) && (value <= 0.4)) {
577  r = 0.0;
578  b = 1.0;
579  g = (value - 0.2) * 5.0;
580  } else if ((value > 0.4) && (value <= 0.6)) {
581  r = 0.0;
582  g = 1.0;
583  b = 1.0 - (value - 0.4) * 5.0;
584  } else if ((value > 0.6) && (value <= 0.8)) {
585  r = (value - 0.6) * 5.0;
586  g = 1.0;
587  b = 0.0;
588  } else if (value > 0.0) {
589  r = 1.0;
590  g = 1.0 - (value - 0.8) * 5.0;
591  b = 0.0;
592  }
593  } else if (color_system_type == 3) {
594  r = g = b = value;
595  }
596 
597  color[0] = r;
598  color[1] = g;
599  color[2] = b;
600  /* --------------end of mapping value to rgb ----------------------- */
601 
602  r = 0.0;
603  g = 0.0;
604  b = 0.0;
605 /* for(j=0;j<num_of_lights;j++) {
606  */
607 #ifdef slow
608  for (i = 0; i < 3; i++) {
609  lp[i] = light_point[j * 3 + i] - p[i];
610  vp[i] = view_point_d[i] - p[i];
611  hp[i] = (lp[i] + vp[i]) / 2.0;
612  }
613  lp_norm = sqrt(SQR(lp[0]) + SQR(lp[1]) + SQR(lp[2]));
614  vp_norm = sqrt(SQR(vp[0]) + SQR(vp[1]) + SQR(vp[2]));
615  hp_norm = sqrt(SQR(hp[0]) + SQR(hp[1]) + SQR(hp[2]));
616  if (fabs(lp_norm) > EPSILON) {
617  for (i = 0; i < 3; i++) lp[i] /= lp_norm;
618  }
619  if (fabs(vp_norm) > EPSILON) {
620  for (i = 0; i < 3; i++) vp[i] /= vp_norm;
621  }
622  if (fabs(hp_norm) > EPSILON) {
623  for (i = 0; i < 3; i++) hp[i] /= hp_norm;
624  }
625  norm = sqrt(SQR(n[0]) + SQR(n[1]) + SQR(n[2]));
626  if (fabs(norm) > EPSILON) {
627  for (i = 0; i < 3; i++) n[i] /= norm;
628  }
629  inprodLN = n[0] * lp[0] + n[1] * lp[1] + n[2] * lp[2];
630  inprodVN = n[0] * vp[0] + n[1] * vp[1] + n[2] * vp[2];
631 
632  inprodHN = n[0] * hp[0] + n[1] * hp[1] + n[2] * hp[2];
633  /* a_current=opacity_decision(in_voxel, vr, vd, grad_minmax, feap_minmax,
634  feai_minmax,
635  dis_minmax, opa_table, mincolor, maxcolor, time_step);
636  cosalpha=sqrt(1.0-pow(inprodLN,2));
637  */
638  cosalpha = inprodLN;
639  costheta = inprodLN * inprodVN -
640  sqrt(1.0 - inprodLN * inprodLN) * sqrt(1.0 - inprodVN * inprodVN);
641 
642  cosalpha = fabs(cosalpha);
643 #endif
644 
645  /* r=color[0]*(ka+kd*cosalpha+ks*pow(costheta,6));
646  g=color[1]*(ka+kd*cosalpha+ks*pow(costheta,6));
647  b=color[2]*(ka+kd*cosalpha+ks*pow(costheta,6));
648  */
649 
650  r += color[0] * k_ads[0] * coff_i;
651  g += color[1] * k_ads[0] * coff_i;
652  b += color[2] * k_ads[0] * coff_i;
653 
654  /*
655  }
656  */
657  r *= a_current;
658  g *= a_current;
659  b *= a_current;
660  if (accum_rgba[3] < 0.99) {
661  accum_rgba[0] += r * (1.0 - accum_rgba[3]);
662  accum_rgba[1] += g * (1.0 - accum_rgba[3]);
663  accum_rgba[2] += b * (1.0 - accum_rgba[3]);
664  accum_rgba[3] += a_current * (1.0 - accum_rgba[3]) * coff_i;
665  }
666 
667  return;
668 }
#define NULL
#define HECMW_calloc(nmemb, size)
Definition: hecmw_malloc.h:21
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
int find_surface_point(double fp[3][3], double point_o[3], double view_point_d[3], double point[3], double n_f[4], double v_f[9], double c_value[3], double *value, double normal[9], int normal_flag, int smooth_shading)
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
#define SQR(x)
int * level