FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_vis_color_composite_sf.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_SF_geom.h"
10 
11 int find_surface_point(double fp[3][3], double point_o[3],
12  double view_point_d[3], double point[3], double n_f[4],
13  double v_f[9], double c_value[3], double *value,
14  double normal[9], int normal_flag, int smooth_shading) {
15  int i, j, n1, n2, flag_inside;
16  double s[3], sina, cosa, ss, t, dis[3], dist, n_norm, v_norm, nv_sign;
17  double tmp_p[7];
18 
19  /*find the equation of the patch first */
20  flag_inside = 0;
21  if (normal_flag == 1) {
22  n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
23  (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
24  n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
25  (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
26  n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
27  (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
28  n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
29  if (fabs(n_norm) > EPSILON) {
30  for (j = 0; j < 3; j++) n_f[j] /= n_norm;
31  }
32  }
33  for (i = 0; i < 3; i++) v_f[i] = view_point_d[i] - point_o[i];
34  v_norm = sqrt(v_f[0] * v_f[0] + v_f[1] * v_f[1] + v_f[2] * v_f[2]);
35  if (fabs(v_norm) > EPSILON) {
36  for (i = 0; i < 3; i++) v_f[i] /= v_norm;
37  }
38  nv_sign = n_f[0] * v_f[0] + n_f[1] * v_f[1] + n_f[2] * v_f[2];
39  if (nv_sign < 0.0) {
40  for (i = 0; i < 3; i++) n_f[i] = -n_f[i];
41  }
42  if (normal_flag == 1)
43  n_f[3] = -n_f[0] * fp[0][0] - n_f[1] * fp[0][1] - n_f[2] * fp[0][2];
44  if (smooth_shading == 0)
45  for (i = 0; i < 3; i++) normal[i] = n_f[i];
46  /*find intersection point*/
47  if (fabs(n_f[0] * (point_o[0] - view_point_d[0]) +
48  n_f[1] * (point_o[1] - view_point_d[1]) +
49  n_f[2] * (point_o[2] - view_point_d[2])) > EPSILON) {
50  t = (-n_f[3] - n_f[0] * view_point_d[0] - n_f[1] * view_point_d[1] -
51  n_f[2] * view_point_d[2]) /
52  (n_f[0] * (point_o[0] - view_point_d[0]) +
53  n_f[1] * (point_o[1] - view_point_d[1]) +
54  n_f[2] * (point_o[2] - view_point_d[2]));
55  if (t > -EPSILON) {
56  for (j = 0; j < 3; j++)
57  point[j] = view_point_d[j] + t * (point_o[j] - view_point_d[j]);
58  /*judge whether it is inside the range of patch by equal area method*/
59  for (j = 0; j < 3; j++) {
60  n1 = j;
61  n2 = j + 1;
62  if (n2 == 3) n2 = 0;
63  if ((fabs(point[0] - fp[n1][0]) < EPSILON) &&
64  (fabs(point[1] - fp[n1][1]) < EPSILON) &&
65  (fabs(point[2] - fp[n1][2]) < EPSILON)) {
66  flag_inside = 1;
67  *value = c_value[n1];
68  if (smooth_shading == 1) {
69  tmp_p[0] = normal[n1 * 3 + 0];
70  tmp_p[1] = normal[n1 * 3 + 1];
71  tmp_p[2] = normal[n1 * 3 + 2];
72  normal[0] = tmp_p[0];
73  normal[1] = tmp_p[1];
74  normal[2] = tmp_p[2];
75  }
76  break;
77  }
78  cosa = ((point[0] - fp[n1][0]) * (fp[n2][0] - fp[n1][0]) +
79  (point[1] - fp[n1][1]) * (fp[n2][1] - fp[n1][1]) +
80  (point[2] - fp[n1][2]) * (fp[n2][2] - fp[n1][2])) /
81  (sqrt(SQR(point[0] - fp[n1][0]) + SQR(point[1] - fp[n1][1]) +
82  SQR(point[2] - fp[n1][2])) *
83  sqrt(SQR(fp[n2][0] - fp[n1][0]) + SQR(fp[n2][1] - fp[n1][1]) +
84  SQR(fp[n2][2] - fp[n1][2])));
85  sina = sqrt(1 - cosa * cosa);
86  s[j] = sqrt(SQR(fp[n2][0] - fp[n1][0]) + SQR(fp[n2][1] - fp[n1][1]) +
87  SQR(fp[n2][2] - fp[n1][2])) *
88  sqrt(SQR(point[0] - fp[n1][0]) + SQR(point[1] - fp[n1][1]) +
89  SQR(point[2] - fp[n1][2])) *
90  sina / 2.0;
91  }
92  if (flag_inside == 0) {
93  cosa = ((fp[2][0] - fp[0][0]) * (fp[1][0] - fp[0][0]) +
94  (fp[2][1] - fp[0][1]) * (fp[1][1] - fp[0][1]) +
95  (fp[2][2] - fp[0][2]) * (fp[1][2] - fp[0][2])) /
96  (sqrt(SQR(fp[2][0] - fp[0][0]) + SQR(fp[2][1] - fp[0][1]) +
97  SQR(fp[2][2] - fp[0][2])) *
98  sqrt(SQR(fp[1][0] - fp[0][0]) + SQR(fp[1][1] - fp[0][1]) +
99  SQR(fp[1][2] - fp[0][2])));
100  sina = sqrt(1 - cosa * cosa);
101  ss = sqrt(SQR(fp[1][0] - fp[0][0]) + SQR(fp[1][1] - fp[0][1]) +
102  SQR(fp[1][2] - fp[0][2])) *
103  sqrt(SQR(fp[2][0] - fp[0][0]) + SQR(fp[2][1] - fp[0][1]) +
104  SQR(fp[2][2] - fp[0][2])) *
105  sina / 2.0;
106  /* if(fabs(ss-s[0]-s[1]-s[2])<EPSILON*ss/100.0)
107  */
108  if (fabs(ss - s[0] - s[1] - s[2]) < EPSILON) {
109  flag_inside = 1;
110  for (j = 0; j < 3; j++)
111  dis[j] = sqrt(SQR(point[0] - fp[j][0]) + SQR(point[1] - fp[j][1]) +
112  SQR(point[2] - fp[j][2]));
113  dist = 1.0 / dis[0] + 1.0 / dis[1] + 1.0 / dis[2];
114  *value = c_value[0] / (dis[0] * dist) + c_value[1] / (dis[1] * dist) +
115  c_value[2] / (dis[2] * dist);
116  if (smooth_shading == 1) {
117  for (j = 0; j < 3; j++) {
118  tmp_p[j] = normal[0 * 3 + j] / (dis[0] * dist) +
119  normal[1 * 3 + j] / (dis[1] * dist) +
120  normal[2 * 3 + j] / (dis[2] * dist);
121  }
122  for (j = 0; j < 3; j++) normal[j] = tmp_p[j];
123  }
124  }
125  }
126  }
127  }
128 
129  return (flag_inside);
130 }
131 
132 int find_point_depth(double fp[3][3], double point_o[3], double view_point_d[3],
133  double n_f[4], double point[3], int normal_flag) {
134  int i, j, flag_inside;
135  double v_f[3], t, n_norm, v_norm, nv_sign;
136 
137  /*find the equation of the patch first */
138  flag_inside = 0;
139  if (normal_flag == 1) {
140  n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
141  (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
142  n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
143  (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
144  n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
145  (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
146  n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
147  if (fabs(n_norm) > EPSILON) {
148  for (j = 0; j < 3; j++) n_f[j] /= n_norm;
149  }
150  }
151  for (i = 0; i < 3; i++) v_f[i] = view_point_d[i] - point_o[i];
152  v_norm = sqrt(v_f[0] * v_f[0] + v_f[1] * v_f[1] + v_f[2] * v_f[2]);
153  if (fabs(v_norm) > EPSILON) {
154  for (i = 0; i < 3; i++) v_f[i] /= v_norm;
155  }
156  nv_sign = n_f[0] * v_f[0] + n_f[1] * v_f[1] + n_f[2] * v_f[2];
157  if (nv_sign < 0.0) {
158  for (i = 0; i < 3; i++) n_f[i] = -n_f[i];
159  }
160  if (normal_flag == 1)
161  n_f[3] = -n_f[0] * fp[0][0] - n_f[1] * fp[0][1] - n_f[2] * fp[0][2];
162 
163  /*find intersection point*/
164  if (fabs(n_f[0] * (point_o[0] - view_point_d[0]) +
165  n_f[1] * (point_o[1] - view_point_d[1]) +
166  n_f[2] * (point_o[2] - view_point_d[2])) > EPSILON) {
167  t = (-n_f[3] - n_f[0] * view_point_d[0] - n_f[1] * view_point_d[1] -
168  n_f[2] * view_point_d[2]) /
169  (n_f[0] * (point_o[0] - view_point_d[0]) +
170  n_f[1] * (point_o[1] - view_point_d[1]) +
171  n_f[2] * (point_o[2] - view_point_d[2]));
172  if (t > -EPSILON) {
173  for (j = 0; j < 3; j++)
174  point[j] = view_point_d[j] + t * (point_o[j] - view_point_d[j]);
175  flag_inside = 1;
176  }
177  }
178 
179  return (flag_inside);
180 }
181 
182 int find2_point_depth(double fp[3][3], double point_o[3],
183  double view_point_d[3], double point[3]) {
184  int i, j, flag_inside;
185  double v_f[3], t, n_norm, v_norm, nv_sign;
186  double n_f[4];
187 
188  /*find the equation of the patch first */
189  flag_inside = 0;
190  n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
191  (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
192  n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
193  (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
194  n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
195  (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
196  n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
197  if (fabs(n_norm) > EPSILON) {
198  for (j = 0; j < 3; j++) n_f[j] /= n_norm;
199  }
200  for (i = 0; i < 3; i++) v_f[i] = view_point_d[i] - point_o[i];
201  v_norm = sqrt(v_f[0] * v_f[0] + v_f[1] * v_f[1] + v_f[2] * v_f[2]);
202  if (fabs(v_norm) > EPSILON) {
203  for (i = 0; i < 3; i++) v_f[i] /= v_norm;
204  }
205  nv_sign = n_f[0] * v_f[0] + n_f[1] * v_f[1] + n_f[2] * v_f[2];
206  if (nv_sign < 0.0) {
207  for (i = 0; i < 3; i++) n_f[i] = -n_f[i];
208  }
209  n_f[3] = -n_f[0] * fp[0][0] - n_f[1] * fp[0][1] - n_f[2] * fp[0][2];
210  /*find intersection point*/
211  if (fabs(n_f[0] * (point_o[0] - view_point_d[0]) +
212  n_f[1] * (point_o[1] - view_point_d[1]) +
213  n_f[2] * (point_o[2] - view_point_d[2])) > EPSILON) {
214  t = (-n_f[3] - n_f[0] * view_point_d[0] - n_f[1] * view_point_d[1] -
215  n_f[2] * view_point_d[2]) /
216  (n_f[0] * (point_o[0] - view_point_d[0]) +
217  n_f[1] * (point_o[1] - view_point_d[1]) +
218  n_f[2] * (point_o[2] - view_point_d[2]));
219  if (t > -EPSILON) {
220  for (j = 0; j < 3; j++)
221  point[j] = view_point_d[j] + t * (point_o[j] - view_point_d[j]);
222  flag_inside = 1;
223  }
224  }
225 
226  return (flag_inside);
227 }
228 
229 void compute_color_sf(double p[3], double value, double n[3],
230  int color_mapping_style, double *interval_point,
231  double view_point_d[3], int interval_mapping_num,
232  int color_system_type, int num_of_lights,
233  double *light_point, double k_ads[3],
234  double accum_rgba[4], double mincolor, double maxcolor,
235  int display_method) {
236  int i, j;
237  double cosalpha, costheta;
238  double color[3];
239  double lp[3], vp[3], lp_norm, vp_norm, norm, hp[3], hp_norm;
240  double inprodLN, inprodVN, inprodHN;
241  double r, g, b;
242 
243  /*--------------map value to rgb ------------------- */
244  if (display_method != 4) {
245  if (color_mapping_style == 1) {
246  if (fabs(maxcolor - mincolor) > EPSILON)
247  value = (value - mincolor) / (maxcolor - mincolor);
248  }
249 
250  if (color_mapping_style == 2) {
251  mincolor = interval_point[0];
252  maxcolor = interval_point[1];
253  if (fabs(maxcolor - mincolor) > EPSILON)
254  value = (value - mincolor) / (maxcolor - mincolor);
255  }
256 
257  if ((color_mapping_style == 3) || (color_mapping_style == 4)) {
258  if (value < interval_point[0])
259  value = 0.0;
260  else if (value > interval_point[interval_mapping_num * 2])
261  value = 1.0;
262  else {
263  for (i = 1; i < interval_mapping_num + 1; i++) {
264  if ((value <= interval_point[i * 2]) &&
265  (value > interval_point[(i - 1) * 2])) {
266  value = (value - interval_point[(i - 1) * 2]) /
267  (interval_point[i * 2] - interval_point[(i - 1) * 2]) *
268  (interval_point[i * 2 + 1] -
269  interval_point[(i - 1) * 2 + 1]) +
270  interval_point[(i - 1) * 2 + 1];
271  break;
272  }
273  }
274  }
275  }
276  }
277 
278  if (color_system_type == 1) {
279  if (value < 0.0) value = 0.0;
280  if (value > 1.0) value = 1.0;
281  if (value <= 0.25) {
282  r = 0.0;
283  g = value * 4.0;
284  b = 1.0;
285  } else if ((value > 0.25) && (value <= 0.5)) {
286  r = 0.0;
287  g = 1.0;
288  b = (0.5 - value) * 4.0;
289  } else if ((value > 0.5) && (value <= 0.75)) {
290  r = (value - 0.5) * 4.0;
291  g = 1.0;
292  b = 0.0;
293  } else if (value > 0.75) {
294  r = 1.0;
295  g = (1.0 - value) * 4.0;
296  b = 0.0;
297  }
298 
299  }
300 
301  else if (color_system_type == 2) {
302  if (value < 0.0) value = 0.0;
303  if (value > 1.0) value = 1.0;
304  if (value <= 0.2) {
305  g = 0.0;
306  b = 1.0;
307  r = (0.2 - value) * 5.0;
308  } else if ((value > 0.2) && (value <= 0.4)) {
309  r = 0.0;
310  b = 1.0;
311  g = (value - 0.2) * 5.0;
312  } else if ((value > 0.4) && (value <= 0.6)) {
313  r = 0.0;
314  g = 1.0;
315  b = 1.0 - (value - 0.4) * 5.0;
316  } else if ((value > 0.6) && (value <= 0.8)) {
317  r = (value - 0.6) * 5.0;
318  g = 1.0;
319  b = 0.0;
320  } else if (value > 0.0) {
321  r = 1.0;
322  g = 1.0 - (value - 0.8) * 5.0;
323  b = 0.0;
324  }
325  } else if (color_system_type == 3) {
326  r = g = b = value;
327  }
328 
329  color[0] = r;
330  color[1] = g;
331  color[2] = b;
332  /* --------------end of mapping value to rgb ----------------------- */
333 
334  r = 0.0;
335  g = 0.0;
336  b = 0.0;
337  for (j = 0; j < num_of_lights; j++) {
338  for (i = 0; i < 3; i++) {
339  lp[i] = light_point[j * 3 + i] - p[i];
340  vp[i] = view_point_d[i] - p[i];
341  hp[i] = (lp[i] + vp[i]) / 2.0;
342  }
343  lp_norm = sqrt(SQR(lp[0]) + SQR(lp[1]) + SQR(lp[2]));
344  vp_norm = sqrt(SQR(vp[0]) + SQR(vp[1]) + SQR(vp[2]));
345  hp_norm = sqrt(SQR(hp[0]) + SQR(hp[1]) + SQR(hp[2]));
346  if (fabs(lp_norm) > EPSILON) {
347  for (i = 0; i < 3; i++) lp[i] /= lp_norm;
348  }
349  if (fabs(vp_norm) > EPSILON) {
350  for (i = 0; i < 3; i++) vp[i] /= vp_norm;
351  }
352  if (fabs(hp_norm) > EPSILON) {
353  for (i = 0; i < 3; i++) hp[i] /= hp_norm;
354  }
355  norm = sqrt(SQR(n[0]) + SQR(n[1]) + SQR(n[2]));
356  if (fabs(norm) > EPSILON) {
357  for (i = 0; i < 3; i++) n[i] /= norm;
358  }
359  inprodLN = n[0] * lp[0] + n[1] * lp[1] + n[2] * lp[2];
360  inprodVN = n[0] * vp[0] + n[1] * vp[1] + n[2] * vp[2];
361 
362  inprodHN = n[0] * hp[0] + n[1] * hp[1] + n[2] * hp[2];
363  /* a_current=opacity_decision(in_voxel, vr, vd, grad_minmax, feap_minmax,
364 feai_minmax,
365  dis_minmax, opa_table, mincolor, maxcolor, time_step);
366 cosalpha=sqrt(1.0-pow(inprodLN,2));
367  */
368  cosalpha = inprodLN;
369  costheta =
370  inprodLN * inprodVN -
371  sqrt(1.0 - inprodLN * inprodLN) * sqrt(1.0 - inprodVN * inprodVN);
372 
373  cosalpha = fabs(cosalpha);
374 
375  if (cosalpha > 0) {
376  r += color[0] * (k_ads[0] + k_ads[1] * cosalpha +
377  k_ads[2] * costheta * costheta * costheta * costheta *
378  costheta * costheta);
379  g += color[1] * (k_ads[0] + k_ads[1] * cosalpha +
380  k_ads[2] * costheta * costheta * costheta * costheta *
381  costheta * costheta);
382  b += color[2] * (k_ads[0] + k_ads[1] * cosalpha +
383  k_ads[2] * costheta * costheta * costheta * costheta *
384  costheta * costheta);
385  }
386 
387  else {
388  r += k_ads[0] * color[0];
389  g += k_ads[0] * color[1];
390  b += k_ads[0] * color[2];
391  }
392  }
393  /* r*=a_current; g*=a_current; b*=a_current;
394  if(accum_rgba[3]<0.99) {
395  accum_rgba[0]+=r*(1.0-accum_rgba[3]);
396  accum_rgba[1]+=g*(1.0-accum_rgba[3]);
397  accum_rgba[2]+=b*(1.0-accum_rgba[3]);
398  accum_rgba[3]+=a_current*(1.0-accum_rgba[3])*coff_i;
399  }
400  */
401  accum_rgba[0] = r;
402  accum_rgba[1] = g;
403  accum_rgba[2] = b;
404  accum_rgba[3] = 1.0;
405  return;
406 }
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_sf(double p[3], double value, double n[3], int color_mapping_style, double *interval_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 accum_rgba[4], double mincolor, double maxcolor, int display_method)
int find2_point_depth(double fp[3][3], double point_o[3], double view_point_d[3], double point[3])
int find_point_depth(double fp[3][3], double point_o[3], double view_point_d[3], double n_f[4], double point[3], int normal_flag)
#define EPSILON
#define SQR(x)