1 /*****************************************************************************
2 * Copyright (c) 2019 FrontISTR Commons
3 * This software is released under the MIT License, see LICENSE.txt
4 *****************************************************************************/
5
6 #include "hecmw_vis_generate_histogram_vr.h"
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "hecmw_vis_bmp.h"
12 #include "hecmw_vis_font_texture.h"
13 #include "hecmw_vis_mem_util.h"
14 #include "hecmw_malloc.h"
15
find_color_minmax_vr(double * var,int * empty_flag,int nx,int ny,int nz,double * mincolor,double * maxcolor)16 void find_color_minmax_vr(double *var, int *empty_flag, int nx, int ny, int nz,
17 double *mincolor, double *maxcolor) {
18 int i;
19 double value;
20
21 for (i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) {
22 if (empty_flag[i] == 1) {
23 value = var[i];
24 if (value < *mincolor) *mincolor = value;
25 if (value > *maxcolor) *maxcolor = value;
26 }
27 }
28 return;
29 }
generate_histogram_graph_vr(double tmincolor,double tmaxcolor,double * var,int * empty_flag,int nx,int ny,int nz,int mynode,int pesize,HECMW_Comm VIS_COMM,int color_system_type)30 void generate_histogram_graph_vr(double tmincolor, double tmaxcolor,
31 double *var, int *empty_flag, int nx, int ny,
32 int nz, int mynode, int pesize,
33 HECMW_Comm VIS_COMM, int color_system_type) {
34 int i, j, k, m;
35 double delta, value, color[3];
36 int count[500], t_count[500], max_number, max_length, start_x, end_x, start_y,
37 end_y;
38 FILE *fp;
39 double *graph;
40 BITMAPFILEHEADER header; /* File header */
41
42 unsigned char r, g, b;
43 int ri, gi, bi;
44 BITMAPINFOHEADER info;
45
46 int start_xs, start_ys;
47 char buf[128];
48 int output7[7][7];
49
50 delta = (tmaxcolor - tmincolor) / 500.0;
51 for (i = 0; i < 500; i++) {
52 count[i] = 0;
53 t_count[i] = 0;
54 }
55 for (i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) {
56 if (empty_flag[i] == 1) {
57 j = (int)((var[i] - tmincolor) / delta);
58 if (j < 0) j = 0;
59 if (j > 499) j = 499;
60 count[j]++;
61 }
62 }
63 if (pesize > 1)
64 HECMW_Allreduce(count, t_count, 500, HECMW_INT, HECMW_SUM, VIS_COMM);
65 else {
66 for (i = 0; i < 500; i++) t_count[i] = count[i];
67 }
68 if (mynode == 0) {
69 fp = fopen("histogram.bmp", "wb");
70 if (fp == NULL)
71 HECMW_vis_print_exit("Cannot generate the histogram output file");
72 graph = (double *)HECMW_calloc(400 * 530 * 3, sizeof(double));
73 if (graph == NULL) HECMW_vis_memory_exit("graph");
74 for (i = 0; i < 400 * 530 * 3; i++) graph[i] = 0.0;
75 max_number = 0;
76 for (i = 0; i < 500; i++) {
77 if (t_count[i] > max_number) max_number = t_count[i];
78 }
79 if (max_number == 0)
80 HECMW_vis_print_exit(
81 "ERROR: HEC-MW-VIS-E2003:Cannot generate histogram graph, the number "
82 "of voxels is 0");
83 max_length = (int)(400 - 30 - 5 - 45 * 1.5);
84 start_x = (int)(5 + 45 * 1.5 + 15);
85 start_y = 15;
86 for (j = 0; j < 500; j++) {
87 end_x =
88 (int)((double)t_count[j] * (double)max_length / (double)max_number +
89 start_x) +
90 2;
91 value = j / 500.0;
92 value2_to_rgb(value, color, color_system_type);
93 for (i = start_x; i < end_x; i++) {
94 graph[((j + 15) * 400 + i) * 3] = color[0];
95 graph[((j + 15) * 400 + i) * 3 + 1] = color[1];
96 graph[((j + 15) * 400 + i) * 3 + 2] = color[2];
97 }
98 }
99 /*start mark scales */
100 start_y = 15;
101 end_y = 515;
102 for (k = 0; k < 11; k++) {
103 value = tmincolor + (tmaxcolor - tmincolor) / 10.0 * k;
104 start_ys = start_y + (int)((double)500.0 / 10 * k) - (int)7 / 2;
105 start_xs = 15;
106 sprintf(buf, "%3.2E", value);
107 for (m = 0; m < 9; m++) {
108 font7_generate(buf[8 - m], output7);
109 for (j = 0; j < 7; j++)
110 for (i = 0; i < 7; i++) {
111 graph[((start_ys + j) * 400 + start_xs - i) * 3] =
112 (double)output7[6 - j][i];
113 graph[((start_ys + j) * 400 + start_xs - i) * 3 + 1] =
114 (double)output7[6 - j][i];
115 graph[((start_ys + j) * 400 + start_xs - i) * 3 + 2] =
116 (double)output7[6 - j][i];
117 }
118 start_xs += 7;
119 if ((value >= 0) && (m == 0)) start_xs -= 7;
120 }
121 if ((k != 0) && (k != 10)) {
122 start_ys += (int)7 / 2;
123
124 for (i = start_x; i < start_x + 5; i++)
125 for (j = 0; j < 3; j++) graph[(start_ys * 400 + i) * 3 + j] = 1.0;
126 }
127 }
128 header.bfSize = 54 + 3 * 400 * 530;
129 #ifdef CONVERSE_ORDER
130 header.bfSize = change_unsigned_int_order(54 + 3 * 400 * 530);
131 #endif
132 header.bfReserved1 = 0;
133 #ifdef CONVERSE_ORDER
134 header.bfReserved1 = change_short_int_order(0);
135 #endif
136
137 header.bfReserved2 = 0;
138 #ifdef CONVERSE_ORDER
139 header.bfReserved2 = change_short_int_order(0);
140 #endif
141
142 header.bfOffBits = 54;
143 #ifdef CONVERSE_ORDER
144 header.bfOffBits = change_unsigned_int_order(54);
145 #endif
146
147 info.biBitCount = 24;
148 #ifdef CONVERSE_ORDER
149 info.biBitCount = change_short_int_order(24);
150 #endif
151
152 info.biSize = 40;
153 #ifdef CONVERSE_ORDER
154 info.biSize = change_unsigned_int_order(40);
155 #endif
156
157 info.biWidth = 400;
158 #ifdef CONVERSE_ORDER
159 info.biWidth = change_int_order(400);
160 #endif
161
162 info.biHeight = 530;
163 #ifdef CONVERSE_ORDER
164 info.biHeight = change_int_order(530);
165 #endif
166
167 info.biSizeImage = 3 * 400 * 530;
168 #ifdef CONVERSE_ORDER
169 info.biSizeImage = change_unsigned_int_order(3 * 400 * 530);
170 #endif
171
172 info.biClrImportant = 0;
173 #ifdef CONVERSE_ORDER
174 info.biClrImportant = change_unsigned_int_order(0);
175 #endif
176
177 info.biClrUsed = 0;
178 #ifdef CONVERSE_ORDER
179 info.biClrUsed = change_unsigned_int_order(0);
180 #endif
181
182 info.biCompression = 0;
183 #ifdef CONVERSE_ORDER
184 info.biCompression = change_unsigned_int_order(0);
185 #endif
186
187 info.biPlanes = 1;
188 #ifdef CONVERSE_ORDER
189 info.biPlanes = change_short_int_order(1);
190 #endif
191
192 info.biXPelsPerMeter = 3780;
193 #ifdef CONVERSE_ORDER
194 info.biXPelsPerMeter = change_int_order(3780);
195 #endif
196
197 info.biYPelsPerMeter = 3780;
198 #ifdef CONVERSE_ORDER
199 info.biYPelsPerMeter = change_int_order(3780);
200 #endif
201
202 putc('B', fp);
203 putc('M', fp);
204 fwrite(&(header.bfSize), sizeof(unsigned int), 1, fp);
205 fwrite(&header.bfReserved1, sizeof(unsigned short int), 1, fp);
206 fwrite(&header.bfReserved2, sizeof(unsigned short int), 1, fp);
207 fwrite(&header.bfOffBits, sizeof(unsigned int), 1, fp);
208 fwrite(&info, 40, 1, fp);
209 for (j = 0; j < 530; j++)
210 for (i = 400 - 1; i >= 0; i--) {
211 ri = (int)(graph[(j * 400 + i) * 3] * 256);
212 gi = (int)(graph[(j * 400 + i) * 3 + 1] * 256);
213 bi = (int)(graph[(j * 400 + i) * 3 + 2] * 256);
214 if (ri < 0) ri = 0;
215 if (ri > 255) ri = 255;
216 if (gi < 0) gi = 0;
217 if (gi > 255) gi = 255;
218 if (bi < 0) bi = 0;
219 if (bi > 255) bi = 255;
220 r = ri;
221 g = gi;
222 b = bi;
223 putc(b, fp);
224 putc(g, fp);
225 putc(r, fp);
226 }
227
228 fclose(fp);
229 HECMW_free(graph);
230 }
231 return;
232 }
233
generate_interval_point_vr(double tmincolor,double tmaxcolor,double * var,int * empty_flag,int nx,int ny,int nz,int mynode,int pesize,HECMW_Comm VIS_COMM,double * interval_point)234 void generate_interval_point_vr(double tmincolor, double tmaxcolor, double *var,
235 int *empty_flag, int nx, int ny, int nz,
236 int mynode, int pesize, HECMW_Comm VIS_COMM,
237 double *interval_point) {
238 int i, j;
239 double delta;
240 int count[500], t_count[500], tmp_count[500], sum_count, interv_count, sum,
241 current_j;
242 delta = (tmaxcolor - tmincolor) / 500.0;
243 for (i = 0; i < 500; i++) {
244 count[i] = 0;
245 t_count[i] = 0;
246 }
247 for (i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) {
248 if (empty_flag[i] == 1) {
249 j = (int)((var[i] - tmincolor) / delta);
250 if (j < 0) j = 0;
251 if (j > 499) j = 499;
252 count[j]++;
253 }
254 }
255 if (pesize > 1)
256 HECMW_Allreduce(count, t_count, 500, HECMW_INT, HECMW_SUM, VIS_COMM);
257 else {
258 for (i = 0; i < 500; i++) t_count[i] = count[i];
259 }
260 sum_count = 0;
261 for (i = 0; i < 500; i++) {
262 sum_count += t_count[i];
263 tmp_count[i] = t_count[i];
264 }
265 interv_count = (int)sum_count / 10;
266 current_j = 0;
267 interval_point[0] = tmincolor;
268 interval_point[1] = 0.0;
269 for (i = 1; i < 10; i++) {
270 interval_point[i * 2 + 1] = (double)i / 10.0;
271 sum = 0;
272 j = current_j;
273 while ((j < 500) && (sum < interv_count)) {
274 sum += t_count[j];
275 j++;
276 }
277 j--;
278 interval_point[i * 2] =
279 ((double)j / 500.0 +
280 1.0 / 500.0 *
281 (1.0 - (double)(sum - interv_count) / (double)tmp_count[j])) *
282 (tmaxcolor - tmincolor) +
283 tmincolor;
284 t_count[j] = sum - interv_count;
285 current_j = j;
286 }
287 interval_point[20] = tmaxcolor;
288 interval_point[21] = 1.0;
289 fprintf(stderr, "The automatic color mapping set is :\n");
290 for (i = 0; i < 11; i++)
291 fprintf(stderr, "%lf %lf \n", interval_point[i * 2],
292 interval_point[i * 2 + 1]);
293 return;
294 }
295
output_histogram_vr(double tmincolor,double tmaxcolor,double * var,int * empty_flag,int nx,int ny,int nz,int mynode,int pesize,HECMW_Comm VIS_COMM)296 void output_histogram_vr(double tmincolor, double tmaxcolor, double *var,
297 int *empty_flag, int nx, int ny, int nz, int mynode,
298 int pesize, HECMW_Comm VIS_COMM) {
299 int i, j;
300 double delta;
301 int count[100], t_count[100];
302 FILE *fp;
303
304 delta = (tmaxcolor - tmincolor) / 100.0;
305 for (i = 0; i < 100; i++) {
306 count[i] = 0;
307 t_count[i] = 0;
308 }
309 for (i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) {
310 if (empty_flag[i] == 1) {
311 j = (int)((var[i] - tmincolor) / delta);
312 if (j < 0) j = 0;
313 if (j > 99) j = 99;
314 count[j]++;
315 }
316 }
317 if (pesize > 1)
318 HECMW_Allreduce(count, t_count, 100, HECMW_INT, HECMW_SUM, VIS_COMM);
319 else {
320 for (i = 0; i < 100; i++) t_count[i] = count[i];
321 }
322 if (mynode == 0) {
323 fp = fopen("histogram.file", "w");
324 if (fp == NULL)
325 HECMW_vis_print_exit("Cannot generate the histogram output file");
326 for (i = 0; i < 100; i++)
327 fprintf(fp, "%d %d -----(%lf --- %lf)\n", i, t_count[i],
328 tmincolor + i * delta, tmincolor + (i + 1) * delta);
329 fclose(fp);
330 }
331 return;
332 }
333
find_minmax_vr(double * voxel_dxyz,double * voxel_orig_xyz,int mynode,double range[6])334 void find_minmax_vr(double *voxel_dxyz, double *voxel_orig_xyz, int mynode,
335 double range[6]) {
336 int i;
337 for (i = 0; i < 6; i++)
338 range[i] = voxel_orig_xyz[mynode * 3 + i / 2] +
339 (i % 2) * voxel_dxyz[mynode * 3 + i / 2];
340 return;
341 }
342
find_dis_minmax(double view_point_d[3],double vertex[24],double dis_minmax[2])343 void find_dis_minmax(double view_point_d[3], double vertex[24],
344 double dis_minmax[2]) {
345 int i;
346 double dis;
347 dis_minmax[0] = dis_minmax[1] =
348 sqrt(SQR(vertex[0]) + SQR(vertex[1]) + SQR(vertex[2]));
349 for (i = 0; i < 8; i++) {
350 dis = sqrt(SQR(vertex[i * 3] - view_point_d[0]) +
351 SQR(vertex[i * 3 + 1] - view_point_d[1]) +
352 SQR(vertex[i * 3 + 2] - view_point_d[2]));
353 if (dis < dis_minmax[0]) dis_minmax[0] = dis;
354 if (dis > dis_minmax[1]) dis_minmax[1] = dis;
355 }
356 return;
357 }
358
find_feap_minmax(int num_of_features,double * fea_point,double mincolor,double maxcolor,double feap_minmax[2])359 void find_feap_minmax(int num_of_features, double *fea_point, double mincolor,
360 double maxcolor, double feap_minmax[2]) {
361 int i, j;
362 double t, mint, color;
363
364 for (i = 0; i <= 255; i++) {
365 color = mincolor + (maxcolor - mincolor) * i / 255.0;
366 mint = 1.0E17;
367 for (j = 0; j < num_of_features; j++) {
368 t = fabs(color - fea_point[j]);
369 if (t < mint) mint = t;
370 }
371 if (mint > feap_minmax[1]) feap_minmax[1] = mint;
372 }
373 feap_minmax[0] = 0.0;
374 return;
375 }
376
find_feai_minmax(int num_of_features,double * fea_point,double mincolor,double maxcolor,double feai_minmax[2])377 void find_feai_minmax(int num_of_features, double *fea_point, double mincolor,
378 double maxcolor, double feai_minmax[2]) {
379 int i, j;
380 double t, t1, t2, mint, color;
381
382 for (i = 0; i <= 255; i++) {
383 color = mincolor + (maxcolor - mincolor) * i / 255.0;
384 mint = 1.0E17;
385 for (j = 0; j < num_of_features; j++) {
386 if ((t >= fea_point[j * 2]) && (t <= fea_point[j * 2 + 1]))
387 t = 0.0;
388 else {
389 t1 = fabs(color - fea_point[j * 2]);
390 t2 = fabs(color - fea_point[j * 2 + 1]);
391 if (t1 < t2)
392 t = t1;
393 else
394 t = t2;
395 }
396 if (t < mint) mint = t;
397 }
398 if (mint > feai_minmax[1]) feai_minmax[1] = mint;
399 }
400 feai_minmax[0] = 0.0;
401 return;
402 }
403