13void find_sample_point(
double in_point[3],
double out_point[3],
int num_sample,
double *sample_point)
18 t=1.0/(double)(num_sample+1);
19 for(i=0;i<num_sample;i++) {
21 sample_point[i*3+j]=in_point[j]+t*(i+1)*(out_point[j]-in_point[j]);
26double value_compute(
double point[3],
int current_ijk[3],
double orig_xyz[3],
double r_dxyz[3],
int r_level[3],
double *var)
29 double value[8], vv[8*3], v;
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];
45 index=(current_ijk[2]+k)*(r_level[0]+1)*(r_level[1]+1)+(current_ijk[1]+j)*(r_level[0]+1)+
47 value[k*4+j*2+i]=var[index];
51 dis[i]=sqrt(
SQR(point[0]-vv[i*3])+
SQR(point[1]-vv[i*3+1])+
SQR(point[2]-vv[i*3+2]));
63 v+=value[i]/(dis[i]*d);
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)
101 int i,
level, min_type;
102 double cp[3], del_l, dist;
105 if(transfer_function_style==1)
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;
111 else if(transfer_function_style==3) {
115 for(i=0;i<num_of_features;i++) {
116 t=fabs(value- fea_point[i*3]);
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;
130 else opacity=opa_value;
132 else if(transfer_function_style==4) {
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];
159 else if(transfer_function_style==5) {
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;
170 else if(transfer_function_style==6) {
176 cp[i]=(current_ijk[i]+0.5)*r_dxyz[i]+orig_xyz[i];
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;
182 else if(transfer_function_style==7) {
185 if(color_mapping_style==1) {
186 if(fabs(maxcolor-mincolor)>
EPSILON)
187 value=(value-mincolor)/(maxcolor-mincolor);
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);
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];
206 if(value>1.0) value=1.0;
207 if(value<0.0) value=0.0;
214 opacity=value/200.0+0.0002;;
215 if(opacity<0) opacity=0.0;
217 else if(transfer_function_style==8) {
218 del_l=(maxcolor-mincolor)/255.0;
219 level=(int)(value-mincolor)/del_l;
222 opacity=opa_table[
level];
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;
239 double tmp_dd, tmp_p[7];
243 if (surf_p1 ==
NULL) {
244 fprintf(stderr,
"There is no enough memory for surf_p\n"), exit(0);
247 for (m = 0; m < vd->
surface[index].
num; m++) {
250 for (i = 0; i < 3; i++)
251 for (j = 0; j < 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]);
261 for (j = 0; j < 3; j++) n_f[j] /= n_norm;
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]);
266 for (i = 0; i < 3; i++) v_f[i] /= v_norm;
268 nv_sign = n_f[0] * v_f[0] + n_f[1] * v_f[1] + n_f[2] * v_f[2];
270 for (i = 0; i < 3; i++) n_f[i] = -n_f[i];
272 n_f[3] = -n_f[0] * fp[0][0] - n_f[1] * fp[0][1] - n_f[2] * fp[0][2];
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]);
283 for (j = 0; j < 3; j++) {
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)) {
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])) *
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])) *
321 if (flag_inside == 0) {
324 if (fabs(ss - s[0] - s[1] - s[2]) <
EPSILON * ss / 100.0)
328 if (flag_inside == 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]));
340 dist = 1.0 / (1.0 / dis[0] + 1.0 / dis[1] + 1.0 / dis[2]);
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;
361 if ((dd ==
NULL) || (o_flag ==
NULL)) {
362 fprintf(stderr,
"There is no enough memory for dd and overlap_flag\n");
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]));
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++) {
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];
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;
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];
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) {
422 double length, coff_i;
425 double value2[8], vv[8 * 3], v;
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];
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];
454 for (m = 0; m < 8; m++) {
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];
463 for (i = 0; i < 8; 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]));
471 for (i = 0; i < 8; i++) d += 1.0 / (dis[i] +
EPSILON);
473 for (i = 0; i < 8; i++) {
474 v += value2[i] / ((dis[i] +
EPSILON) * d);
481#ifdef transfer_change
492 a_current = opa_value;
493 if (transfer_function_style == 3) {
494 a_current = opa_value;
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];
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];
514 if (color_mapping_style == 1) {
515 if (fabs(maxcolor - mincolor) >
EPSILON)
516 value = (value - mincolor) / (maxcolor - mincolor);
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);
526 if (color_mapping_style == 3) {
527 if (value < interval_point[0])
529 else if (value > interval_point[interval_mapping_num * 2])
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];
546 if (color_system_type == 1) {
547 if (value < 0.0) value = 0.0;
548 if (value > 1.0) value = 1.0;
553 }
else if ((value > 0.25) && (value <= 0.5)) {
556 b = (0.5 - value) * 4.0;
557 }
else if ((value > 0.5) && (value <= 0.75)) {
558 r = (value - 0.5) * 4.0;
561 }
else if (value > 0.75) {
563 g = (1.0 - value) * 4.0;
569 else if (color_system_type == 2) {
570 if (value < 0.0) value = 0.0;
571 if (value > 1.0) value = 1.0;
575 r = (0.2 - value) * 5.0;
576 }
else if ((value > 0.2) && (value <= 0.4)) {
579 g = (value - 0.2) * 5.0;
580 }
else if ((value > 0.4) && (value <= 0.6)) {
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;
588 }
else if (value > 0.0) {
590 g = 1.0 - (value - 0.8) * 5.0;
593 }
else if (color_system_type == 3) {
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;
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]));
617 for (i = 0; i < 3; i++) lp[i] /= lp_norm;
620 for (i = 0; i < 3; i++) vp[i] /= vp_norm;
623 for (i = 0; i < 3; i++) hp[i] /= hp_norm;
625 norm = sqrt(
SQR(n[0]) +
SQR(n[1]) +
SQR(n[2]));
627 for (i = 0; i < 3; i++) n[i] /= norm;
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];
632 inprodHN = n[0] * hp[0] + n[1] * hp[1] + n[2] * hp[2];
639 costheta = inprodLN * inprodVN -
640 sqrt(1.0 - inprodLN * inprodLN) * sqrt(1.0 - inprodVN * inprodVN);
642 cosalpha = fabs(cosalpha);
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;
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;
#define HECMW_calloc(nmemb, size)
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)