1
2 /****************************************************************************
3 *
4 * MODULE: r.out.vtk
5 *
6 * AUTHOR(S): Original author
7 * Soeren Gebbert soerengebbert@gmx.de
8 * 08 23 2005 Berlin
9 * PURPOSE: Converts raster maps into the VTK-Ascii format
10 *
11 * COPYRIGHT: (C) 2005 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <grass/gis.h>
23 #include <grass/raster.h>
24 #include <grass/glocale.h>
25 #include <grass/config.h>
26 #include "globaldefs.h"
27
28 /*local prototypes */
29 static double get_raster_value_as_double(int maptype, void *ptr,
30 double nullval);
31
32
33 /* ************************************************************************* */
34 /* Get the value of the current raster pointer as double ******************* */
35 /* ************************************************************************* */
get_raster_value_as_double(int MapType,void * ptr,double nullval)36 double get_raster_value_as_double(int MapType, void *ptr, double nullval)
37 {
38 double val = nullval;
39
40 if (MapType == CELL_TYPE) {
41 if (Rast_is_null_value(ptr, MapType)) {
42 val = nullval;
43 }
44 else {
45 val = *(CELL *) ptr;
46 }
47 }
48 if (MapType == FCELL_TYPE) {
49 if (Rast_is_null_value(ptr, MapType)) {
50 val = nullval;
51 }
52 else {
53 val = *(FCELL *) ptr;
54 }
55 }
56 if (MapType == DCELL_TYPE) {
57 if (Rast_is_null_value(ptr, MapType)) {
58 val = nullval;
59 }
60 else {
61 val = *(DCELL *) ptr;
62 }
63 }
64
65 return val;
66 }
67
68 /* ************************************************************************* */
69 /* Write the default VTK Header, Elevation is not supportet *************** */
70 /* ************************************************************************* */
71 void
write_vtk_normal_header(FILE * fp,struct Cell_head region,double elevation,int type)72 write_vtk_normal_header(FILE * fp, struct Cell_head region, double elevation,
73 int type)
74 {
75 G_debug(3, _("write_vtk_normal_header: Writing VTK-Header"));
76
77 /*Simple vtk ASCII header */
78 fprintf(fp, "# vtk DataFile Version 3.0\n");
79 fprintf(fp, "GRASS GIS 7 Export\n");
80 fprintf(fp, "ASCII\n");
81 fprintf(fp, "DATASET STRUCTURED_POINTS\n"); /*We are using the structured point dataset. */
82
83 if (type)
84 fprintf(fp, "DIMENSIONS %i %i %i\n", region.cols, region.rows, 1);
85 else
86 fprintf(fp, "DIMENSIONS %i %i %i\n", region.cols + 1, region.rows + 1,
87 1);
88
89 fprintf(fp, "SPACING %lf %lf %lf\n", region.ew_res, region.ns_res, 0.0);
90
91 if (type)
92 fprintf(fp, "ORIGIN %lf %lf %lf\n",
93 (region.west + region.ew_res / 2) - x_extent,
94 (region.south + region.ns_res / 2) - y_extent, elevation);
95 else
96 fprintf(fp, "ORIGIN %lf %lf %lf\n", region.west - x_extent,
97 region.south - y_extent, elevation);
98
99 }
100
101
102 /* ************************************************************************* */
103 /* Write the Elevation VTK Header, Elevation is supportet ****************** */
104 /* ************************************************************************* */
write_vtk_structured_elevation_header(FILE * fp,struct Cell_head region)105 void write_vtk_structured_elevation_header(FILE * fp, struct Cell_head region)
106 {
107 G_debug(3,
108 _("write_vtk_structured_elevation_header: Writing VTK-Header"));
109
110 /*Simple vtk ASCII header */
111 fprintf(fp, "# vtk DataFile Version 3.0\n");
112 fprintf(fp, "GRASS GIS 7 Export\n");
113 fprintf(fp, "ASCII\n");
114 fprintf(fp, "DATASET STRUCTURED_GRID\n"); /*We are using the structured grid dataset. */
115 fprintf(fp, "DIMENSIONS %i %i %i\n", region.cols, region.rows, 1);
116 fprintf(fp, "POINTS %i float\n", region.cols * region.rows); /*The Coordinates */
117 }
118
119 /* ************************************************************************* */
120 /* Write the Rectilinear Elevtaion VTK Header, Elevation is supportet ****** */
121 /* ************************************************************************* */
write_vtk_polygonal_elevation_header(FILE * fp,struct Cell_head region)122 void write_vtk_polygonal_elevation_header(FILE * fp, struct Cell_head region)
123 {
124 G_debug(3, _("write_vtk_polygonal_elevation_header: Writing VTK-Header"));
125
126 /*Simple vtk ASCII header */
127 fprintf(fp, "# vtk DataFile Version 3.0\n");
128 fprintf(fp, "GRASS GIS 7 Export\n");
129 fprintf(fp, "ASCII\n");
130 fprintf(fp, "DATASET POLYDATA\n"); /*We are using polydataset. */
131 fprintf(fp, "POINTS %i float\n", region.cols * region.rows);
132 }
133
134 /* ************************************************************************* */
135 /* Write the CellDataHeader ************************************************ */
136 /* ************************************************************************* */
write_vtk_celldata_header(FILE * fp,struct Cell_head region)137 void write_vtk_celldata_header(FILE * fp, struct Cell_head region)
138 {
139 G_debug(3, _("write_vtk_celldata_header: Writing VTK-DataHeader"));
140 fprintf(fp, "CELL_DATA %i\n", region.cols * region.rows);
141 }
142
143 /* ************************************************************************* */
144 /* Write the PointDataHeader ************************************************ */
145 /* ************************************************************************* */
write_vtk_pointdata_header(FILE * fp,struct Cell_head region)146 void write_vtk_pointdata_header(FILE * fp, struct Cell_head region)
147 {
148 G_debug(3, _("writeVTKPointHeader: Writing VTK-DataHeader"));
149 fprintf(fp, "POINT_DATA %i\n", region.cols * region.rows);
150 }
151
152
153 /* ************************************************************************* */
154 /* Write the VTK Structured Coordinates ************************************ */
155 /* ************************************************************************* */
156 void
write_vtk_structured_coordinates(int fd,FILE * fp,char * varname,struct Cell_head region,int out_type,char * null_value,double scale,int dp)157 write_vtk_structured_coordinates(int fd, FILE * fp, char *varname,
158 struct Cell_head region, int out_type,
159 char *null_value, double scale, int dp)
160 {
161 int ncols = region.cols;
162 int nrows = region.rows;
163 int row, col;
164 int rowcount = 0, colcount = 0;
165 double nspos = 0.0, ewpos = 0.0;
166 double nullvalue, value;
167 void *ptr, *raster;
168
169 G_debug(3, _("write_vtk_structured_coordinates: Writing Coordinates"));
170
171 /*the nullvalue */
172 if (!sscanf(null_value, "%lf", &nullvalue)) {
173 G_warning("Null value is not valid, using 0 instead.");
174 nullvalue = 0;
175 }
176
177 raster = Rast_allocate_buf(out_type);
178
179 rowcount = 0;
180 for (row = nrows - 1; row >= 0; row--) {
181 colcount = 0;
182 G_percent((row - nrows) * (-1), nrows, 2);
183
184 Rast_get_row(fd, raster, row, out_type);
185
186 nspos = region.ns_res / 2 + region.south + rowcount * region.ns_res;
187 nspos -= y_extent;
188
189 for (col = 0, ptr = raster; col < ncols;
190 col++, ptr = G_incr_void_ptr(ptr, Rast_cell_size(out_type))) {
191 ewpos =
192 region.ew_res / 2 + region.west + colcount * region.ew_res;
193 ewpos -= x_extent;
194
195 value = get_raster_value_as_double(out_type, ptr, nullvalue);
196 fprintf(fp, "%.*f %.*f %.*f\n", dp, ewpos, dp, nspos, dp,
197 value * scale);
198
199 colcount++;
200 }
201 rowcount++;
202 }
203 return;
204 }
205
206 /* ************************************************************************* */
207 /* Write Polygonal Coordinates ********************************************* */
208 /* ************************************************************************* */
209 void
write_vtk_polygonal_coordinates(int fd,FILE * fp,char * varname,struct Cell_head region,int out_type,char * null_value,double scale,int polytype,int dp)210 write_vtk_polygonal_coordinates(int fd, FILE * fp, char *varname,
211 struct Cell_head region, int out_type,
212 char *null_value, double scale, int polytype,
213 int dp)
214 {
215 int ncols = region.cols;
216 int nrows = region.rows;
217 int row, col;
218 int rowcount = 0, colcount = 0;
219 double nspos = 0.0, ewpos = 0.0;
220 double value, nullvalue;
221 void *ptr, *raster;
222 int i, j, count;
223
224 G_debug(3,
225 _("write_vtk_polygonal_coordinates: Writing VTK Polygonal data"));
226
227 /*the nullvalue */
228 if (!sscanf(null_value, "%lf", &nullvalue)) {
229 G_warning("Null value is not valid, using 0 instead.");
230 nullvalue = 0;
231 }
232
233
234 /*First we are writing the coordinate points, the elevation cell is only used for the z coordinate */
235 raster = Rast_allocate_buf(out_type);
236
237 rowcount = 0;
238 for (row = nrows - 1; row >= 0; row--) {
239 colcount = 0;
240 G_percent((row - nrows) * (-1), nrows, 10);
241
242 Rast_get_row(fd, raster, row, out_type);
243
244 nspos = region.ns_res / 2 + region.south + rowcount * region.ns_res;
245 nspos -= y_extent;
246
247 for (col = 0, ptr = raster; col < ncols;
248 col++, ptr = G_incr_void_ptr(ptr, Rast_cell_size(out_type))) {
249 ewpos =
250 region.ew_res / 2 + region.west + colcount * region.ew_res;
251 ewpos -= x_extent;
252
253 value = get_raster_value_as_double(out_type, ptr, nullvalue);
254 fprintf(fp, "%.*f %.*f %.*f\n", dp, ewpos, dp, nspos, dp,
255 value * scale);
256
257 colcount++;
258 }
259 rowcount++;
260 }
261
262 /*Now we need to write the Connectivity between the points */
263
264 if (polytype == QUADS) { /*The default */
265 /*If Datafiltering should be supportet, we use Polygons to represent the grid */
266 fprintf(fp, "POLYGONS %i %i\n", (region.rows - 1) * (region.cols - 1),
267 5 * (region.rows - 1) * (region.cols - 1));
268
269 /*We creat a grid of quads, the corners of the quads are the datapoints */
270 for (i = 0; i < region.rows - 1; i++) {
271 for (j = 0; j < region.cols - 1; j++) {
272 fprintf(fp, "4 %i %i %i %i \n", i * region.cols + j,
273 i * region.cols + j + 1,
274 (i + 1) * region.cols + j + 1,
275 (i + 1) * region.cols + j);
276 }
277 }
278 }
279
280
281 if (polytype == TRIANGLE_STRIPS) {
282 /*TriangleStrips, take a look ate www.vtk.org for the definition of triangle_strips in vtk */
283 /*If we use triangelestrips, the number of points per strip is equal to the double number of cols we have */
284 fprintf(fp, "TRIANGLE_STRIPS %i %i\n", region.rows - 1,
285 (region.rows - 1) + (region.rows - 1) * (2 * region.cols));
286
287 /*For every Row-1 make a strip */
288 for (i = 0; i < region.rows - 1; i++) {
289
290 /*Number of Points */
291 fprintf(fp, "%i ", region.cols * 2);
292
293 /*Write the points */
294 for (j = 0; j < region.cols; j++) {
295 fprintf(fp, "%i %i ", i * region.cols + j,
296 (i + 1) * region.cols + j);
297 }
298 fprintf(fp, "\n");
299 }
300 }
301
302 if (polytype == VERTICES) {
303 /*If the Data should be Triangulated by VTK, we use vertices to represent the grid */
304 fprintf(fp, "VERTICES %i %i\n", region.rows,
305 region.rows + (region.rows) * (region.cols));
306
307 count = 0;
308 for (i = 0; i < region.rows; i++) {
309 fprintf(fp, "%i ", (region.cols));
310 for (j = 0; j < region.cols; j++) {
311 fprintf(fp, "%i ", count);
312 count++;
313 }
314 fprintf(fp, "\n");
315 }
316 }
317
318 return;
319 }
320
321 /* ************************************************************************* */
322 /* Write the VTK Data ****************************************************** */
323 /* ************************************************************************* */
324 void
write_vtk_data(int fd,FILE * fp,char * varname,struct Cell_head region,int out_type,char * null_value,int dp)325 write_vtk_data(int fd, FILE * fp, char *varname, struct Cell_head region,
326 int out_type, char *null_value, int dp)
327 {
328 int ncols = region.cols;
329 int nrows = region.rows;
330 int row, col;
331 double value, nullvalue;
332 void *ptr, *raster;
333
334 G_debug(3, _("write_vtk_data: Writing VTK-Data"));
335
336 /*the nullvalue */
337 if (!sscanf(null_value, "%lf", &nullvalue)) {
338 G_warning("Null value is not valid, using 0 instead.");
339 nullvalue = 0;
340 }
341
342 fprintf(fp, "SCALARS %s float 1\n", varname);
343 fprintf(fp, "LOOKUP_TABLE default\n");
344
345
346 raster = Rast_allocate_buf(out_type);
347
348 for (row = nrows - 1; row >= 0; row--) {
349 G_percent((row - nrows) * (-1), nrows, 10);
350
351 Rast_get_row(fd, raster, row, out_type);
352
353 for (col = 0, ptr = raster; col < ncols;
354 col++, ptr = G_incr_void_ptr(ptr, Rast_cell_size(out_type))) {
355
356 value = get_raster_value_as_double(out_type, ptr, nullvalue);
357 fprintf(fp, "%.*f ", dp, value);
358
359 }
360 fprintf(fp, "\n");
361 }
362 return;
363
364 }
365
366
367
368 /* ************************************************************************* */
369 /* Write the VTK RGB Image Data ******************************************** */
370 /* ************************************************************************* */
371 void
write_vtk_rgb_image_data(int redfd,int greenfd,int bluefd,FILE * fp,const char * varname,struct Cell_head region,int out_type,int dp)372 write_vtk_rgb_image_data(int redfd, int greenfd, int bluefd, FILE * fp,
373 const char *varname, struct Cell_head region,
374 int out_type, int dp)
375 {
376 int ncols = region.cols;
377 int nrows = region.rows;
378 int row, col;
379 void *redptr, *redraster;
380 void *greenptr, *greenraster;
381 void *blueptr, *blueraster;
382 double r = 0.0, g = 0.0, b = 0.0;
383
384 G_debug(3, _("write_vtk_rgb_image_data: Writing VTK-ImageData"));
385
386 fprintf(fp, "COLOR_SCALARS %s 3\n", varname);
387
388 redraster = Rast_allocate_buf(out_type);
389 greenraster = Rast_allocate_buf(out_type);
390 blueraster = Rast_allocate_buf(out_type);
391
392 for (row = nrows - 1; row >= 0; row--) {
393 G_percent((row - nrows) * (-1), nrows, 10);
394
395 Rast_get_row(redfd, redraster, row, out_type);
396 Rast_get_row(greenfd, greenraster, row, out_type);
397 Rast_get_row(bluefd, blueraster, row, out_type);
398
399 for (col = 0, redptr = redraster, greenptr = greenraster, blueptr =
400 blueraster; col < ncols;
401 col++, redptr =
402 G_incr_void_ptr(redptr, Rast_cell_size(out_type)), greenptr =
403 G_incr_void_ptr(greenptr, Rast_cell_size(out_type)), blueptr =
404 G_incr_void_ptr(blueptr, Rast_cell_size(out_type))) {
405
406 r = get_raster_value_as_double(out_type, redptr, 0.0);
407 g = get_raster_value_as_double(out_type, greenptr, 0.0);
408 b = get_raster_value_as_double(out_type, blueptr, 0.0);
409
410 /*Test of valuerange, the data should be 1 byte gray values */
411 if (r > 255 || g > 255 || b > 255 || r < 0 || g < 0 || b < 0) {
412 G_warning(_("Wrong map values! Values should in between 0 and 255!\n"));
413 fprintf(fp, "0 0 0 \n");
414 }
415 else {
416 fprintf(fp, "%.*f %.*f %.*f \n", dp, r / 255, dp, g / 255, dp,
417 b / 255);
418 }
419 }
420 }
421 return;
422
423 }
424
425
426 /* ************************************************************************* */
427 /* Write the VTK Vector Data *********************************************** */
428 /* ************************************************************************* */
429 void
write_vtk_vector_data(int xfd,int yfd,int zfd,FILE * fp,const char * varname,struct Cell_head region,int out_type,int dp)430 write_vtk_vector_data(int xfd, int yfd, int zfd, FILE * fp,
431 const char *varname, struct Cell_head region,
432 int out_type, int dp)
433 {
434 int ncols = region.cols;
435 int nrows = region.rows;
436 int row, col;
437 void *xptr, *xraster;
438 void *yptr, *yraster;
439 void *zptr, *zraster;
440 double x = 0.0, y = 0.0, z = 0.0;
441
442 G_debug(3, _("write_vtk_vector_data: Writing VTK-vector data"));
443
444 fprintf(fp, "VECTORS %s float\n", varname);
445
446 xraster = Rast_allocate_buf(out_type);
447 yraster = Rast_allocate_buf(out_type);
448 zraster = Rast_allocate_buf(out_type);
449
450 for (row = nrows - 1; row >= 0; row--) {
451 G_percent((row - nrows) * (-1), nrows, 10);
452
453 Rast_get_row(xfd, xraster, row, out_type);
454 Rast_get_row(yfd, yraster, row, out_type);
455 Rast_get_row(zfd, zraster, row, out_type);
456
457 for (col = 0, xptr = xraster, yptr = yraster, zptr =
458 zraster; col < ncols;
459 col++, xptr =
460 G_incr_void_ptr(xptr, Rast_cell_size(out_type)), yptr =
461 G_incr_void_ptr(yptr, Rast_cell_size(out_type)), zptr =
462 G_incr_void_ptr(zptr, Rast_cell_size(out_type))) {
463
464 x = get_raster_value_as_double(out_type, xptr, 0.0);
465 y = get_raster_value_as_double(out_type, yptr, 0.0);
466 z = get_raster_value_as_double(out_type, zptr, 0.0);
467
468 fprintf(fp, "%.*f %.*f %.*f \n", dp, x, dp, y, dp, z);
469 }
470 }
471 return;
472
473 }
474