FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
11int 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
132int 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
182int 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
229void 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,
364feai_minmax,
365 dis_minmax, opa_table, mincolor, maxcolor, time_step);
366cosalpha=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)