Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

gvl_calc.c

Last change on this file was 73852, checked in by neteler, 6 years ago
bulk fixing of typos in comments (using tools/fix_typos.sh)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id` Property svn:mime-type set to `text/x-csrc`
File size: 23.3 KB

Line
1	/*!
2	\file lib/ogsf/gvl_calc.c
3
4	\brief OGSF library - loading and manipulating volumes (lower level functions)
5
6	GRASS OpenGL gsurf OGSF Library
7
8	(C) 1999-2008 by the GRASS Development Team
9
10	This program is free software under the
11	GNU General Public License (>=v2).
12	Read the file COPYING that comes with GRASS
13	for details.
14
15	\author Tomas Paudits (February 2004)
16	\author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
17	*/
18
19	#include <math.h>
20
21	#include <grass/gis.h>
22	#include <grass/ogsf.h>
23
24	#include "rgbpack.h"
25	#include "mc33_table.h"
26
27	/*!
28	\brief memory buffer for writing
29	*/
30	#define BUFFER_SIZE 1000000
31
32	/* USEFUL MACROS */
33
34	/* interp. */
35	#define LINTERP(d,a,b) (a + d * (b - a))
36	#define TINTERP(d,v) ((v[0](1.-d[0])(1.-d[1])*(1.-d[2])) +\
37	(v[1]d[0](1.-d[1])*(1.-d[2])) + \
38	(v[2]d[0]d[1]*(1.-d[2])) + \
39	(v[3](1.-d[0])d[1]*(1.-d[2])) + \
40	(v[4](1.-d[0])(1.-d[1])*d[2]) + \
41	(v[5]d[0](1.-d[1])*d[2]) + \
42	(v[6]d[0]d[1]*d[2]) + \
43	(v[7](1.-d[0])d[1]*d[2]))
44
45	#define FOR_VAR i_for
46	#define FOR_0_TO_N(n, cmd) { int FOR_VAR; for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) {cmd;} }
47
48	/*!
49	\brief writing and reading isosurface data
50	*/
51	#define WRITE(c) gvl_write_char(dbuff->ndx_new++, &(dbuff->new), c)
52	#define READ() gvl_read_char(dbuff->ndx_old++, dbuff->old)
53	#define SKIP(n) dbuff->ndx_old = dbuff->ndx_old + n
54
55	/*!
56	\brief check and set data descriptor
57	*/
58	#define IS_IN_DATA(att) ((isosurf->data_desc >> att) & 1)
59	#define SET_IN_DATA(att) isosurf->data_desc = (isosurf->data_desc \| (1 << att))
60
61	typedef struct
62	{
63	unsigned char *old;
64	unsigned char *new;
65	int ndx_old;
66	int ndx_new;
67	int num_zero;
68	} data_buffer;
69
70	int mc33_process_cube(int c_ndx, float *v);
71
72	/* global variables */
73	int Rows, Cols, Depths;
74	double ResX, ResY, ResZ;
75
76	/************************************************************************/
77	/* ISOSURFACES */
78
79	/*!
80	\brief Write cube index
81
82	\param ndx
83	\param dbuff
84	*/
85	void iso_w_cndx(int ndx, data_buffer * dbuff)
86	{
87	/* cube don't contains polys */
88	if (ndx == -1) {
89	if (dbuff->num_zero == 0) {
90	WRITE(0);
91	dbuff->num_zero++;
92	}
93	else if (dbuff->num_zero == 254) {
94	WRITE(dbuff->num_zero + 1);
95	dbuff->num_zero = 0;
96	}
97	else {
98	dbuff->num_zero++;
99	}
100	}
101	else { /* isosurface cube */
102	if (dbuff->num_zero == 0) {
103	WRITE((ndx / 256) + 1);
104	WRITE(ndx % 256);
105	}
106	else {
107	WRITE(dbuff->num_zero);
108	dbuff->num_zero = 0;
109	WRITE((ndx / 256) + 1);
110	WRITE(ndx % 256);
111	}
112	}
113	}
114
115	/*!
116	\brief Read cube index
117
118	\param dbuff
119	*/
120	int iso_r_cndx(data_buffer * dbuff)
121	{
122	int ndx, ndx2;
123
124	if (dbuff->num_zero != 0) {
125	dbuff->num_zero--;
126	ndx = -1;
127	}
128	else {
129	WRITE(ndx = READ());
130	if (ndx == 0) {
131	WRITE(dbuff->num_zero = READ());
132	dbuff->num_zero--;
133	ndx = -1;
134	}
135	else {
136	WRITE(ndx2 = READ());
137	ndx = (ndx - 1) * 256 + ndx2;
138	}
139	}
140
141	return ndx;
142	}
143
144	/*!
145	\brief Get value from data input
146
147	\param isosurf
148	\param desc
149	\param x,y,z
150	\param[out] value
151
152	\return 0
153	\return ?
154	*/
155	int iso_get_cube_value(geovol_isosurf * isosurf, int desc, int x, int y,
156	int z, float *v)
157	{
158	double d;
159	geovol_file *vf;
160	int type, ret = 1;
161
162	/* get volume file from attribute handle */
163	vf = gvl_file_get_volfile(isosurf->att[desc].hfile);
164	type = gvl_file_get_data_type(vf);
165
166	/* get value from volume file */
167	if (type == VOL_DTYPE_FLOAT) {
168	gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
169	(int)(z * ResZ), v);
170	}
171	else if (type == VOL_DTYPE_DOUBLE) {
172	gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
173	(int)(z * ResZ), &d);
174	*v = (float)d;
175	}
176	else {
177	return 0;
178	}
179
180	/* null check */
181	if (gvl_file_is_null_value(vf, v))
182	ret = 0;
183
184	/* adjust data */
185	switch (desc) {
186	case (ATT_TOPO):
187	v = (v) - isosurf->att[desc].constant;
188	break;
189	case (ATT_MASK):
190	if (isosurf->att[desc].constant)
191	ret = !ret;
192	break;
193	}
194
195	return ret;
196	}
197
198	/*!
199	\brief Get volume file values range
200
201	\param isosurf
202	\param desc
203	\param[out] min
204	\param[out] max
205	*/
206	void iso_get_range(geovol_isosurf * isosurf, int desc, double *min,
207	double *max)
208	{
209	gvl_file_get_min_max(gvl_file_get_volfile(isosurf->att[desc].hfile), min,
210	max);
211	}
212
213	/*!
214	\brief Read values for cube
215
216	\param isosurf
217	\param desc
218	\param x,y,z
219	\param[out] v
220
221	\return
222	*/
223	int iso_get_cube_values(geovol_isosurf * isosurf, int desc, int x, int y,
224	int z, float *v)
225	{
226	int p, ret = 1;
227
228	for (p = 0; p < 8; ++p) {
229	if (iso_get_cube_value
230	(isosurf, desc, x + ((p ^ (p >> 1)) & 1), y + ((p >> 1) & 1),
231	z + ((p >> 2) & 1), &v[p]) == 0) {
232	ret = 0;
233	}
234	}
235
236	return ret;
237	}
238
239	/*!
240	\brief Calculate cube grads
241
242	\param isosurf
243	\param x,y,z
244	\param grad
245	*/
246	void iso_get_cube_grads(geovol_isosurf * isosurf, int x, int y, int z,
247	float (*grad)[3])
248	{
249	float v[3];
250	int i, j, k, p;
251
252	for (p = 0; p < 8; ++p) {
253	i = x + ((p ^ (p >> 1)) & 1);
254	j = y + ((p >> 1) & 1);
255	k = z + ((p >> 2) & 1);
256
257	/* x */
258	if (i == 0) {
259	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
260	iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
261	grad[p][0] = v[2] - v[1];
262	}
263	else {
264	if (i == (Cols - 1)) {
265	iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
266	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
267	grad[p][0] = v[1] - v[0];
268	}
269	else {
270	iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
271	iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
272	grad[p][0] = (v[2] - v[0]) / 2;
273	}
274	}
275
276	/* y */
277	if (j == 0) {
278	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
279	iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
280	grad[p][1] = v[2] - v[1];
281	}
282	else {
283	if (j == (Rows - 1)) {
284	iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
285	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
286	grad[p][1] = v[1] - v[0];
287	}
288	else {
289	iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
290	iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
291	grad[p][1] = (v[2] - v[0]) / 2;
292	}
293	}
294
295	/* z */
296	if (k == 0) {
297	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
298	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
299	grad[p][2] = v[2] - v[1];
300	}
301	else {
302	if (k == (Depths - 1)) {
303	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
304	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
305	grad[p][2] = v[1] - v[0];
306	}
307	else {
308	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
309	iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
310	grad[p][2] = (v[2] - v[0]) / 2;
311	}
312	}
313	}
314	}
315
316	/*!
317	\brief Process cube
318
319	\param isosurf
320	\param x,y,z
321	\param dbuff
322	*/
323	void iso_calc_cube(geovol_isosurf * isosurf, int x, int y, int z,
324	data_buffer * dbuff)
325	{
326	int i, c_ndx;
327	int crnt, v1, v2, c;
328	float val[MAX_ATTS][8], grad[8][3];
329	float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
330	double min, max;
331
332	if (isosurf->att[ATT_TOPO].changed) {
333	/* read topo values, if there are NULL values then return */
334	if (!iso_get_cube_values(isosurf, ATT_TOPO, x, y, z, val[ATT_TOPO])) {
335	iso_w_cndx(-1, dbuff);
336	return;
337	}
338
339	/* mask */
340	if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
341	if (!iso_get_cube_values
342	(isosurf, ATT_MASK, x, y, z, val[ATT_MASK])) {
343	iso_w_cndx(-1, dbuff);
344	return;
345	}
346	}
347
348	/* index to precalculated table */
349	c_ndx = 0;
350	for (i = 0; i < 8; i++) {
351	if (val[ATT_TOPO][i] > 0)
352	c_ndx \|= 1 << i;
353	}
354	c_ndx = mc33_process_cube(c_ndx, val[ATT_TOPO]);
355
356	iso_w_cndx(c_ndx, dbuff);
357
358	if (c_ndx == -1)
359	return;
360
361	/* calc cube grads */
362	iso_get_cube_grads(isosurf, x, y, z, grad);
363
364	}
365	else {
366	/* read cube index */
367	if ((c_ndx = iso_r_cndx(dbuff)) == -1)
368	return;
369	}
370
371	/* get color values */
372	if (isosurf->att[ATT_COLOR].changed &&
373	isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
374	iso_get_cube_values(isosurf, ATT_COLOR, x, y, z, val[ATT_COLOR]);
375	}
376
377	/* get transparency values */
378	if (isosurf->att[ATT_TRANSP].changed &&
379	isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
380	iso_get_cube_values(isosurf, ATT_TRANSP, x, y, z, val[ATT_TRANSP]);
381	}
382
383	/* get shine values */
384	if (isosurf->att[ATT_SHINE].changed &&
385	isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
386	iso_get_cube_values(isosurf, ATT_SHINE, x, y, z, val[ATT_SHINE]);
387	}
388
389	/* get emit values */
390	if (isosurf->att[ATT_EMIT].changed &&
391	isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
392	iso_get_cube_values(isosurf, ATT_EMIT, x, y, z, val[ATT_EMIT]);
393	}
394
395	FOR_0_TO_N(3, d_sum[FOR_VAR] = 0.;
396	n_sum[FOR_VAR] = 0.);
397
398	/* loop in edges */
399	for (i = 0; i < cell_table[c_ndx].nedges; i++) {
400	/* get edge number */
401	crnt = cell_table[c_ndx].edges[i];
402
403	/* set topo */
404	if (isosurf->att[ATT_TOPO].changed) {
405	/* interior vertex */
406	if (crnt == 12) {
407	FOR_0_TO_N(3,
408	WRITE((d3[FOR_VAR] =
409	d_sum[FOR_VAR] /
410	((float)(cell_table[c_ndx].nedges))) *
411	255));
412	GS_v3norm(n_sum);
413	FOR_0_TO_N(3,
414	WRITE((n_sum[FOR_VAR] /
415	((float)(cell_table[c_ndx].nedges)) +
416	1.) * 127));
417	/* edge vertex */
418	}
419	else {
420	/* set egdes verts */
421	v1 = edge_vert[crnt][0];
422	v2 = edge_vert[crnt][1];
423
424	/* calc intersection point - edge and isosurf */
425	d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] -
426	val[ATT_TOPO][v2]);
427
428	d_sum[edge_vert_pos[crnt][0]] += d;
429	d_sum[edge_vert_pos[crnt][1]] += edge_vert_pos[crnt][2];
430	d_sum[edge_vert_pos[crnt][3]] += edge_vert_pos[crnt][4];
431
432	WRITE(d * 255);
433
434	/* set normal for intersect. point */
435	FOR_0_TO_N(3, n[FOR_VAR] =
436	LINTERP(d, grad[v1][FOR_VAR], grad[v2][FOR_VAR]));
437	GS_v3norm(n);
438	FOR_0_TO_N(3, n_sum[FOR_VAR] += n[FOR_VAR]);
439	FOR_0_TO_N(3, WRITE((n[FOR_VAR] + 1.) * 127));
440	}
441	}
442	else {
443	/* read x,y,z of intersection point in cube coords */
444	if (crnt == 12) {
445	WRITE(c = READ());
446	d3[0] = ((float)c) / 255.0;
447	WRITE(c = READ());
448	d3[1] = ((float)c) / 255.0;
449	WRITE(c = READ());
450	d3[2] = ((float)c) / 255.0;
451	}
452	else {
453	/* set egdes verts */
454	v1 = edge_vert[crnt][0];
455	v2 = edge_vert[crnt][1];
456
457	WRITE(c = READ());
458	d = ((float)c) / 255.0;
459	}
460
461	/* set normals */
462	FOR_0_TO_N(3, WRITE(READ()));
463	}
464
465	/* set color */
466	if (isosurf->att[ATT_COLOR].changed &&
467	isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
468	if (crnt == 12) {
469	tv = TINTERP(d3, val[ATT_COLOR]);
470	}
471	else {
472	tv = LINTERP(d, val[ATT_COLOR][v1], val[ATT_COLOR][v2]);
473	}
474
475	c = Gvl_get_color_for_value(isosurf->att[ATT_COLOR].att_data,
476	&tv);
477
478	WRITE(c & RED_MASK);
479	WRITE((c & GRN_MASK) >> 8);
480	WRITE((c & BLU_MASK) >> 16);
481
482	if (IS_IN_DATA(ATT_COLOR))
483	SKIP(3);
484	}
485	else {
486	if (isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
487	FOR_0_TO_N(3, WRITE(READ()));
488	}
489	else {
490	if (IS_IN_DATA(ATT_COLOR))
491	SKIP(3);
492	}
493	}
494
495	/* set transparency */
496	if (isosurf->att[ATT_TRANSP].changed &&
497	isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
498	if (crnt == 12) {
499	tv = TINTERP(d3, val[ATT_TRANSP]);
500	}
501	else {
502	tv = LINTERP(d, val[ATT_TRANSP][v1], val[ATT_TRANSP][v2]);
503	}
504
505	iso_get_range(isosurf, ATT_TRANSP, &min, &max);
506	c = (min != max) ? 255 - (tv - min) / (max - min) * 255 : 0;
507
508	WRITE(c);
509	if (IS_IN_DATA(ATT_TRANSP))
510	SKIP(1);
511	}
512	else {
513	if (isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
514	WRITE(READ());
515	}
516	else {
517	if (IS_IN_DATA(ATT_TRANSP))
518	SKIP(1);
519	}
520	}
521
522	/* set shin */
523	if (isosurf->att[ATT_SHINE].changed &&
524	isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
525	if (crnt == 12) {
526	tv = TINTERP(d3, val[ATT_SHINE]);
527	}
528	else {
529	tv = LINTERP(d, val[ATT_SHINE][v1], val[ATT_SHINE][v2]);
530	}
531
532	iso_get_range(isosurf, ATT_SHINE, &min, &max);
533	c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
534
535	WRITE(c);
536	if (IS_IN_DATA(ATT_SHINE))
537	SKIP(1);
538	}
539	else {
540	if (isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
541	WRITE(READ());
542	}
543	else {
544	if (IS_IN_DATA(ATT_SHINE))
545	SKIP(1);
546	}
547	}
548
549	/* set emit */
550	if (isosurf->att[ATT_EMIT].changed &&
551	isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
552	if (crnt == 12) {
553	tv = TINTERP(d3, val[ATT_EMIT]);
554	}
555	else {
556	tv = LINTERP(d, val[ATT_EMIT][v1], val[ATT_EMIT][v2]);
557	}
558
559	iso_get_range(isosurf, ATT_EMIT, &min, &max);
560	c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
561
562	WRITE(c);
563	if (IS_IN_DATA(ATT_SHINE))
564	SKIP(1);
565	}
566	else {
567	if (isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
568	WRITE(READ());
569	}
570	else {
571	if (IS_IN_DATA(ATT_EMIT))
572	SKIP(1);
573	}
574	}
575	}
576	}
577
578	/*!
579	\brief Fill data structure with computed isosurfaces polygons
580
581	\param gvol pointer to geovol struct
582
583	\return 1
584	*/
585	int gvl_isosurf_calc(geovol * gvol)
586	{
587	int x, y, z;
588	int i, a, read;
589	geovol_file *vf;
590	geovol_isosurf *isosurf;
591
592	data_buffer *dbuff;
593	int *need_update, need_update_global;
594
595	dbuff = G_malloc(gvol->n_isosurfs * sizeof(data_buffer));
596	need_update = G_malloc(gvol->n_isosurfs * sizeof(int));
597
598	/* flag - changed any isosurface */
599	need_update_global = 0;
600
601	/* initialize */
602	for (i = 0; i < gvol->n_isosurfs; i++) {
603	isosurf = gvol->isosurf[i];
604
605	/* initialize read/write buffers */
606	dbuff[i].old = NULL;
607	dbuff[i].new = NULL;
608	dbuff[i].ndx_old = 0;
609	dbuff[i].ndx_new = 0;
610	dbuff[i].num_zero = 0;
611
612	need_update[i] = 0;
613	for (a = 1; a < MAX_ATTS; a++) {
614	if (isosurf->att[a].changed) {
615	read = 0;
616	/* changed to map attribute */
617	if (isosurf->att[a].att_src == MAP_ATT) {
618	vf = gvl_file_get_volfile(isosurf->att[a].hfile);
619	read = 1;
620	}
621	/* changed threshold value */
622	if (a == ATT_TOPO) {
623	isosurf->att[a].hfile = gvol->hfile;
624	vf = gvl_file_get_volfile(gvol->hfile);
625	read = 1;
626	}
627	/* initialize reading in selected mode */
628	if (read) {
629	gvl_file_set_mode(vf, 3);
630	gvl_file_start_read(vf);
631	}
632
633	/* set update flag - isosurface will be calc */
634	if (read \|\| IS_IN_DATA(a)) {
635	need_update[i] = 1;
636	need_update_global = 1;
637	}
638	}
639	}
640
641	if (need_update[i]) {
642	/* set data buffer */
643	dbuff[i].old = isosurf->data;
644	}
645	}
646
647	/* calculate if only some isosurface changed */
648	if (need_update_global) {
649
650	ResX = gvol->isosurf_x_mod;
651	ResY = gvol->isosurf_y_mod;
652	ResZ = gvol->isosurf_z_mod;
653
654	Cols = gvol->cols / ResX;
655	Rows = gvol->rows / ResY;
656	Depths = gvol->depths / ResZ;
657
658	/* calc isosurface - marching cubes - start */
659
660	for (z = 0; z < Depths - 1; z++) {
661	for (y = 0; y < Rows - 1; y++) {
662	for (x = 0; x < Cols - 1; x++) {
663	for (i = 0; i < gvol->n_isosurfs; i++) {
664	/* recalculate only changed isosurfaces */
665	if (need_update[i]) {
666	iso_calc_cube(gvol->isosurf[i], x, y, z,
667	&dbuff[i]);
668	}
669	}
670	}
671	}
672	}
673
674	}
675	/* end */
676
677	/* deinitialize */
678	for (i = 0; i < gvol->n_isosurfs; i++) {
679	isosurf = gvol->isosurf[i];
680
681	/* set new isosurface data */
682	if (need_update[i]) {
683	if (dbuff[i].num_zero != 0)
684	gvl_write_char(dbuff[i].ndx_new++, &(dbuff[i].new),
685	dbuff[i].num_zero);
686
687	if (dbuff[i].old == isosurf->data)
688	dbuff[i].old = NULL;
689	G_free(isosurf->data);
690	gvl_align_data(dbuff[i].ndx_new, &(dbuff[i].new));
691	isosurf->data = dbuff[i].new;
692	isosurf->data_desc = 0;
693	}
694
695	for (a = 1; a < MAX_ATTS; a++) {
696	if (isosurf->att[a].changed) {
697	read = 0;
698	/* changed map attribute */
699	if (isosurf->att[a].att_src == MAP_ATT) {
700	vf = gvl_file_get_volfile(isosurf->att[a].hfile);
701	read = 1;
702	}
703	/* changed threshold value */
704	if (a == ATT_TOPO) {
705	isosurf->att[a].hfile = gvol->hfile;
706	vf = gvl_file_get_volfile(gvol->hfile);
707	read = 1;
708	}
709	/* deinitialize reading */
710	if (read) {
711	gvl_file_end_read(vf);
712
713	/* set data description */
714	SET_IN_DATA(a);
715	}
716	isosurf->att[a].changed = 0;
717	}
718	else if (isosurf->att[a].att_src == MAP_ATT) {
719	/* set data description */
720	SET_IN_DATA(a);
721	}
722	}
723	}
724
725	/* TODO: G_free() dbuff and need_update ??? */
726
727	return (1);
728	}
729
730	/*!
731	\brief ADD
732
733	\param pos
734	\param data
735	\param c
736	*/
737	void gvl_write_char(int pos, unsigned char **data, unsigned char c)
738	{
739	/* check to need allocation memory */
740	if ((pos % BUFFER_SIZE) == 0) {
741	data = (char )G_realloc(*data,
742	sizeof(char) * ((pos / BUFFER_SIZE) +
743	1) * BUFFER_SIZE);
744	if (!(*data)) {
745	return;
746	}
747
748	G_debug(3,
749	"gvl_write_char(): reallocate memory for pos : %d to : %lu B",
750	pos, sizeof(char) * ((pos / BUFFER_SIZE) + 1) * BUFFER_SIZE);
751	}
752
753	(*data)[pos] = c;
754	}
755
756	/*!
757	\brief Read char
758
759	\param pos position index
760	\param data data buffer
761
762	\return char on success
763	\return NULL on failure
764	*/
765	unsigned char gvl_read_char(int pos, const unsigned char *data)
766	{
767	if (!data)
768	return '\0';
769
770	return data[pos];
771	}
772
773	/*!
774	\brief Append data to buffer
775
776	\param pos position index
777	\param data data buffer
778	*/
779	void gvl_align_data(int pos, unsigned char **data)
780	{
781	unsigned char p = data;
782
783
784	/* realloc memory to fit in data length */
785	p = (unsigned char )G_realloc(p, sizeof(unsigned char) pos); /* G_fatal_error */
786	if (!p) {
787	return;
788	}
789
790	G_debug(3, "gvl_align_data(): reallocate memory finally to : %d B", pos);
791
792	if (pos == 0)
793	p = NULL;
794
795	*data = p;
796
797	return;
798	}
799
800	/************************************************************************/
801	/* SLICES */
802
803	/************************************************************************/
804
805	#define DISTANCE_2(x1, y1, x2, y2) sqrt((x1-x2)(x1-x2)+(y1-y2)(y1-y2))
806
807	#define SLICE_MODE_INTERP_NO 0
808	#define SLICE_MODE_INTERP_YES 1
809
810	/*!
811	\brief Get volume value
812
813	\param gvl pointer to geovol struct
814	\param x,y,z
815
816	\return value
817	*/
818	float slice_get_value(geovol * gvl, int x, int y, int z)
819	{
820	static double d;
821	static geovol_file *vf;
822	static int type;
823	static float value;
824
825	if (x < 0 \|\| y < 0 \|\| z < 0 \|\| (x > gvl->cols - 1) \|\| (y > gvl->rows - 1)
826	\|\| (z > gvl->depths - 1))
827	return 0.;
828
829	/* get volume file from attribute handle */
830	vf = gvl_file_get_volfile(gvl->hfile);
831	type = gvl_file_get_data_type(vf);
832
833	/* get value from volume file */
834	if (type == VOL_DTYPE_FLOAT) {
835	gvl_file_get_value(vf, x, y, z, &value);
836	}
837	else if (type == VOL_DTYPE_DOUBLE) {
838	gvl_file_get_value(vf, x, y, z, &d);
839	value = (float)d;
840	}
841	else {
842	return 0.;
843	}
844
845	return value;
846	}
847
848	/*!
849	\brief Calculate slices
850
851	\param gvl pointer to geovol struct
852	\param ndx_slc
853	\param colors
854
855	\return 1
856	*/
857	int slice_calc(geovol * gvl, int ndx_slc, void *colors)
858	{
859	int cols, rows, c, r;
860	int i, j, k, pos, color;
861	int p_x, p_y, *p_z;
862	float p_ex, p_ey, *p_ez;
863	float value, v[8];
864	float x, y, z, ei, ej, ek, stepx, stepy, stepz;
865	float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
866
867	geovol_slice *slice;
868	geovol_file *vf;
869
870	slice = gvl->slice[ndx_slc];
871
872	/* set mods, pointer to x, y, z step value */
873	if (slice->dir == X) {
874	modx = ResY;
875	mody = ResZ;
876	modz = ResX;
877	p_x = &k;
878	p_y = &i;
879	p_z = &j;
880	p_ex = &ek;
881	p_ey = &ei;
882	p_ez = &ej;
883	}
884	else if (slice->dir == Y) {
885	modx = ResX;
886	mody = ResZ;
887	modz = ResY;
888	p_x = &i;
889	p_y = &k;
890	p_z = &j;
891	p_ex = &ei;
892	p_ey = &ek;
893	p_ez = &ej;
894	}
895	else {
896	modx = ResX;
897	mody = ResY;
898	modz = ResZ;
899	p_x = &i;
900	p_y = &j;
901	p_z = &k;
902	p_ex = &ei;
903	p_ey = &ej;
904	p_ez = &ek;
905	}
906
907	/* distance between slice def. points */
908	distxy = DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
909	distz = fabsf(slice->z2 - slice->z1);
910
911	/* distance between slice def points is zero - nothing to do */
912	if (distxy == 0. \|\| distz == 0.) {
913	return (1);
914	}
915
916	/* start reading volume file */
917	vf = gvl_file_get_volfile(gvl->hfile);
918	gvl_file_set_mode(vf, 3);
919	gvl_file_start_read(vf);
920
921	/* set xy resolution */
922	modxy =
923	DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
924	(slice->y2 - slice->y1) / distxy * mody, 0., 0.);
925
926	/* cols/rows of slice */
927	f_cols = distxy / modxy;
928	cols = f_cols > (int)f_cols ? (int)f_cols + 1 : (int)f_cols;
929
930	f_rows = distz / modz;
931	rows = f_rows > (int)f_rows ? (int)f_rows + 1 : (int)f_rows;
932
933	/* set x,y initially to first slice point */
934	x = slice->x1;
935	y = slice->y1;
936
937	/* set x,y step */
938	stepx = (slice->x2 - slice->x1) / f_cols;
939	stepy = (slice->y2 - slice->y1) / f_cols;
940	stepz = (slice->z2 - slice->z1) / f_rows;
941
942	/* set position in slice data */
943	pos = 0;
944
945	/* loop in slice cols */
946	for (c = 0; c < cols + 1; c++) {
947
948	/* convert x, y to integer - index in grid */
949	i = (int)x;
950	j = (int)y;
951
952	/* distance between index and real position */
953	ei = x - (float)i;
954	ej = y - (float)j;
955
956	/* set z to slice z1 point */
957	z = slice->z1;
958
959	/* loop in slice rows */
960	for (r = 0; r < rows + 1; r++) {
961
962	/* distance between index and real position */
963	k = (int)z;
964	ek = z - (float)k;
965
966	/* get interpolated value */
967	if (slice->mode == SLICE_MODE_INTERP_YES) {
968	/* get grid values */
969	v[0] = slice_get_value(gvl, p_x, p_y, *p_z);
970	v[1] = slice_get_value(gvl, p_x + 1, p_y, *p_z);
971	v[2] = slice_get_value(gvl, p_x, p_y + 1, *p_z);
972	v[3] = slice_get_value(gvl, p_x + 1, p_y + 1, *p_z);
973
974	v[4] = slice_get_value(gvl, p_x, p_y, *p_z + 1);
975	v[5] = slice_get_value(gvl, p_x + 1, p_y, *p_z + 1);
976	v[6] = slice_get_value(gvl, p_x, p_y + 1, *p_z + 1);
977	v[7] = slice_get_value(gvl, p_x + 1, p_y + 1, *p_z + 1);
978
979	/* get interpolated value */
980	value = v[0] * (1. - p_ex) (1. - p_ey) (1. - *p_ez)
981	+ v[1] * (p_ex) (1. - p_ey) (1. - *p_ez)
982	+ v[2] * (1. - p_ex) (p_ey) (1. - *p_ez)
983	+ v[3] * (p_ex) (p_ey) (1. - *p_ez)
984	+ v[4] * (1. - p_ex) (1. - p_ey) (*p_ez)
985	+ v[5] * (p_ex) (1. - p_ey) (*p_ez)
986	+ v[6] * (1. - p_ex) (p_ey) (*p_ez)
987	+ v[7] * (p_ex) (p_ey) (*p_ez);
988
989	/* no interp value */
990	}
991	else {
992	value = slice_get_value(gvl, p_x, p_y, *p_z);
993	}
994
995	/* translate value to color */
996	color = Gvl_get_color_for_value(colors, &value);
997
998	/* write color to slice data */
999	gvl_write_char(pos++, &(slice->data), color & RED_MASK);
1000	gvl_write_char(pos++, &(slice->data), (color & GRN_MASK) >> 8);
1001	gvl_write_char(pos++, &(slice->data), (color & BLU_MASK) >> 16);
1002
1003	/* step in z */
1004	if (r + 1 > f_rows) {
1005	z += stepz * (f_rows - (float)r);
1006	}
1007	else {
1008	z += stepz;
1009	}
1010	}
1011
1012	/* step in x,y */
1013	if (c + 1 > f_cols) {
1014	x += stepx * (f_cols - (float)c);
1015	y += stepy * (f_cols - (float)c);
1016	}
1017	else {
1018	x += stepx;
1019	y += stepy;
1020	}
1021	}
1022
1023	/* end reading volume file */
1024	gvl_file_end_read(vf);
1025	gvl_align_data(pos, &(slice->data));
1026
1027	return (1);
1028	}
1029
1030	/*!
1031	\brief Calculate slices for given volume set
1032
1033	\param gvol pointer to geovol struct
1034
1035	\return 1
1036	*/
1037	int gvl_slices_calc(geovol *gvol)
1038	{
1039	int i;
1040	void *colors;
1041
1042	G_debug(5, "gvl_slices_calc(): id=%d", gvol->gvol_id);
1043
1044	/* set current resolution */
1045	ResX = gvol->slice_x_mod;
1046	ResY = gvol->slice_y_mod;
1047	ResZ = gvol->slice_z_mod;
1048
1049	/* set current num of cols, rows, depths */
1050	Cols = gvol->cols / ResX;
1051	Rows = gvol->rows / ResY;
1052	Depths = gvol->depths / ResZ;
1053
1054	/* load colors for geovol file */
1055	Gvl_load_colors_data(&colors, gvl_file_get_name(gvol->hfile));
1056
1057	/* calc changed slices */
1058	for (i = 0; i < gvol->n_slices; i++) {
1059	if (gvol->slice[i]->changed) {
1060	slice_calc(gvol, i, colors);
1061
1062	/* set changed flag */
1063	gvol->slice[i]->changed = 0;
1064	}
1065	}
1066
1067	/* free color */
1068	Gvl_unload_colors_data(colors);
1069
1070	return (1);
1071	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: grass/trunk/lib/ogsf/gvl_calc.c

Download in other formats: