/* r.out.pov v1.0 beta convert grass raster to tga height-field for POVray. Persistence of Vision (POV) raytracer can use a height-field defined in a Targa (.TGA) image file format where the RGB pixel values are 24 bits (3 bytes). A 16 bit unsigned integer height-field value is assigned as follows: RED high byte GREEN low byte BLUE empty This code is based on code by W.D. Kirby who wrote DEM2XYZ in 1996. Author: Klaus Meyer, GEUM.tec GbR, eMail: GEUM.tec@geum.de Date: July 1998 */ #include #include #include #include #include #include void writeHeader(FILE *outf); void processProfiles(int inputFile, FILE *outputF); #define YMAX 65536 /* max length scan line */ #define XMAX 65536 /* max # of scan lines */ #define MIN(x,y) (((x) < (y)) ? (x) : (y)) #define MAX(x,y) (((x) > (y)) ? (x) : (y)) #define SW 0 #define NW 1 #define NE 2 #define SE 3 int base[XMAX]; /* array of base elevations */ double defaultElev = 0; /* elevation for empty points */ int width, height; /* height-field image width and height */ int hfType = 0; /* height-field type */ double hfBias, hfNorm = 1.0; double verticalScale=1.0; /* to stretch or shrink elevations */ double minValue=50000., maxValue=-50000.; char mapLabel[145]; int DEMlevel, elevationPattern, groundSystem, groundZone; double projectParams[15]; int planeUnitOfMeasure, elevUnitOfMeasure, polygonSizes; double groundCoords[4][2], elevBounds[2], localRotation; int accuracyCode; double spatialResolution[3]; int profileDimension[2]; int firstRow, lastRow; int wcount = 0, hcount; double deltaY; char inText[24], **junk; double eastMost, westMost, southMost, northMost; int eastMostSample, westMostSample, southMostSample, northMostSample; int rowCount, columnCount, r, c; long int cellCount = 0 ; int rowStr,rowEnd, colStr, colEnd; int main(int argc, char *argv[]) { struct Cell_head region; struct Range range; CELL range_min, range_max; FILE *outf; char *outfilename; CELL *cell; char *name; char *mapset; int fd; int nrows, ncols; double bias; char fmt[20]; struct GModule *module; struct { struct Option *map ; struct Option *tga ; struct Option *hftype ; struct Option *bias ; struct Option *scaleFactor ; } parm; G_gisinit (argv[0]); module = G_define_module(); module->keywords = _("raster"); module->description = _("Converts a raster map layer into a height-field file for POVRAY."); /* Define the different options */ parm.map = G_define_option() ; parm.map->key = "map"; parm.map->type = TYPE_STRING; parm.map->required = YES; parm.map->gisprompt = "old,cell,raster" ; parm.map->description= _("Name of an existing raster map") ; parm.tga = G_define_option() ; parm.tga->key = "tga"; parm.tga->type = TYPE_STRING; parm.tga->required = YES; parm.tga->description= _("Name of output povray file (TGA height field file)") ; parm.hftype = G_define_option() ; parm.hftype->key = "hftype"; parm.hftype->type = TYPE_INTEGER; parm.hftype->required = NO; parm.hftype->description= _("Height-field type (0=actual heights 1=normalized)"); parm.bias = G_define_option() ; parm.bias->key = "bias"; parm.bias->type = TYPE_DOUBLE; parm.bias->required = NO; parm.bias->description= _("Elevation bias"); parm.scaleFactor = G_define_option() ; parm.scaleFactor->key = "scale"; parm.scaleFactor->type = TYPE_DOUBLE; parm.scaleFactor->required = NO; parm.scaleFactor->description= _("Vertical scaling factor"); if (G_parser(argc, argv)) exit(EXIT_FAILURE); strcpy (fmt, "%ld "); if (parm.hftype->answer != NULL) { sscanf (parm.hftype->answer, "%d", &hfType); if (hfType > 0) sprintf (fmt, "%%%dld ", hfType); else sprintf (fmt, "%%%dld", hfType); } strcpy(fmt, "%lf "); if (parm.bias->answer != NULL) { sscanf (parm.bias->answer, "%lf", &bias); if (bias > 0) sprintf (fmt, "%%%dlf ", bias); else sprintf (fmt, "%%%dlf", bias); } strcpy(fmt, "%lf "); if (parm.scaleFactor->answer != NULL) { sscanf (parm.scaleFactor->answer, "%lf", &verticalScale); if (verticalScale > 0) sprintf (fmt, "%%%dlf ", verticalScale); else sprintf (fmt, "%%%dlf", verticalScale); } name = parm.map->answer; mapset = G_find_cell2 (name, ""); if (mapset == NULL) G_fatal_error(_("Raster map <%s> not found"), name); fd = G_open_cell_old (name, mapset); if (fd < 0) G_fatal_error(_("Unable to open raster map <%s>"), name); outfilename = parm.tga->answer; if (outfilename == NULL) G_fatal_error(_("Invalid output filename <%s>"), outfilename); if (NULL == (outf = fopen (outfilename, "wb"))) G_fatal_error(_("Unable to open output file <%s>"), outfilename); cell = G_allocate_cell_buf(); nrows = G_window_rows(); ncols = G_window_cols(); if (nrows > YMAX || ncols > XMAX) G_fatal_error(_("Raster map is too big! Exceeds %d columns or %d rows"), XMAX, YMAX); columnCount = ncols; rowCount = nrows; width = ncols; height = nrows; /* SW, NW, NE, SE corners */ G_get_window(®ion); eastMost = region.east; westMost = region.west; northMost = region.north; southMost = region.south; G_init_range(&range); G_read_range(name, mapset, &range); G_get_range_min_max (&range, &range_min, &range_max); if (range.min < 0 || range.max < 0) G_warning(_("Negative elevation values in input")); elevBounds[0] = range.min; elevBounds[1] = range.max; /* Normalize using max value of unsigned short integer (16 bit) */ if (hfType == 1) hfNorm = 65535.0 / (elevBounds[1]+hfBias); spatialResolution[0] = region.ew_res; spatialResolution[1] = region.ew_res; spatialResolution[2] = region.ns_res; /* write TGA image file header */ (void) writeHeader(outf); /* have to read everything */ (void) processProfiles(fd,outf); fclose (outf); G_close_cell(fd); exit(EXIT_SUCCESS); } void writeHeader (FILE *outputF) { int i; /* Write TGA image header */ for (i = 0; i < 10; i++) /* 00, 00, 02, then 7 00's... */ if (i == 2) putc((short)i, outputF); else putc(0, outputF); putc(0, outputF); /* y origin set to "First_Line" */ putc(0, outputF); putc((short)(width % 256), outputF); /* write width and height */ putc((short)(width / 256), outputF); putc((short)(height % 256), outputF); putc((short)(height / 256), outputF); putc(24, outputF); /* 24 bits/pixel (16 million colors!) */ putc(32, outputF); /* Bitmask, pertinent bit: top-down raster */ } /* * read profiles */ void processProfiles (int inputFile, FILE *outputF) { CELL *cell; int c, r; double tempFloat; cell = G_allocate_cell_buf(); for (r = 0; r < rowCount; r++) { if (G_get_map_row (inputFile, cell, r) < 0) exit(1); /* break; */ for (c = 0; c < columnCount; c++) { tempFloat = ((float) cell[c] * verticalScale) + hfBias; /* Normalize */ tempFloat *= hfNorm; if(tempFloat < 0.0) tempFloat = 0.0; if(tempFloat > 65535.0) tempFloat = 65535.0; if(tempFloat>maxValue) maxValue = tempFloat; if(tempFloat