FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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>
10#include "hecmw_malloc.h"
11
12#if 0
13void 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
26double 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/*
72void 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
92double 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
230int 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++)
387o_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
405void 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;
485if(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