[GRASS-SVN] r51526 - in grass-addons/grass7/raster: . r.vol.dem
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 24 17:23:22 EDT 2012
Author: neteler
Date: 2012-04-24 14:23:22 -0700 (Tue, 24 Apr 2012)
New Revision: 51526
Added:
grass-addons/grass7/raster/r.vol.dem/
grass-addons/grass7/raster/r.vol.dem/Makefile
grass-addons/grass7/raster/r.vol.dem/chull.c
grass-addons/grass7/raster/r.vol.dem/chull.h
grass-addons/grass7/raster/r.vol.dem/globals.h
grass-addons/grass7/raster/r.vol.dem/macros.h
grass-addons/grass7/raster/r.vol.dem/main.c
Log:
Benjamin Ducke: new (few G7 updates lacking)
Added: grass-addons/grass7/raster/r.vol.dem/Makefile
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/Makefile 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,4 @@
+MODULE_TOPDIR =.. /..PGM = r.vol.dem LIBES = $(GISLIB) $(DATETIMELIB)
+EXTRA_CFLAGS = -I.. /.. / include include $(MODULE_TOPDIR) / include / Make / Module.make default:
+
+cmd
Added: grass-addons/grass7/raster/r.vol.dem/chull.c
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/chull.c (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/chull.c 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,1283 @@
+/*
+ This code is described in "Computational Geometry in C" (Second Edition),
+ Chapter 4. It is not written to be comprehensible without the
+ explanation in that book.
+
+ Input: 3n integer coordinates for the points.
+ Output: the 3D convex hull, in postscript with embedded comments
+ showing the vertices and faces.
+
+ Compile: gcc -o chull chull.c
+
+ Written by Joseph O'Rourke, with contributions by
+ Kristy Anderson, John Kutcher, Catherine Schevon, Susan Weller.
+ Last modified: March 1998
+ Questions to orourke at cs.smith.edu.
+ --------------------------------------------------------------------
+ This code is Copyright 1998 by Joseph O'Rourke. It may be freely
+ redistributed in its entirety provided that this copyright notice is
+ not removed.
+ --------------------------------------------------------------------
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <grass/gis.h>
+
+#include "globals.h"
+
+/*Define Boolean type */
+typedef enum
+{ BFALSE, BTRUE } bool;
+
+/* Define vertex indices. */
+#define X 0
+#define Y 1
+#define Z 2
+
+/* Define structures for vertices, edges and faces */
+typedef struct tVertexStructure tsVertex;
+typedef tsVertex *tVertex;
+
+typedef struct tEdgeStructure tsEdge;
+typedef tsEdge *tEdge;
+
+typedef struct tFaceStructure tsFace;
+typedef tsFace *tFace;
+
+struct tVertexStructure
+{
+ long int v[3];
+ int vnum;
+ tEdge duplicate; /* pointer to incident cone edge (or NULL) */
+ bool onhull; /* T iff point on hull. */
+ bool mark; /* T iff point already processed. */
+ tVertex next, prev;
+};
+
+struct tEdgeStructure
+{
+ tFace adjface[2];
+ tVertex endpts[2];
+ tFace newface; /* pointer to incident cone face. */
+ bool delete; /* T iff edge should be delete. */
+ tEdge next, prev;
+};
+
+struct tFaceStructure
+{
+ tEdge edge[3];
+ tVertex vertex[3];
+ bool visible; /* T iff face visible from new point. */
+ tFace next, prev;
+};
+
+/* Define flags */
+#define ONHULL BTRUE
+#define REMOVED BTRUE
+#define VISIBLE BTRUE
+#define PROCESSED BTRUE
+#define SAFE 2000000000 /* Range of safe coord values. */
+
+/* Global variable definitions */
+tVertex vertices = NULL;
+tEdge edges = NULL;
+tFace faces = NULL;
+bool debug = BFALSE;
+bool check = BFALSE;
+
+/* Function declarations */
+tVertex MakeNullVertex(void);
+void ReadVertices(long int *px, long int *py, long int *pz, int num_points);
+void Print(void);
+void Dump(FILE * tmpfile);
+void SubVec(long int a[3], long int b[3], long int c[3]);
+int DoubleTriangle(void);
+void ConstructHull(void);
+bool AddOne(tVertex p);
+int VolumeSign(tFace f, tVertex p);
+int Volumei(tFace f, tVertex p);
+tFace MakeConeFace(tEdge e, tVertex p);
+void MakeCcw(tFace f, tEdge e, tVertex p);
+tEdge MakeNullEdge(void);
+tFace MakeNullFace(void);
+tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace f);
+void CleanUp(void);
+void CleanEdges(void);
+void CleanFaces(void);
+void CleanVertices(void);
+bool Collinear(tVertex a, tVertex b, tVertex c);
+void CheckEuler(long int V, long int E, long int F);
+void PrintPoint(tVertex p);
+void Checks(void);
+void Consistency(void);
+void Convexity(void);
+void PrintOut(tVertex v);
+void PrintVertices(void);
+void PrintEdges(void);
+void PrintFaces(void);
+
+#include "macros.h"
+
+
+/*
+
+ Release all memory allocated for edges, faces and vertices
+
+ */
+void freeMem(void)
+{
+ tEdge e; /* Primary index into edge list. */
+ tFace f; /* Primary pointer into face list. */
+ tVertex v;
+ tEdge te; /* Temporary edge pointer. */
+ tFace tf; /* Temporary face pointer. */
+ tVertex tv; /* Temporary vertex pointer. */
+
+ if (DEBUG) {
+ fprintf(stdout, "FREE MEM:\n");
+ fflush(stdout);
+ }
+
+ if (DEBUG) {
+ fprintf(stdout, " EDGES:\n");
+ fflush(stdout);
+ }
+ e = edges;
+ do {
+ te = e;
+ e = e->next;
+ DELETE(edges, te);
+ } while (e != edges);
+
+ if (DEBUG) {
+ fprintf(stdout, " FACES:\n");
+ fflush(stdout);
+ }
+ f = faces;
+ do {
+ tf = f;
+ f = f->next;
+ DELETE(faces, tf);
+ } while (f != faces);
+
+ if (DEBUG) {
+ fprintf(stdout, " VERTICES:\n");
+ fflush(stdout);
+ }
+ v = vertices;
+ do {
+ tv = v;
+ v = v->next;
+ DELETE(vertices, tv);
+ } while (v != vertices);
+
+ FREE(te);
+ FREE(tf);
+ FREE(tv);
+
+ DELETE(edges, e);
+ DELETE(faces, f);
+ DELETE(vertices, v);
+
+ FREE(edges);
+ FREE(faces);
+ FREE(vertices);
+
+ if (DEBUG) {
+ fprintf(stdout, "MEM FREE'D!\n");
+ fflush(stdout);
+ }
+
+}
+
+
+/*-------------------------------------------------------------------*/
+int make3DHull(long int *px, long int *py, long int *pz, int num_points,
+ FILE * tmpfile)
+{
+ int error;
+
+ check = BFALSE;
+ debug = BFALSE;
+
+ if (DEBUG > 1)
+ check = BTRUE;
+ if (DEBUG > 2)
+ debug = BTRUE;
+
+ ReadVertices(px, py, pz, num_points);
+
+ error = DoubleTriangle();
+ if (error < 0) {
+ G_warning("All points of this layer are in the same voxel plane.\n");
+ freeMem();
+ return (error);
+ }
+
+ ConstructHull();
+
+ Dump(tmpfile);
+
+ freeMem();
+
+ return (0);
+}
+
+/*---------------------------------------------------------------------
+MakeNullVertex: Makes a vertex, nulls out fields.
+---------------------------------------------------------------------*/
+tVertex MakeNullVertex(void)
+{
+ tVertex v;
+
+ NEW(v, tsVertex);
+ v->duplicate = NULL;
+ v->onhull = !ONHULL;
+ v->mark = !PROCESSED;
+ ADD(vertices, v);
+
+ return v;
+}
+
+/*---------------------------------------------------------------------
+ReadVertices: Reads in the vertices, and links them into a circular
+list with MakeNullVertex. There is no need for the # of vertices to be
+the first line: the function looks for EOF instead. Sets the global
+variable vertices via the ADD macro.
+---------------------------------------------------------------------*/
+void ReadVertices(long int *px, long int *py, long int *pz, int num_points)
+{
+ tVertex v;
+ int vnum = 0;
+ int i;
+
+ for (i = 0; i < num_points; i++) {
+ v = MakeNullVertex();
+ v->v[X] = px[i];
+ v->v[Y] = py[i];
+ v->v[Z] = pz[i];
+ v->vnum = vnum++;
+ if ((abs(px[i]) > SAFE) || (abs(py[i]) > SAFE) || (abs(pz[i]) > SAFE)) {
+ printf
+ ("Coordinate of vertex below might be too large: run with -c flag\n");
+ PrintPoint(v);
+ }
+ }
+}
+
+/*---------------------------------------------------------------------
+Dump: Dumps out the vertices and the faces to a file.
+Uses the vnum indices corresponding to the order in which the vertices
+were input.
+Output is in GRASS ASCII file format.
+---------------------------------------------------------------------*/
+void Dump(FILE * tmpfile)
+{
+ /* Pointers to vertices, edges, faces. */
+ tVertex v;
+ tEdge e;
+ tFace f;
+ double dx, dy, dz;
+ long int a[3], b[3]; /* used to compute normal vector */
+
+ /* Counters for Euler's formula. */
+ long int V = 0, E = 0, F = 0;
+
+ /* Note: lowercase==pointer, uppercase==counter. */
+ long int cat;
+
+ f = faces;
+ do {
+ ++F;
+ f = f->next;
+ } while (f != faces);
+
+ /* GRASS 6 map header */
+ fprintf(tmpfile, "ORGANIZATION: \n");
+ fprintf(tmpfile, "DIGIT DATE: \n");
+ fprintf(tmpfile, "DIGIT NAME: \n");
+ fprintf(tmpfile, "MAP NAME: \n");
+ fprintf(tmpfile, "MAP DATE: \n");
+ fprintf(tmpfile, "MAP SCALE: 10000\n");
+ fprintf(tmpfile, "OTHER INFO: %li faces.\n", F);
+ fprintf(tmpfile, "ZONE: 0\n");
+ fprintf(tmpfile, "MAP THRESH: 0.5\n");
+ fprintf(tmpfile, "VERTI:\n");
+
+ cat = 1;
+
+ /* putting all faces into one object does not produce nice output in NVIZ !
+ If we decided to dump all layers into a single file, we would also need
+ a cat for each object to be passed from the main program.
+ */
+
+ /*
+ printf("F %li 1\n", F*4);
+ do {
+ dx = ((double) ( f->vertex[0]->v[X] ) / (1000));
+ dy = ((double) ( f->vertex[0]->v[Y] ) / (1000));
+ dz = ((double) ( f->vertex[0]->v[Z] ) / (1000));
+ printf(" %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double) ( f->vertex[1]->v[X] ) / (1000));
+ dy = ((double) ( f->vertex[1]->v[Y] ) / (1000));
+ dz = ((double) ( f->vertex[1]->v[Z] ) / (1000));
+ printf(" %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double) ( f->vertex[2]->v[X] ) / (1000));
+ dy = ((double) ( f->vertex[2]->v[Y] ) / (1000));
+ dz = ((double) ( f->vertex[2]->v[Z] ) / (1000));
+ printf(" %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double) ( f->vertex[0]->v[X] ) / (1000));
+ dy = ((double) ( f->vertex[0]->v[Y] ) / (1000));
+ dz = ((double) ( f->vertex[0]->v[Z] ) / (1000));
+ printf(" %.3f %.3f %.3f\n", dx, dy, dz);
+ f = f->next;
+ } while ( f != faces );
+ printf(" %li 1\n", cat);
+ */
+
+ do {
+ fprintf(tmpfile, "F 4 1\n");
+ dx = ((double)(f->vertex[0]->v[X]) / (PRECISION));
+ dy = ((double)(f->vertex[0]->v[Y]) / (PRECISION));
+ dz = ((double)(f->vertex[0]->v[Z]) / (PRECISION));
+ fprintf(tmpfile, " %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double)(f->vertex[1]->v[X]) / (PRECISION));
+ dy = ((double)(f->vertex[1]->v[Y]) / (PRECISION));
+ dz = ((double)(f->vertex[1]->v[Z]) / (PRECISION));
+ fprintf(tmpfile, " %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double)(f->vertex[2]->v[X]) / (PRECISION));
+ dy = ((double)(f->vertex[2]->v[Y]) / (PRECISION));
+ dz = ((double)(f->vertex[2]->v[Z]) / (PRECISION));
+ fprintf(tmpfile, " %.3f %.3f %.3f\n", dx, dy, dz);
+ dx = ((double)(f->vertex[0]->v[X]) / (PRECISION));
+ dy = ((double)(f->vertex[0]->v[Y]) / (PRECISION));
+ dz = ((double)(f->vertex[0]->v[Z]) / (PRECISION));
+ fprintf(tmpfile, " %.3f %.3f %.3f\n", dx, dy, dz);
+ fprintf(tmpfile, " %li 1\n", cat);
+ cat++;
+ f = f->next;
+ } while (f != faces);
+
+ if (DEBUG > 0) {
+ fprintf(stdout, "3D Convex hull check (Euler):\n");
+ check = BTRUE;
+ CheckEuler(V, E, F);
+ }
+
+
+
+}
+
+
+
+/*---------------------------------------------------------------------
+Print: Prints out the vertices and the faces. Uses the vnum indices
+corresponding to the order in which the vertices were input.
+Output is in PostScript format.
+---------------------------------------------------------------------*/
+void Print(void)
+{
+ /* Pointers to vertices, edges, faces. */
+ tVertex v;
+ tEdge e;
+ tFace f;
+ long int xmin, ymin, xmax, ymax;
+ long int a[3], b[3]; /* used to compute normal vector */
+
+ /* Counters for Euler's formula. */
+ long int V = 0, E = 0, F = 0;
+
+ /* Note: lowercase==pointer, uppercase==counter. */
+
+ /*-- find X min & max --*/
+ v = vertices;
+ xmin = xmax = v->v[X];
+ do {
+ if (v->v[X] > xmax)
+ xmax = v->v[X];
+ else if (v->v[X] < xmin)
+ xmin = v->v[X];
+ v = v->next;
+ } while (v != vertices);
+
+ /*-- find Y min & max --*/
+ v = vertices;
+ ymin = ymax = v->v[Y];
+ do {
+ if (v->v[Y] > ymax)
+ ymax = v->v[Y];
+ else if (v->v[Y] < ymin)
+ ymin = v->v[Y];
+ v = v->next;
+ } while (v != vertices);
+
+ /* PostScript header */
+ printf("%%!PS\n");
+ printf("%%%%BoundingBox: %li %li %li %li\n", xmin, ymin, xmax, ymax);
+ printf(".00 .00 setlinewidth\n");
+ printf("%li %li translate\n", -xmin + 72, -ymin + 72);
+ /* The +72 shifts the figure one inch from the lower left corner */
+
+ /* Vertices. */
+ v = vertices;
+ do {
+ if (v->mark)
+ V++;
+ v = v->next;
+ } while (v != vertices);
+ printf("\n%%%% Vertices:\tV = %li\n", V);
+ printf("%%%% index:\tx\ty\tz\n");
+ do {
+ printf("%%%% %5d:\t%li\t%li\t%li\n",
+ v->vnum, v->v[X], v->v[Y], v->v[Z]);
+ v = v->next;
+ } while (v != vertices);
+
+ /* Faces. */
+ /* visible faces are printed as PS output */
+ f = faces;
+ do {
+ ++F;
+ f = f->next;
+ } while (f != faces);
+ printf("\n%%%% Faces:\tF = %li\n", F);
+ printf("%%%% Visible faces only: \n");
+ do {
+ /* Print face only if it is visible: if normal vector >= 0 */
+ SubVec(f->vertex[1]->v, f->vertex[0]->v, a);
+ SubVec(f->vertex[2]->v, f->vertex[1]->v, b);
+ if ((a[0] * b[1] - a[1] * b[0]) >= 0) {
+ printf("%%%% vnums: %d %d %d\n",
+ f->vertex[0]->vnum,
+ f->vertex[1]->vnum, f->vertex[2]->vnum);
+ printf("newpath\n");
+ printf("%li\t%li\tmoveto\n",
+ f->vertex[0]->v[X], f->vertex[0]->v[Y]);
+ printf("%li\t%li\tlineto\n",
+ f->vertex[1]->v[X], f->vertex[1]->v[Y]);
+ printf("%li\t%li\tlineto\n",
+ f->vertex[2]->v[X], f->vertex[2]->v[Y]);
+ printf("closepath stroke\n\n");
+ }
+ f = f->next;
+ } while (f != faces);
+
+ /* prints a list of all faces */
+ printf("%%%% List of all faces: \n");
+ printf("%%%%\tv0\tv1\tv2\t(vertex indices)\n");
+ do {
+ printf("%%%%\t%d\t%d\t%d\n",
+ f->vertex[0]->vnum, f->vertex[1]->vnum, f->vertex[2]->vnum);
+ f = f->next;
+ } while (f != faces);
+
+ /* Edges. */
+ e = edges;
+ do {
+ E++;
+ e = e->next;
+ } while (e != edges);
+ printf("\n%%%% Edges:\tE = %li\n", E);
+ /* Edges not printed out (but easily added). */
+
+ printf("\nshowpage\n\n");
+
+ if (DEBUG > 0) {
+ fprintf(stdout, "3D Convex hull check (Euler):\n");
+ check = BTRUE;
+ CheckEuler(V, E, F);
+ }
+
+}
+
+/*---------------------------------------------------------------------
+SubVec: Computes a - b and puts it into c.
+---------------------------------------------------------------------*/
+void SubVec(long int a[3], long int b[3], long int c[3])
+{
+ long int i;
+
+ for (i = 0; i < 2; i++)
+ c[i] = a[i] - b[i];
+
+}
+
+/*---------------------------------------------------------------------
+ DoubleTriangle builds the initial double triangle. It first finds 3
+ noncollinear points and makes two faces out of them, in opposite order.
+ It then finds a fourth point that is not coplanar with that face. The
+ vertices are stored in the face structure in counterclockwise order so
+ that the volume between the face and the point is negative. Lastly, the
+ 3 newfaces to the fourth point are constructed and the data structures
+ are cleaned up.
+---------------------------------------------------------------------*/
+
+/* RETURN: 0 if OK */
+/* -1 if all points collinear */
+/* -2 if all points coplanar */
+
+int DoubleTriangle(void)
+{
+ tVertex v0, v1, v2, v3, t;
+ tFace f0, f1 = NULL;
+ tEdge e0, e1, e2, s;
+ long int vol;
+
+ /* Find 3 noncollinear points. */
+ v0 = vertices;
+ while (Collinear(v0, v0->next, v0->next->next)) {
+ if ((v0 = v0->next) == vertices) {
+ if (debug) {
+ printf("DoubleTriangle: All points are collinear!\n");
+ }
+ return (-1);
+ }
+ }
+ v1 = v0->next;
+ v2 = v1->next;
+
+ /* Mark the vertices as processed. */
+ v0->mark = PROCESSED;
+ v1->mark = PROCESSED;
+ v2->mark = PROCESSED;
+
+ /* Create the two "twin" faces. */
+ f0 = MakeFace(v0, v1, v2, f1);
+ f1 = MakeFace(v2, v1, v0, f0);
+
+ /* Link adjacent face fields. */
+ f0->edge[0]->adjface[1] = f1;
+ f0->edge[1]->adjface[1] = f1;
+ f0->edge[2]->adjface[1] = f1;
+ f1->edge[0]->adjface[1] = f0;
+ f1->edge[1]->adjface[1] = f0;
+ f1->edge[2]->adjface[1] = f0;
+
+ /* Find a fourth, noncoplanar point to form tetrahedron. */
+ v3 = v2->next;
+ vol = VolumeSign(f0, v3);
+ while (!vol) {
+ if ((v3 = v3->next) == v0) {
+ if (debug) {
+ printf("DoubleTriangle: All points are coplanar!\n");
+ }
+ return (-2);
+ }
+ vol = VolumeSign(f0, v3);
+ }
+
+ /* Insure that v3 will be the first added. */
+ vertices = v3;
+ if (debug) {
+ fprintf(stderr,
+ "DoubleTriangle: finished. Head repositioned at v3.\n");
+ PrintOut(vertices);
+ }
+
+
+ return (0);
+}
+
+
+/*---------------------------------------------------------------------
+ConstructHull adds the vertices to the hull one at a time. The hull
+vertices are those in the list marked as onhull.
+---------------------------------------------------------------------*/
+void ConstructHull(void)
+{
+ tVertex v, vnext;
+ long int vol;
+ bool changed; /* T if addition changes hull; not used. */
+ int i;
+ int numVertices;
+
+ if (VERBOSE) {
+ fprintf(stdout, " Constructing 3D hull: \n");
+ }
+
+
+ v = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ v = vnext;
+ i++;
+ } while (v != vertices);
+ numVertices = i;
+
+ v = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ if (!v->mark) {
+ v->mark = PROCESSED;
+ changed = AddOne(v);
+ CleanUp();
+
+ if (check) {
+ fprintf(stderr, "ConstructHull: After Add of %d & Cleanup:\n",
+ v->vnum);
+ Checks();
+ }
+ if (debug)
+ PrintOut(v);
+ }
+ v = vnext;
+ i++;
+ if (VERBOSE) {
+ G_percent(i, numVertices, 1);
+ }
+ } while (v != vertices);
+
+ fflush(stdout);
+
+}
+
+/*---------------------------------------------------------------------
+AddOne is passed a vertex. It first determines all faces visible from
+that point. If none are visible then the point is marked as not
+onhull. Next is a loop over edges. If both faces adjacent to an edge
+are visible, then the edge is marked for deletion. If just one of the
+adjacent faces is visible then a new face is constructed.
+---------------------------------------------------------------------*/
+bool AddOne(tVertex p)
+{
+ tFace f;
+ tEdge e, temp;
+ long int vol;
+ bool vis = BFALSE;
+
+ if (debug) {
+ fprintf(stderr, "AddOne: starting to add v%d.\n", p->vnum);
+ PrintOut(vertices);
+ }
+
+ /* Mark faces visible from p. */
+ f = faces;
+ do {
+ vol = VolumeSign(f, p);
+ if (debug)
+ fprintf(stderr,
+ "faddr: %6x paddr: %6x Vol = %li\n", (unsigned int)f,
+ (unsigned int)p, vol);
+ if (vol < 0) {
+ f->visible = VISIBLE;
+ vis = BTRUE;
+ }
+ f = f->next;
+ } while (f != faces);
+
+ /* If no faces are visible from p, then p is inside the hull. */
+ if (!vis) {
+ p->onhull = !ONHULL;
+ return BFALSE;
+ }
+
+ /* Mark edges in interior of visible region for deletion.
+ Erect a newface based on each border edge. */
+ e = edges;
+ do {
+ temp = e->next;
+ if (e->adjface[0]->visible && e->adjface[1]->visible)
+ /* e interior: mark for deletion. */
+ e->delete = REMOVED;
+ else if (e->adjface[0]->visible || e->adjface[1]->visible)
+ /* e border: make a new face. */
+ e->newface = MakeConeFace(e, p);
+ e = temp;
+ } while (e != edges);
+ return BTRUE;
+}
+
+/*---------------------------------------------------------------------
+VolumeSign returns the sign of the volume of the tetrahedron determined by f
+and p. VolumeSign is +1 iff p is on the negative side of f,
+where the positive side is determined by the rh-rule. So the volume
+is positive if the ccw normal to f points outside the tetrahedron.
+The final fewer-multiplications form is due to Bob Williamson.
+---------------------------------------------------------------------*/
+int VolumeSign(tFace f, tVertex p)
+{
+ double vol;
+ long int voli;
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+
+ ax = f->vertex[0]->v[X] - p->v[X];
+ ay = f->vertex[0]->v[Y] - p->v[Y];
+ az = f->vertex[0]->v[Z] - p->v[Z];
+ bx = f->vertex[1]->v[X] - p->v[X];
+ by = f->vertex[1]->v[Y] - p->v[Y];
+ bz = f->vertex[1]->v[Z] - p->v[Z];
+ cx = f->vertex[2]->v[X] - p->v[X];
+ cy = f->vertex[2]->v[Y] - p->v[Y];
+ cz = f->vertex[2]->v[Z] - p->v[Z];
+
+ vol = ax * (by * cz - bz * cy)
+ + ay * (bz * cx - bx * cz)
+ + az * (bx * cy - by * cx);
+
+ if (debug)
+ fprintf(stderr,
+ "Face=%6x; Vertex=%d: vol(int) = %li, vol(double) = %.f\n",
+ (unsigned int)f, p->vnum, voli, vol);
+
+ /* The volume should be an integer. */
+ if (vol > 0.5)
+ return 1;
+ else if (vol < -0.5)
+ return -1;
+ else
+ return 0;
+}
+
+/*---------------------------------------------------------------------*/
+int Volumei(tFace f, tVertex p)
+{
+ double vol;
+ long int voli;
+ double ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz;
+
+ ax = f->vertex[0]->v[X] - p->v[X];
+ ay = f->vertex[0]->v[Y] - p->v[Y];
+ az = f->vertex[0]->v[Z] - p->v[Z];
+ bx = f->vertex[1]->v[X] - p->v[X];
+ by = f->vertex[1]->v[Y] - p->v[Y];
+ bz = f->vertex[1]->v[Z] - p->v[Z];
+ cx = f->vertex[2]->v[X] - p->v[X];
+ cy = f->vertex[2]->v[Y] - p->v[Y];
+ cz = f->vertex[2]->v[Z] - p->v[Z];
+
+ vol = (ax * (by * cz - bz * cy)
+ + ay * (bz * cx - bx * cz)
+ + az * (bx * cy - by * cx));
+
+ if (debug)
+ fprintf(stderr,
+ "Face=%6x; Vertex=%d: vol(int) = %li, vol(double) = %.f\n",
+ (unsigned int)f, p->vnum, voli, vol);
+
+ /* The volume should be an integer. */
+ if (vol > 0.5)
+ return 1;
+ else if (vol < -0.5)
+ return -1;
+ else
+ return 0;
+}
+
+
+/*-------------------------------------------------------------------*/
+void PrintPoint(tVertex p)
+{
+ int i;
+
+ for (i = 0; i < 3; i++)
+ printf("\t%li", p->v[i]);
+ putchar('\n');
+}
+
+/*---------------------------------------------------------------------
+MakeConeFace makes a new face and two new edges between the
+edge and the point that are passed to it. It returns a pointer to
+the new face.
+---------------------------------------------------------------------*/
+tFace MakeConeFace(tEdge e, tVertex p)
+{
+ tEdge new_edge[2];
+ tFace new_face;
+ int i, j;
+
+ /* Make two new edges (if don't already exist). */
+ for (i = 0; i < 2; ++i)
+ /* If the edge exists, copy it into new_edge. */
+ if (!(new_edge[i] = e->endpts[i]->duplicate)) {
+ /* Otherwise (duplicate is NULL), MakeNullEdge. */
+ new_edge[i] = MakeNullEdge();
+ new_edge[i]->endpts[0] = e->endpts[i];
+ new_edge[i]->endpts[1] = p;
+ e->endpts[i]->duplicate = new_edge[i];
+ }
+
+ /* Make the new face. */
+ new_face = MakeNullFace();
+ new_face->edge[0] = e;
+ new_face->edge[1] = new_edge[0];
+ new_face->edge[2] = new_edge[1];
+ MakeCcw(new_face, e, p);
+
+ /* Set the adjacent face pointers. */
+ for (i = 0; i < 2; ++i)
+ for (j = 0; j < 2; ++j)
+ /* Only one NULL link should be set to new_face. */
+ if (!new_edge[i]->adjface[j]) {
+ new_edge[i]->adjface[j] = new_face;
+ break;
+ }
+
+ return new_face;
+}
+
+/*---------------------------------------------------------------------
+MakeCcw puts the vertices in the face structure in counterclock wise
+order. We want to store the vertices in the same
+order as in the visible face. The third vertex is always p.
+---------------------------------------------------------------------*/
+void MakeCcw(tFace f, tEdge e, tVertex p)
+{
+ tFace fv; /* The visible face adjacent to e */
+ int i; /* Index of e->endpoint[0] in fv. */
+ tEdge s; /* Temporary, for swapping */
+
+ if (e->adjface[0]->visible)
+ fv = e->adjface[0];
+ else
+ fv = e->adjface[1];
+
+ /* Set vertex[0] & [1] of f to have the same orientation
+ as do the corresponding vertices of fv. */
+ for (i = 0; fv->vertex[i] != e->endpts[0]; ++i) ;
+ /* Orient f the same as fv. */
+ if (fv->vertex[(i + 1) % 3] != e->endpts[1]) {
+ f->vertex[0] = e->endpts[1];
+ f->vertex[1] = e->endpts[0];
+ }
+ else {
+ f->vertex[0] = e->endpts[0];
+ f->vertex[1] = e->endpts[1];
+ SWAP(s, f->edge[1], f->edge[2]);
+ }
+ /* This swap is tricky. e is edge[0]. edge[1] is based on endpt[0],
+ edge[2] on endpt[1]. So if e is oriented "forwards," we
+ need to move edge[1] to follow [0], because it precedes. */
+
+ f->vertex[2] = p;
+}
+
+/*---------------------------------------------------------------------
+MakeNullEdge creates a new cell and initializes all pointers to NULL
+and sets all flags to off. It returns a pointer to the empty cell.
+---------------------------------------------------------------------*/
+tEdge MakeNullEdge(void)
+{
+ tEdge e;
+
+ NEW(e, tsEdge);
+ e->adjface[0] = e->adjface[1] = e->newface = NULL;
+ e->endpts[0] = e->endpts[1] = NULL;
+ e->delete = !REMOVED;
+ ADD(edges, e);
+ return e;
+}
+
+/*--------------------------------------------------------------------
+MakeNullFace creates a new face structure and initializes all of its
+flags to NULL and sets all the flags to off. It returns a pointer
+to the empty cell.
+---------------------------------------------------------------------*/
+tFace MakeNullFace(void)
+{
+ tFace f;
+ int i;
+
+ NEW(f, tsFace);
+ for (i = 0; i < 3; ++i) {
+ f->edge[i] = NULL;
+ f->vertex[i] = NULL;
+ }
+ f->visible = !VISIBLE;
+ ADD(faces, f);
+ return f;
+}
+
+/*---------------------------------------------------------------------
+MakeFace creates a new face structure from three vertices (in ccw
+order). It returns a pointer to the face.
+---------------------------------------------------------------------*/
+tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace fold)
+{
+ tFace f;
+ tEdge e0, e1, e2;
+
+ /* Create edges of the initial triangle. */
+ if (!fold) {
+ e0 = MakeNullEdge();
+ e1 = MakeNullEdge();
+ e2 = MakeNullEdge();
+ }
+ else { /* Copy from fold, in reverse order. */
+ e0 = fold->edge[2];
+ e1 = fold->edge[1];
+ e2 = fold->edge[0];
+ }
+ e0->endpts[0] = v0;
+ e0->endpts[1] = v1;
+ e1->endpts[0] = v1;
+ e1->endpts[1] = v2;
+ e2->endpts[0] = v2;
+ e2->endpts[1] = v0;
+
+ /* Create face for triangle. */
+ f = MakeNullFace();
+ f->edge[0] = e0;
+ f->edge[1] = e1;
+ f->edge[2] = e2;
+ f->vertex[0] = v0;
+ f->vertex[1] = v1;
+ f->vertex[2] = v2;
+
+ /* Link edges to face. */
+ e0->adjface[0] = e1->adjface[0] = e2->adjface[0] = f;
+
+ return f;
+}
+
+/*---------------------------------------------------------------------
+CleanUp goes through each data structure list and clears all
+flags and NULLs out some pointers. The order of processing
+(edges, faces, vertices) is important.
+---------------------------------------------------------------------*/
+void CleanUp(void)
+{
+ CleanEdges();
+ CleanFaces();
+ CleanVertices();
+}
+
+/*---------------------------------------------------------------------
+CleanEdges runs through the edge list and cleans up the structure.
+If there is a newface then it will put that face in place of the
+visible face and NULL out newface. It also deletes so marked edges.
+---------------------------------------------------------------------*/
+void CleanEdges(void)
+{
+ tEdge e; /* Primary index into edge list. */
+ tEdge t; /* Temporary edge pointer. */
+
+ /* Integrate the newface's into the data structure. */
+ /* Check every edge. */
+ e = edges;
+ do {
+ if (e->newface) {
+ if (e->adjface[0]->visible)
+ e->adjface[0] = e->newface;
+ else
+ e->adjface[1] = e->newface;
+ e->newface = NULL;
+ }
+ e = e->next;
+ } while (e != edges);
+
+ /* Delete any edges marked for deletion. */
+ while (edges && edges->delete) {
+ e = edges;
+ DELETE(edges, e);
+ }
+ e = edges->next;
+ do {
+ if (e->delete) {
+ t = e;
+ e = e->next;
+ DELETE(edges, t);
+ }
+ else
+ e = e->next;
+ } while (e != edges);
+}
+
+/*---------------------------------------------------------------------
+CleanFaces runs through the face list and deletes any face marked visible.
+---------------------------------------------------------------------*/
+void CleanFaces(void)
+{
+ tFace f; /* Primary pointer into face list. */
+ tFace t; /* Temporary pointer, for deleting. */
+
+
+ while (faces && faces->visible) {
+ f = faces;
+ DELETE(faces, f);
+ }
+ f = faces->next;
+ do {
+ if (f->visible) {
+ t = f;
+ f = f->next;
+ DELETE(faces, t);
+ }
+ else
+ f = f->next;
+ } while (f != faces);
+}
+
+/*---------------------------------------------------------------------
+CleanVertices runs through the vertex list and deletes the
+vertices that are marked as processed but are not incident to any
+undeleted edges.
+---------------------------------------------------------------------*/
+void CleanVertices(void)
+{
+ tEdge e;
+ tVertex v, t;
+
+ /* Mark all vertices incident to some undeleted edge as on the hull. */
+ e = edges;
+ do {
+ e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
+ e = e->next;
+ } while (e != edges);
+
+ /* Delete all vertices that have been processed but
+ are not on the hull. */
+ while (vertices && vertices->mark && !vertices->onhull) {
+ v = vertices;
+ DELETE(vertices, v);
+ }
+ v = vertices->next;
+ do {
+ if (v->mark && !v->onhull) {
+ t = v;
+ v = v->next;
+ DELETE(vertices, t)
+ }
+ else
+ v = v->next;
+ } while (v != vertices);
+
+ /* Reset flags. */
+ v = vertices;
+ do {
+ v->duplicate = NULL;
+ v->onhull = !ONHULL;
+ v = v->next;
+ } while (v != vertices);
+}
+
+/*---------------------------------------------------------------------
+Collinear checks to see if the three points given are collinear,
+by checking to see if each element of the cross product is zero.
+---------------------------------------------------------------------*/
+bool Collinear(tVertex a, tVertex b, tVertex c)
+{
+ return
+ (c->v[Z] - a->v[Z]) * (b->v[Y] - a->v[Y]) -
+ (b->v[Z] - a->v[Z]) * (c->v[Y] - a->v[Y]) == 0
+ && (b->v[Z] - a->v[Z]) * (c->v[X] - a->v[X]) -
+ (b->v[X] - a->v[X]) * (c->v[Z] - a->v[Z]) == 0
+ && (b->v[X] - a->v[X]) * (c->v[Y] - a->v[Y]) -
+ (b->v[Y] - a->v[Y]) * (c->v[X] - a->v[X]) == 0;
+}
+
+/*---------------------------------------------------------------------
+Consistency runs through the edge list and checks that all
+adjacent faces have their endpoints in opposite order. This verifies
+that the vertices are in counterclockwise order.
+---------------------------------------------------------------------*/
+void Consistency(void)
+{
+ register tEdge e;
+ register int i, j;
+
+ e = edges;
+
+ do {
+ /* find index of endpoint[0] in adjacent face[0] */
+ for (i = 0; e->adjface[0]->vertex[i] != e->endpts[0]; ++i) ;
+
+ /* find index of endpoint[0] in adjacent face[1] */
+ for (j = 0; e->adjface[1]->vertex[j] != e->endpts[0]; ++j) ;
+
+ /* check if the endpoints occur in opposite order */
+ if (!(e->adjface[0]->vertex[(i + 1) % 3] ==
+ e->adjface[1]->vertex[(j + 2) % 3] ||
+ e->adjface[0]->vertex[(i + 2) % 3] ==
+ e->adjface[1]->vertex[(j + 1) % 3]))
+ break;
+ e = e->next;
+
+ } while (e != edges);
+
+ if (e != edges)
+ fprintf(stderr, "Checks: edges are NOT consistent.\n");
+ else
+ fprintf(stderr, "Checks: edges consistent.\n");
+
+}
+
+/*---------------------------------------------------------------------
+Convexity checks that the volume between every face and every
+point is negative. This shows that each point is inside every face
+and therefore the hull is convex.
+---------------------------------------------------------------------*/
+void Convexity(void)
+{
+ register tFace f;
+ register tVertex v;
+ long int vol;
+
+ f = faces;
+
+ do {
+ v = vertices;
+ do {
+ if (v->mark) {
+ vol = VolumeSign(f, v);
+ if (vol < 0)
+ break;
+ }
+ v = v->next;
+ } while (v != vertices);
+
+ f = f->next;
+
+ } while (f != faces);
+
+ if (f != faces)
+ fprintf(stderr, "Checks: NOT convex.\n");
+ else if (check)
+ fprintf(stderr, "Checks: convex.\n");
+}
+
+/*---------------------------------------------------------------------
+CheckEuler checks Euler's relation, as well as its implications when
+all faces are known to be triangles. Only prints positive information
+when debug is true, but always prints negative information.
+---------------------------------------------------------------------*/
+void CheckEuler(long int V, long int E, long int F)
+{
+ if (check)
+ fprintf(stderr, "Checks: V, E, F = %li %li %li:\t", V, E, F);
+
+ if ((V - E + F) != 2)
+ fprintf(stderr, "Checks: V-E+F != 2\n");
+ else if (check)
+ fprintf(stderr, "V-E+F = 2\t");
+
+
+ if (F != (2 * V - 4))
+ fprintf(stderr, "Checks: F=%li != 2V-4=%li; V=%li\n",
+ F, 2 * V - 4, V);
+ else if (check)
+ fprintf(stderr, "F = 2V-4\t");
+
+ if ((2 * E) != (3 * F))
+ fprintf(stderr, "Checks: 2E=%li != 3F=%li; E=%li, F=%li\n",
+ 2 * E, 3 * F, E, F);
+ else if (check)
+ fprintf(stderr, "2E = 3F\n");
+}
+
+/*-------------------------------------------------------------------*/
+void Checks(void)
+{
+ tVertex v;
+ tEdge e;
+ tFace f;
+ long int V = 0, E = 0, F = 0;
+
+ Consistency();
+ Convexity();
+ if (v = vertices)
+ do {
+ if (v->mark)
+ V++;
+ v = v->next;
+ } while (v != vertices);
+ if (e = edges)
+ do {
+ E++;
+ e = e->next;
+ } while (e != edges);
+ if (f = faces)
+ do {
+ F++;
+ f = f->next;
+ } while (f != faces);
+ CheckEuler(V, E, F);
+}
+
+
+/*===================================================================
+These functions are used whenever the debug flag is set.
+They print out the entire contents of each data structure.
+Printing is to standard error. To grab the output in a file in the csh,
+use this:
+ chull < i.file >&! o.file
+=====================================================================*/
+
+/*-------------------------------------------------------------------*/
+void PrintOut(tVertex v)
+{
+ fprintf(stderr, "\nHead vertex %d = %6x :\n", v->vnum, (unsigned int)v);
+ PrintVertices();
+ PrintEdges();
+ PrintFaces();
+}
+
+/*-------------------------------------------------------------------*/
+void PrintVertices(void)
+{
+ tVertex temp;
+
+ temp = vertices;
+ fprintf(stderr, "Vertex List\n");
+ if (vertices)
+ do {
+ fprintf(stderr, " addr %6x\t", (unsigned int)vertices);
+ fprintf(stderr, " vnum %4d", vertices->vnum);
+ fprintf(stderr, " (%6li,%6li,%6li)", vertices->v[X],
+ vertices->v[Y], vertices->v[Z]);
+ fprintf(stderr, " active:%3d", vertices->onhull);
+ fprintf(stderr, " dup:%5x", (unsigned int)vertices->duplicate);
+ fprintf(stderr, " mark:%2d\n", vertices->mark);
+ vertices = vertices->next;
+ } while (vertices != temp);
+
+}
+
+/*-------------------------------------------------------------------*/
+void PrintEdges(void)
+{
+ tEdge temp;
+ int i;
+
+ temp = edges;
+ fprintf(stderr, "Edge List\n");
+ if (edges)
+ do {
+ fprintf(stderr, " addr: %6x\t", (unsigned int)edges);
+ fprintf(stderr, "adj: ");
+ for (i = 0; i < 2; ++i)
+ fprintf(stderr, "%6x", (unsigned int)edges->adjface[i]);
+ fprintf(stderr, " endpts:");
+ for (i = 0; i < 2; ++i)
+ fprintf(stderr, "%4d", edges->endpts[i]->vnum);
+ fprintf(stderr, " del:%3d\n", edges->delete);
+ edges = edges->next;
+ } while (edges != temp);
+
+}
+
+/*-------------------------------------------------------------------*/
+void PrintFaces(void)
+{
+ int i;
+ tFace temp;
+
+ temp = faces;
+ fprintf(stderr, "Face List\n");
+ if (faces)
+ do {
+ fprintf(stderr, " addr: %6x\t", (unsigned int)faces);
+ fprintf(stderr, " edges:");
+ for (i = 0; i < 3; ++i)
+ fprintf(stderr, "%6x", (unsigned int)faces->edge[i]);
+ fprintf(stderr, " vert:");
+ for (i = 0; i < 3; ++i)
+ fprintf(stderr, "%4d", faces->vertex[i]->vnum);
+ fprintf(stderr, " vis: %d\n", faces->visible);
+ faces = faces->next;
+ } while (faces != temp);
+
+}
Property changes on: grass-addons/grass7/raster/r.vol.dem/chull.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.vol.dem/chull.h
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/chull.h (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/chull.h 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,2 @@
+int make3DHull(long int *X, long int *Y, long int *Z, int num_points,
+ FILE * tmpfile);
Property changes on: grass-addons/grass7/raster/r.vol.dem/chull.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.vol.dem/globals.h
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/globals.h (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/globals.h 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,40 @@
+#ifndef GLOBAL_H
+#define GLOBAL_H
+
+#define PROGVERSION 0.55
+#define PROGNAME "r.vol.dem"
+
+/* if DEBUG > 1 we will dump even more details */
+#define DEBUG 0
+
+/* NULL valued cells */
+#define NULLVALUE -1
+/* masked cells */
+#define MASKVALUE -2
+/* unset cells */
+#define UNSET -3
+
+/* value to represent NULL in VTK files */
+#define DNULLVALUE -99999.99
+
+#define TMPFILE "voxeltmp.tmp"
+
+#ifdef MAIN
+/* number of decimal places with which coordinates are stored */
+double PRECISION;
+
+/* verbose output includes progress display */
+int VERBOSE;
+
+/* number of DEMs in the input */
+int NSLICES;
+
+#else
+
+extern double PRECISION;
+extern int VERBOSE;
+extern int NSLICES;
+
+#endif /* MAIN */
+
+#endif /* GLOBAL_H */
Property changes on: grass-addons/grass7/raster/r.vol.dem/globals.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.vol.dem/macros.h
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/macros.h (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/macros.h 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,39 @@
+
+/*====================================================================
+ macros.h
+
+ macros used to access data structures and perform quick tests.
+
+ ====================================================================*/
+
+/* general-purpose macros */
+#define SWAP(t,x,y) { t = x; x = y; y = t; }
+
+#define NEW(p,type) if ((p=(type *) malloc (sizeof(type))) == NULL) {\
+ printf ("Out of Memory!\n");\
+ exit(0);\
+ }
+
+#define FREE(p) if (p) { free ((char *) p); p = NULL; }
+
+
+#define ADD( head, p ) if ( head ) { \
+ p->next = head; \
+ p->prev = head->prev; \
+ head->prev = p; \
+ p->prev->next = p; \
+ } \
+ else { \
+ head = p; \
+ head->next = head->prev = p; \
+ }
+
+#define DELETE( head, p ) if ( head ) { \
+ if ( head == head->next ) \
+ head = NULL; \
+ else if ( p == head ) \
+ head = head->next; \
+ p->next->prev = p->prev; \
+ p->prev->next = p->next; \
+ FREE( p ); \
+ }
Property changes on: grass-addons/grass7/raster/r.vol.dem/macros.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.vol.dem/main.c
===================================================================
--- grass-addons/grass7/raster/r.vol.dem/main.c (rev 0)
+++ grass-addons/grass7/raster/r.vol.dem/main.c 2012-04-24 21:23:22 UTC (rev 51526)
@@ -0,0 +1,1445 @@
+/* r.vol.dem */
+/* Purpose: Interpolate a voxel model from a series of DEMs by flood filling the voxel space in between. */
+
+
+/* TODO:
+
+ BUG: output_ascii_points (): when option -f is given, the v.in.ascii module aborts with an error message!
+ strangely, if you copy and paste the line from the tmpfile, which is exactly the same syntax,
+ it works! (???)
+
+ - fix warnings and license issues for chull.c (provide complete copy of the original sources)
+
+ - tempfiles
+ this program leaves lots of them around right now ...
+
+ - cell counts: AVOID counting cells twice for layers that have the same label!
+
+ - USER MAN PAGE:
+
+ filling up or down does not make a difference in volume, unless -f option is also given!
+ BUT: it makes a difference in the labeling of layers:
+ When filling down, the bottom DEM
+ is interpreted as the last boundary of a layer and is thus assigned an individual label
+ When fillin up, the top DEM
+ is interpreted as the last boundary of a layer and is thus assigned an individual label
+ Where this is undesirable, boundaries can be dissolved by giving them identical labels using the labels= option!
+
+ Visible differences in Paraview are due to rounding effects in the software. Displaying points in NVIZ will show you
+ that both data sets are the same [check if this ist true: same number of points in both cases?]
+
+ Vector hulls are currently very rough approximations with lots of shortcomings
+
+ Layers are numbered from bottom to top : 0 .. n
+
+ DEMs have to be given from bottom to top as well. Order matters!
+
+ - LICENSE ISSUES with convex hull 3D code !(?): include at least a full version of original chull.c and add some comments!
+ - also replace warinng messages (e.g. all points colinear) with own versions.
+
+ - make clean and get rid of compiler warnings
+
+ in a first pass, that can also be disabled by the user
+ - check for overlaps
+
+ in a second pass, to be optionally enabled by user:
+ - create a 0/1 mask file for each DEM to mask
+ out interpolation for edges that are not straight
+
+ - fuzzy layer boundaries (!)
+
+ - in mode -f revert interpolation direction once at the end to make sure everything gets filled (?)!
+ - maybe as a separate option -r if it will be good for anything
+
+ - standard options for voxels (data type ...)
+
+ - automatically adjust 3D region to fit full DEM data (x,y,z)
+
+ - output a map that shows areas with topologic problems. User can overlay all DEMs and this map and query for problems.
+ Map will contain highest ID of the DEM(s)
+ that caused a problem here. After clearing all problems, user might have to generate another map with errors from lower ID DEMs ...
+ [DONE BUT UNTESTED]
+
+ - DNULLVALUE needs to be adjustable, also many other vars in globals.h: PRECISION, TMPFILE, ...
+
+ - EXPORTS: switch all export functions away from ASCII bridge:
+ 3D vector hulls output:
+ - export 3D vector hulls into one GRASS bin file and create proper attribute table
+ - 3D vector hulls don't overlap nicely, because the points are in the voxel centers.
+ all hull points have to be extruded half a voxel size into normal direction!
+ - layers with all points on one voxel plane do not get turned into a hull currently!
+ - RRR:GGG:BBB color option for each layer
+ voxel output:
+ - add history information
+ - a category number for each DEM (3d raster cats?)
+ - RRR:GGG:BBB color option for each layer
+
+ - surface smoothing for VTK file output
+ smoothing surfaces for VTK-output with top= and bottom= options (r3.out.vtk) smoothes all surfaces by degrees!
+
+ - more progress indicators in export functions! [lots done but test what it looks like]
+
+ - implement quite operation: most messages don't respect VERBOSE! [lots done, but needs testing]
+
+ - v.in.ascii now takes --o for overwrite. Since when?
+
+ */
+
+#define MAIN
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include <string.h>
+#include <math.h>
+#include <time.h>
+
+#include <grass/gis.h>
+#include <grass/raster3d.h>
+#include <grass/glocale.h>
+
+#include "chull.h"
+
+#include "globals.h"
+
+/* ================================================================================
+
+ GLOBAL DECLARATIONS
+
+ =================================================================================== */
+
+
+/* by making this global, all func's will have access to module options and flags */
+struct
+{
+ struct Option *maps; /* raster DEM(s) to work on */
+ struct Option *output; /* output rast3d map */
+ struct Option *values; /* specify the value that each layer represents */
+ struct Option *colors; /* RRR:GGG:BBB triplets for each layer */
+ struct Option *error; /* raster map to represent topology errors in input DEMs */
+ struct Option *algorithm; /* choose interpolation algorithm */
+}
+parm;
+struct
+{
+ struct Flag *countcells; /* calculate 3D cell stats for each layer */
+ struct Flag *fillnull; /* fill NULL areas in input DEMs from above or below */
+ struct Flag *hull; /* export a 3D hull polygon for every layer */
+ struct Flag *grasspts; /* produce a GRASS points vector file */
+ struct Flag *points; /* produce VTK point data instead of cells */
+ struct Flag *quiet; /* no status display on screen */
+ struct Flag *skip; /* start interpolation even if topology errors were detected */
+ struct Flag *vtk; /* generate a VTK file */
+ struct Flag *zadjust; /* adjust z range of current region to fit DEM data */
+}
+flag;
+
+
+
+/* ================================================================================
+
+ UTILITY FUNCTIONS
+
+ =================================================================================== */
+
+/* returns 1 if the integer is a DEM ID, not a MASKED or NULL value or UNSET,
+ 0 otherwise */
+int isValue(int val)
+{
+ if (val == UNSET) {
+ return (0);
+ }
+ if (val == NULLVALUE) {
+ return (0);
+ }
+ if (val == MASKVALUE) {
+ return (0);
+ }
+ return (1);
+}
+
+int isUnset(int val)
+{
+ if (val == UNSET) {
+ return (1);
+ }
+ return (0);
+}
+
+int isMask(int val)
+{
+ if (val == MASKVALUE) {
+ return (1);
+ }
+ return (0);
+}
+
+int isNull(int val)
+{
+ if (val == NULLVALUE) {
+ return (1);
+ }
+ return (0);
+}
+
+
+/* returns 1 if the string contains a valid RRR:GGG:BBB color triplet,
+ 0 otherwise
+
+ IF the string is valid, the R, G and B levels (0-255) will be
+ stored at the int pointers (NULL otherwise)
+ */
+
+int isRGB(char *triplet, unsigned short int *R, unsigned short int *G,
+ unsigned short int *B)
+{
+ R = NULL;
+ G = NULL;
+ B = NULL;
+ return (0);
+}
+
+
+/* ================================================================================
+
+ Z-RANGE AND TOPOLOGY CHECKING
+
+ =================================================================================== */
+
+
+/* checks if a DEM value intersects a voxel at a given depth */
+/* depth goes from 0 (bottom voxel) to slices-1 (top voxel) */
+/* returns 1 if DEM intersects, 0 if not */
+
+/* if a DEM intersects the top surface of a voxel, that counts as */
+/* a cut, the bottom does not. */
+int cut(DCELL val, int depth, struct Cell_head window)
+{
+ if (DEBUG > 2) {
+ fprintf(stdout, "CUT: %.3f at depth %i", (double)val, depth);
+ }
+ if ((val > (window.bottom + depth * window.tb_res)) &&
+ (val < (window.bottom + (depth + 1) * window.tb_res))) {
+ if (DEBUG > 2) {
+ fprintf(stdout, " = YES\n");
+ }
+ return (1);
+ }
+ if (val == window.bottom + (depth + 1) * window.tb_res) {
+ if (DEBUG > 1) {
+ fprintf(stdout, " = YES\n");
+ }
+ return (1);
+ }
+ if (DEBUG > 2) {
+ fprintf(stdout, " = NO\n");
+ }
+ return (0);
+}
+
+
+/* Adjust the ACTIVE region's Z range so that it will be big enough for the
+ highest and lowest point in the input DEMs.
+ This setting will stay in effect even after r.vol.dem exits!
+ */
+void adjust_zrange(void)
+{
+ fprintf(stderr, "DEBUG: Adjustment of Z range not implemented yet.\n");
+ return;
+}
+
+
+/* check if DEMs are within z-range of current 3d region */
+void check_zrange(int *fd, DCELL ** dcell, int nfiles)
+{
+
+ struct Cell_head window;
+
+ int i;
+ int row, col, nrows, ncols;
+
+ double top, bottom;
+ int topDEM, bottomDEM;
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ G_get_window(&window);
+
+ if (flag.zadjust->answer) {
+ adjust_zrange();
+ fprintf(stderr, "DEBUG: Insert a return() statement here.\n");
+ }
+
+ if (VERBOSE)
+ fprintf(stdout, " Z-range checking, pass 1: \n");
+
+ /* get maximum and minimum elevation from all DEMs */
+ row = 0;
+ if (Rast_get_d_row(fd[0], dcell[0], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."), row,
+ parm.maps->answers[0]);
+ }
+ top = dcell[0][0];
+ topDEM = 0;
+ for (i = 0; i < nfiles; i++) {
+ for (row = 0; row < nrows; row++) {
+ if (Rast_get_d_row(fd[i], dcell[i], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."),
+ row, parm.maps->answers[i]);
+ }
+ for (col = 0; col < ncols; col++) {
+ if (dcell[i][col] > top) {
+ top = dcell[i][col];
+ topDEM = i;
+ }
+ }
+ }
+
+ G_percent(i, (nfiles - 1), 1);
+ }
+ fflush(stdout);
+ if (DEBUG)
+ fprintf(stdout, "TOP = %.2f in %s\n", top,
+ parm.maps->answers[topDEM]);
+
+ if (VERBOSE)
+ fprintf(stdout, " z-range-checking, pass 2: \n");
+
+ row = 0;
+ if (Rast_get_d_row(fd[0], dcell[0], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."), row,
+ parm.maps->answers[0]);
+ }
+ bottom = dcell[0][0];
+ bottomDEM = 0;
+ for (i = 0; i < nfiles; i++) {
+ for (row = 0; row < nrows; row++) {
+ if (Rast_get_d_row(fd[i], dcell[i], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."),
+ row, parm.maps->answers[i]);
+ }
+ for (col = 0; col < ncols; col++) {
+ if (dcell[i][col] < bottom) {
+ bottom = dcell[i][col];
+ bottomDEM = i;
+ }
+ }
+ }
+ G_percent(i, (nfiles - 1), 1);
+ }
+ fflush(stdout);
+ if (DEBUG)
+ fprintf(stdout, "BOTTOM = %.2f in %s\n", bottom,
+ parm.maps->answers[bottomDEM]);
+
+ if (top > window.top) {
+ G_fatal_error(_("Highest DEM value (%.3f in %s) outside extent of current 3d region (top=%.3f) "),
+ top, parm.maps->answers[topDEM], window.top);
+ }
+
+ if (bottom < window.bottom) {
+ G_fatal_error(_("Lowest DEM value (%.3f in %s) outside extent of current 3d region (bottom=%.3f) "),
+ bottom, parm.maps->answers[bottomDEM], window.bottom);
+ }
+
+ return;
+}
+
+
+void check_topologies(int *fd, DCELL ** dcell, int nfiles)
+{
+
+ int i, j;
+ struct Cell_head window;
+ int row, col, nrows, ncols;
+ int nerrors;
+
+ CELL **c_error = NULL;
+ int f_error = 0;
+
+ if (VERBOSE)
+ fprintf(stdout, " Checking DEM topologies: \n");
+
+ nerrors = 0;
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* does the user want a topology error map? */
+ /* if so, we need to alloc mem and create a new CELL map */
+ if (parm.error->answer != NULL) {
+
+ /* KILL ME */
+ fprintf(stderr,
+ "DEBUG: Topology error mapping is an untested feature.\n");
+ /* KILL ME */
+
+ c_error = (CELL **) G_malloc(nrows * sizeof(CELL *));
+ for (row = 0; row < nrows; row++) {
+ c_error[row] = Rast_allocate_c_buf();
+ }
+ f_error = Rast_open_new(parm.error->answer, CELL_TYPE);
+ if (f_error < 0) {
+ G_fatal_error
+ ("Could not create new CELL raster map for mapping topology errors.");
+ }
+ /* initialise error map */
+ for (row = 0; row < nrows; row++) {
+ Rast_set_c_null_value(&c_error[row][0], ncols);
+ }
+ }
+
+ for (row = 0; row < nrows; row++) {
+ for (i = 0; i < nfiles; i++) {
+ if (Rast_get_d_row(fd[i], dcell[i], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."),
+ row, parm.maps->answers[i]);
+ }
+ }
+ for (col = 0; col < ncols; col++) {
+ for (i = 0; i < nfiles; i++) {
+ /* check for undershoots from DEMs above */
+ for (j = i + 1; j < nfiles; j++) {
+ if (dcell[j][col] <= dcell[i][col]) {
+ /* UNDERSHOOT */
+
+ G_warning(_("Z value in %s undershoots Z value in %s at cells %i,%i, \ncoordinates: E=%f, N=%f.\n"),
+ parm.maps->answers[j],
+ parm.maps->answers[i], col, row,
+ Rast_col_to_easting((double)col, &window),
+ Rast_row_to_northing((double)row, &window)
+ );
+ nerrors++;
+ if (parm.error->answer != NULL) {
+ c_error[row][col] = j;
+ }
+ }
+ }
+ /* check for overshoots from DEMs below */
+ for (j = i - 1; j >= 0; j--) {
+ if (dcell[j][col] >= dcell[i][col]) {
+ /* OVERSHOOT! */
+ G_warning(_("Z value in %s overshoots Z value in %s at cells %i,%i, \ncoordinates: E=%f, N=%f.\n"),
+ parm.maps->answers[j],
+ parm.maps->answers[i], col, row,
+ Rast_col_to_easting((double)col, &window),
+ Rast_row_to_northing((double)row, &window)
+ );
+
+ nerrors++;
+ if (parm.error->answer != NULL) {
+ c_error[row][col] = j;
+ }
+ }
+ }
+ }
+ }
+
+ if (VERBOSE)
+ G_percent(row, (nrows - 1), 2);
+ }
+ fflush(stdout);
+ if (parm.error->answer != NULL) {
+ for (row = 0; row < nrows; row++) {
+ Rast_put_c_row(f_error, c_error[row]);
+ }
+ Rast_close(f_error);
+ }
+ if (nerrors > 0) {
+ G_fatal_error(_("Found %i errors in DEM topology! Aborting."),
+ nerrors);
+ }
+
+ return;
+}
+
+
+
+/* ================================================================================
+
+ CORE INTERPOLATION ROUTINES
+
+ =================================================================================== */
+
+
+void interpolate_up(int *flist, DCELL ** dcell, int col, int nfiles)
+{
+
+ struct Cell_head window;
+ int fillval;
+ double level;
+ double stopat;
+ int i, j, k;
+ int DEMAbove;
+ int highestDEMValue;
+
+ Rast_get_window(&window);
+
+ /* DEMs must be provided in this order: bottom,...,top ... */
+ for (i = 0; i < NSLICES; i++) {
+ /* ... the interpolated column, however, has its lowest value */
+ /* at the bottom (last element of array), so that order needs to be switched! */
+ flist[(NSLICES - 1) - i] = NULLVALUE; /* NULLVALUE = no DEM cuts at this point */
+ for (j = 0; j < nfiles; j++) {
+ if (!Rast_is_d_null_value(&dcell[j][col])) { /* SKIP NULL VALUED CELLS */
+ if (cut(dcell[j][col], i, window)) {
+ flist[(NSLICES - 1) - i] = j;
+ }
+ }
+ }
+ }
+
+ /* interpolate one column by 'flood filling' DEM values up from the bottom */
+ fillval = flist[0];
+
+ /* first we need to know where to stop the interpolation ! */
+
+ /* find the highest DEM that has a data value at this point! */
+ highestDEMValue = -1;
+ i = (nfiles - 1);
+ while ((i >= 0) && (highestDEMValue == -1)) {
+ if (!Rast_is_d_null_value(&dcell[i][col])) {
+ highestDEMValue = i;
+ }
+ i--;
+ }
+ i++;
+
+ if (highestDEMValue == -1) {
+ /* there is no DEM with any data at this point ! */
+ stopat = window.bottom - 1; /* this will prevent any interpolation from happening here ! */
+ }
+ else {
+ /* we found a DEM with a value ... */
+ if (highestDEMValue == (nfiles - 1)) {
+ /* ... and it is the top one, so we just set this as interpolation stop */
+ stopat = dcell[(nfiles - 1)][col];
+ }
+ else {
+ /* ... but it is not the top one, so we need to decide what to do */
+ if (flag.fillnull->answer) {
+ /* this flag forces us to keep filling (all the way up to the top of the current 3D region) */
+ stopat = window.top;
+ }
+ else {
+ /* we will only fill up to the highest DEM with a value */
+ /* but we make sure to fill up at least one slice ! */
+ stopat =
+ dcell[highestDEMValue][col] + (window.tb_res +
+ (window.tb_res * 0.5));
+ }
+ }
+ }
+
+ if (DEBUG > 1) {
+ fprintf(stderr, " STOPAT = %.2f\n", stopat);
+ }
+
+ for (i = (NSLICES - 1); i >= 0; i--) {
+ if (isValue(flist[i])) {
+ /* check, if there is another DEM cutting point above this slice */
+ /* if not, only keep filling, if flag -f has been set ! */
+ if (flag.fillnull->answer) {
+ fillval = flist[i];
+ }
+ else {
+ DEMAbove = 0;
+ for (k = i; k >= 0; k--) {
+ if (isValue(flist[k])) {
+ DEMAbove = 1;
+ }
+ }
+ if (DEMAbove) {
+ fillval = flist[i];
+ }
+ }
+ }
+ else {
+ /* check if we are still below top level */
+ level = window.bottom + (((NSLICES - 1) - i) * window.tb_res);
+ if (level < stopat) {
+ flist[i] = fillval;
+ }
+ else {
+ flist[i] = NULLVALUE;
+ }
+ }
+ }
+
+ return;
+
+}
+
+
+void interpolate_down(int *flist, DCELL ** dcell, int col, int nfiles)
+{
+
+ struct Cell_head window;
+ int fillval;
+ double level;
+ double stopat = 0.0;
+ int i, j, k;
+ int DEMBelow;
+ int lowestDEMValue;
+
+ Rast_get_window(&window);
+
+ /* DEMs must be provided in this order: bottom,...,top ... */
+ for (i = 0; i < NSLICES; i++) {
+ /* ... the interpolated column, however, has its lowest value */
+ /* at the bottom, so that order needs to be switched! */
+ flist[(NSLICES - 1) - i] = NULLVALUE; /* NULLVALUE = no DEM cuts at this point */
+ for (j = 0; j < nfiles; j++) {
+ if (!Rast_is_d_null_value(&dcell[j][col])) { /* SKIP NULL VALUED CELLS */
+ if (cut(dcell[j][col], i, window)) {
+ flist[(NSLICES - 1) - i] = j;
+ }
+ }
+ }
+ }
+
+ /* interpolate one column by 'flood filling' DEM values down from the top */
+ fillval = flist[0];
+
+ /* first we need to know where to stop the interpolation ! */
+
+ /* find the lowest DEM that has a data value at this point! */
+ lowestDEMValue = -1;
+ i = 0;
+ while ((i < nfiles) && (lowestDEMValue == -1)) {
+ if (!Rast_is_d_null_value(&dcell[i][col])) {
+ lowestDEMValue = i;
+ }
+ i++;
+ }
+ i--;
+
+ if (lowestDEMValue == -1) {
+ /* there is no DEM with any data at this point ! */
+ stopat = window.top + 1; /* this will prevent any interpolation from happening here ! */
+ }
+ else {
+ /* we found a DEM with a value ... */
+ if (lowestDEMValue == 0) {
+ /* ... and it is the bottom one, so we just set this as interpolation stop */
+ stopat = dcell[0][col];
+ }
+ else {
+ /* ... but it is not the bottom one, so we need to decide what to do */
+ if (flag.fillnull->answer) {
+ /* this flag forces us to keep filling (all the way down to the bottom of the current 3D region) */
+ stopat = window.bottom;
+ }
+ else {
+ /* we will only fill down to the lowest DEM with a value */
+ /* but we make sure to fill down at least one slice ! */
+ stopat =
+ dcell[lowestDEMValue][col] - (window.tb_res +
+ (window.tb_res * 0.5));
+ }
+ }
+ }
+
+ if (DEBUG > 1) {
+ fprintf(stderr, " STOPAT = %.2f\n", stopat);
+ }
+
+ for (i = 1; i < NSLICES; i++) {
+ if (isValue(flist[i])) {
+ /* check, if there is another DEM cutting point below this slice */
+ /* if not, only keep filling, if flag -f has been set ! */
+ if (flag.fillnull->answer) {
+ fillval = flist[i];
+ }
+ else {
+ DEMBelow = 0;
+ for (k = i; k <= NSLICES; k++) {
+ if (isValue(flist[k])) {
+ DEMBelow = 1;
+ }
+ }
+ if (DEMBelow) {
+ fillval = flist[i];
+ }
+ }
+ }
+ else {
+ /* check if we are still above bottom level */
+ level = window.bottom + (((NSLICES - 1) - i) * window.tb_res);
+ if (level >= stopat) {
+ flist[i] = fillval;
+ }
+ else {
+ flist[i] = NULLVALUE;
+ }
+ }
+ }
+
+ return;
+}
+
+
+void interpolate(int *fd, DCELL ** dcell, int nfiles, int nvalues,
+ double ***ascii)
+{
+
+ struct Cell_head window;
+
+ int *flist;
+
+ int i;
+ int row, col, nrows, ncols;
+
+ Rast_get_window(&window);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ flist = G_malloc(sizeof(int) * NSLICES);
+
+ if (VERBOSE)
+ fprintf(stdout, "\nInterpolating: \n");
+
+ for (row = 0; row < nrows; row++) {
+ /* read one row from all layers */
+ for (i = 0; i < nfiles; i++) {
+ if (Rast_get_d_row(fd[i], dcell[i], row) < 0) {
+ G_fatal_error(_("Could not read row number %i from map %s."),
+ row, parm.maps->answers[i]);
+ }
+ }
+ for (col = 0; col < ncols; col++) {
+
+ if (!strcmp(parm.algorithm->answer, "up")) {
+ interpolate_up(flist, dcell, col, nfiles);
+ }
+ if (!strcmp(parm.algorithm->answer, "down")) {
+ interpolate_down(flist, dcell, col, nfiles);
+ }
+
+ for (i = 1; i < NSLICES; i++) {
+ if (isValue(flist[i])) {
+ if (nvalues == 0) {
+ /* either enumerate layers from 0 .. n */
+ ascii[i][row][col] = (double)(flist[i]);
+ }
+ else {
+ /* OR write user-provided values for layers? */
+ ascii[i][row][col] =
+ atof(parm.values->answers[flist[i]]);
+ }
+ }
+ else {
+ /* write a NULL valued voxel */
+ ascii[i][row][col] = DNULLVALUE;
+ }
+ }
+ if (DEBUG > 1) {
+ fprintf(stdout, "col %i (", col);
+ for (i = 0; i < NSLICES; i++) {
+ if (isValue(flist[i])) {
+ fprintf(stdout, "%i ", flist[i]);
+ }
+ if (isNull(flist[i])) {
+ fprintf(stdout, "N ");
+ }
+ if (isMask(flist[i])) {
+ fprintf(stdout, "M ");
+ }
+ }
+ fprintf(stdout, ")\n");
+ }
+ }
+
+ if (VERBOSE) {
+ G_percent(row, nrows - 1, 2);
+ }
+ if (DEBUG > 1)
+ fprintf(stdout, "row=%i\n", row);
+ }
+ fprintf(stdout, "\n");
+ fflush(stdout);
+}
+
+
+/* ================================================================================
+
+ DATA OUTPUT AND EXPORT
+
+ =================================================================================== */
+
+void count_3d_cells(double ***ascii, int nfiles, int nvalues)
+{
+ struct Cell_head window;
+ int a, i, j, k;
+ unsigned long int sum;
+ int nrows, ncols;
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ fprintf(stdout, "\nCell counts: \n");
+
+ for (a = 0; a < nfiles; a++) {
+ sum = 0;
+ for (i = 0; i < NSLICES; i++) {
+ for (j = 0; j < nrows; j++) {
+ for (k = 0; k < ncols; k++) {
+ if (nvalues == 0) {
+ if (a == (int)ascii[i][j][k]) {
+ sum++;
+ }
+ }
+ else {
+ if (atof(parm.values->answers[a]) == ascii[i][j][k]) {
+ sum++;
+ }
+ }
+ }
+ }
+ }
+ if (nvalues == 0) {
+ fprintf(stdout, " Layer %i has %li 3D cells\n", a, sum);
+ }
+ else {
+ fprintf(stdout, " Layer %i (label=%f) has %li 3D cells\n", a,
+ atof(parm.values->answers[a]), sum);
+ }
+ }
+}
+
+
+void output_ascii(double ***ascii)
+{
+
+ struct Cell_head window;
+ int i, j, k;
+ int nrows, ncols;
+ FILE *tmpfile;
+
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* open a tmp file for writing ascii output data */
+ tmpfile = fopen(TMPFILE, "w+");
+ if (tmpfile == NULL) {
+ G_fatal_error(_("Could not create temporary file for writing ASCII output (%s)"),
+ TMPFILE);
+ }
+ /* write output to ASCII file */
+ fprintf(tmpfile, "north: %f\n", window.north);
+ fprintf(tmpfile, "south: %f\n", window.south);
+ fprintf(tmpfile, "east: %f\n", window.east);
+ fprintf(tmpfile, "west: %f\n", window.west);
+ fprintf(tmpfile, "top: %f\n", window.top);
+ fprintf(tmpfile, "bottom: %f\n", window.bottom);
+ fprintf(tmpfile, "rows: %i\n", nrows);
+ fprintf(tmpfile, "cols: %i\n", ncols);
+ fprintf(tmpfile, "levels:%i\n", NSLICES);
+ for (i = 0; i < NSLICES; i++) {
+ for (j = 0; j < nrows; j++) {
+ for (k = 0; k < ncols; k++) {
+ /* below is the correct data order for Paraview and GRASS rast3d ASCII files ! */
+ fprintf(tmpfile, "%f ",
+ ascii[(NSLICES - 1) - i][(nrows - 1) - j][k]);
+ }
+ fprintf(tmpfile, "\n");
+ }
+ }
+ fclose(tmpfile);
+}
+
+
+void output_ascii_points(double ***ascii)
+{
+
+ struct Cell_head window;
+ int i, j, k;
+ int nrows, ncols;
+ FILE *tmpfile;
+ char tmp[2048];
+ double x, y, z, w;
+ int error;
+
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* open a tmp file for writing ascii output data */
+ sprintf(tmp, "%s.points.txt", parm.output->answer);
+ tmpfile = fopen(tmp, "w+");
+ if (tmpfile == NULL) {
+ G_fatal_error(_("Could not create temporary file for writing ASCII output (%s)"),
+ tmp);
+ }
+
+ fprintf(tmpfile, "# Simple 3D vector points file.\n");
+ fprintf(tmpfile, "# This file was created by %s version %.2f \n",
+ PROGNAME, PROGVERSION);
+ fprintf(tmpfile, "# Tab delimited point data: x, y, z, value. \n");
+ fprintf(tmpfile, "# Import into GRASS GIS using this command:\n");
+ fprintf(tmpfile,
+ "# v.in.ascii in=%s.points.txt out=%s -z x=1 y=2 z=3 cat=0 columns='x double, y double, z double, w double' fs=tab\n",
+ parm.output->answer, parm.output->answer);
+
+ /* write output to ASCII file */
+ for (i = 0; i < NSLICES; i++) {
+ for (j = 0; j < nrows; j++) {
+ for (k = 0; k < ncols; k++) {
+ /* coordinates for points will be centers of 3D raster cells ! */
+ /* only produce points at non-NULL locations ! */
+ w = ascii[i][(nrows - 1) - j][k];
+ if (w != DNULLVALUE) {
+ x = window.west + (window.ew_res * k) +
+ (window.ew_res * 0.5);
+ y = window.south + (window.ns_res * j) +
+ (window.ns_res * 0.5);
+ z = window.top - (window.tb_res * i) -
+ (window.tb_res * 0.5);
+ fprintf(tmpfile, "%.6f\t%.6f\t%.6f\t%.6f\n", x, y, z, w);
+ }
+ }
+ }
+ }
+
+ if (VERBOSE) {
+ fprintf(stdout, "\nRunning v.in.ascii to create 3D points map: \n");
+ }
+
+ /* TODO: KILL THIS AND REPLACE WITH NEAT --o OPTION CHECK */
+ /* remove output map, if exists */
+ sprintf(tmp, "g.remove vect=%s", parm.output->answer);
+ error = system(tmp);
+ if (error != EXIT_SUCCESS) {
+ G_fatal_error("Error running command: %s", tmp);
+ }
+
+ /* call v.in.ascii to make GRASS native vector map */
+ sprintf(tmp,
+ "v.in.ascii in=%s.points.txt out=%s -z x=1 y=2 z=3 cat=0 columns='x double, y double, z double, w double' fs=tab\n",
+ parm.output->answer, parm.output->answer);
+
+ error = system(tmp);
+ if (error != EXIT_SUCCESS) {
+ G_fatal_error("Error running command: %s", tmp);
+ }
+
+ if (VERBOSE)
+ fprintf(stdout, " DONE \n");
+
+ fclose(tmpfile);
+}
+
+
+void output_ascii_hulls(double ***ascii, int nvalues, int nfiles)
+{
+
+ struct Cell_head window;
+ int i, j, k, l;
+ double label;
+ int num_points;
+ int nrows, ncols;
+ FILE *tmpfile;
+ char tmp[2048];
+ double w;
+ long int *px;
+ long int *py;
+ long int *pz;
+ int curPoint;
+
+ int error;
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+
+ /* MAIN LOOP START */
+ /* do this vor every layer, i.e. all 3D cells with the same label */
+ for (l = 0; l < nfiles; l++) {
+
+ /* open tmp file for writing ascii output data */
+ sprintf(tmp, "%s.hull.%i.txt", parm.output->answer, l);
+ tmpfile = fopen(tmp, "w+");
+ if (tmpfile == NULL) {
+ G_fatal_error(_("Could not create temporary file for writing ASCII output (%s)"),
+ tmp);
+ }
+
+ if (nvalues == 0) {
+ /* either enumerate layers from 0 .. (nfiles-1) */
+ label = (double)l;
+ }
+ else {
+ /* OR use user-provided values for layers? */
+ label = atof(parm.values->answers[l]);
+ }
+
+ num_points = 0;
+ /* create arrays large enough to store all points of the current layer */
+ for (i = 0; i < NSLICES; i++) {
+ for (j = 0; j < nrows; j++) {
+ for (k = 0; k < ncols; k++) {
+ w = ascii[i][j][k];
+ if (w == label) {
+ num_points++;
+ }
+ }
+ }
+ }
+
+ if (DEBUG) {
+ fprintf(stderr, "EXPORT HULLS: %i points with label %.6f\n",
+ num_points, label);
+ }
+
+ px = (long int *)malloc(num_points * sizeof(long int));
+ py = (long int *)malloc(num_points * sizeof(long int));
+ pz = (long int *)malloc(num_points * sizeof(long int));
+
+ /* pick points with current label from ASCII array */
+
+ curPoint = 0;
+ for (i = 0; i < NSLICES; i++) {
+ for (j = 0; j < nrows; j++) {
+ for (k = 0; k < ncols; k++) {
+ /* coordinates for points will be centers of 3D raster cells ! */
+ /* only produce points at non-NULL locations ! */
+ /* we have to store them as int values for performance reasons, so we will multiply by PRECISION */
+ /* PRECISION = 1000 = 10^3 = 3 decimal places precision! */
+ w = ascii[i][(nrows - 1) - j][k];
+ if (w == label) {
+ px[curPoint] =
+ ((long int)window.west + (window.ew_res * k)
+ + (window.ew_res * 0.5)) * (PRECISION);
+ py[curPoint] =
+ ((long int)window.south + (window.ns_res * j)
+ + (window.ns_res * 0.5)) * (PRECISION);
+ pz[curPoint] =
+ ((long int)window.top - (window.tb_res * i)
+ - (window.tb_res * 0.5)) * (PRECISION);
+ curPoint++;
+ }
+ }
+ }
+ }
+
+ if (VERBOSE) {
+ fprintf(stdout, "\nExporting 3D convex hull for layer %i: \n", l);
+ }
+ /* make 3D hull */
+ error = make3DHull(px, py, pz, num_points, tmpfile);
+
+ if (error < 0) {
+
+ fprintf(stdout,
+ "DEBUG: Simple planar hulls not implemented yet.\n");
+ fprintf(stdout, "DEBUG: This layer will not be exported.\n");
+
+ /* all points in the same voxel plane. Let's output a very simple hull with 8 vertices */
+
+ /* ONCE THIS IS IMPLEMENTED: call v.in.ascii for these hulls, too! */
+
+ if (l == 0) {
+ /* this is the bottom DEM: we will extend it half a slice (tbres) down */
+ fprintf(tmpfile,
+ "DEBUG: Simple planar hulls not implemented yet.\n");
+ fprintf(tmpfile, "DEBUG: This layer was not exported.\n");
+ }
+ else {
+ if (l == (nfiles - 1)) {
+ /* this is the top DEM: we will extend it half a slice (tbres) up */
+ fprintf(tmpfile,
+ "DEBUG: Simple planar hulls not implemented yet.\n");
+ fprintf(stdout, "DEBUG: This layer was not exported.\n");
+ }
+ else {
+ /* this DEM is somewhare in between, its vector representation will be planar! */
+ fprintf(tmpfile,
+ "DEBUG: Simple planar hulls not implemented yet.\n");
+ fprintf(stdout, "DEBUG: This layer was not exported.\n");
+ }
+ }
+ }
+
+ /* free memory for points */
+ free(px);
+ free(py);
+ free(pz);
+
+ fclose(tmpfile);
+
+ /* call v.in.ascii to create a native GRASS vector map */
+ /* DEBUG: remove once planar hulls are implemented */
+ if (error > -1) {
+ if (VERBOSE) {
+ fprintf(stdout,
+ "\nRunning v.in.ascii to create 3D points map: \n");
+ }
+ sprintf(tmp,
+ "v.in.ascii -z in=%s.hull.%i.txt out=%s_hull_%i form=standard",
+ parm.output->answer, l, parm.output->answer, l);
+
+ error = system(tmp);
+ if (error != EXIT_SUCCESS) {
+ G_fatal_error("Error running command: %s", tmp);
+ }
+
+ if (VERBOSE)
+ fprintf(stdout, " DONE \n");
+ }
+ }
+ /* MAIN LOOP END */
+}
+
+
+
+void output_vtk()
+{
+
+ char sys[2048];
+ int error;
+
+ if (flag.vtk->answer) {
+
+ /* write a VTK file for visualisation in paraview etc. */
+ if (VERBOSE)
+ fprintf(stdout, "\nCreating VTK file (%s.vtk): \n",
+ parm.output->answer);
+
+ if (flag.points->answer) {
+ sprintf(sys, "r3.out.vtk -p in=%s out=%s.vtk null=%f",
+ parm.output->answer, parm.output->answer, DNULLVALUE);
+
+ }
+ else {
+ sprintf(sys, "r3.out.vtk in=%s out=%s.vtk null=%f",
+ parm.output->answer, parm.output->answer, DNULLVALUE);
+ }
+ if (DEBUG > 0) {
+ fprintf(stdout, "%s\n", sys);
+ }
+ error = system(sys);
+ if (error != 0) {
+ G_fatal_error("Error running command: %s", sys);
+ }
+
+ }
+}
+
+
+/* ================================================================================
+
+ MAIN
+
+ =================================================================================== */
+
+
+int main(int argc, char *argv[])
+{
+ struct GModule *module;
+
+ struct Cell_head window;
+
+ int i, j, k;
+ int nrows, ncols;
+
+ int *fd;
+ int nfiles, nvalues, ncolors;
+ double dslices;
+ char *name, *mapset;
+ DCELL **dcell;
+
+ double ***ascii; /* 3D array that holds values to print into ASCII file */
+
+ unsigned short int *R;
+ unsigned short int *G;
+ unsigned short int *B;
+
+ char sys[2048];
+
+ module = G_define_module();
+ module->description =
+ "Creates a 3D raster model (voxels) from a series of raster DEMs";
+ /* do not pause after a warning message was displayed */
+ G_sleep_on_error(0);
+
+ /* DEFINE OPTIONS AND FLAGS */
+ /* input raster map */
+ parm.maps = G_define_standard_option(G_OPT_R_MAPS);
+ parm.maps->key = "map";
+ parm.maps->type = TYPE_STRING;
+ parm.maps->required = YES;
+ parm.maps->description =
+ "Input DEMs (at least 2) in raster format. Bottom DEM first";
+ parm.maps->gisprompt = "old,cell,raster";
+ parm.maps->multiple = YES;
+
+ /* Output map name */
+ parm.output = G_define_option();
+ parm.output->key = "output";
+ parm.output->type = TYPE_STRING;
+ parm.output->required = YES;
+ parm.output->description = "Output rast3d map name";
+
+ /* optional: specify a value for each layer */
+ parm.values = G_define_option();
+ parm.values->key = "labels";
+ parm.values->type = TYPE_DOUBLE;
+ parm.values->required = NO;
+ parm.values->description = "List of label values, one for each 3D layer";
+
+ /* optional: specify an RRR:GGG:BBB color triplet for each layer */
+ parm.colors = G_define_option();
+ parm.colors->key = "colors";
+ parm.colors->type = TYPE_STRING;
+ parm.colors->required = NO;
+ parm.colors->description =
+ "List of RRR:GGG:BBB color triplets, one for each layer";
+
+ /* optional: Raster map to represent topology errors in input DEMs */
+ parm.error = G_define_option();
+ parm.error->key = "errormap";
+ parm.error->type = TYPE_STRING;
+ parm.error->required = NO;
+ parm.error->description =
+ "Raster map to represent topology errors in input DEMs";
+
+ /* algorithm to use for flood filling */
+ parm.algorithm = G_define_option();
+ parm.algorithm->key = "algorithm";
+ parm.algorithm->type = TYPE_STRING;
+ parm.algorithm->required = NO;
+ parm.algorithm->options = "up,down";
+ parm.algorithm->answer = "up";
+ parm.algorithm->description = "3D flood fill algorithm to use";
+
+ /* 3D cell stats? */
+ flag.countcells = G_define_flag();
+ flag.countcells->key = 'c';
+ flag.countcells->description = "Calculate 3D cell counts for each layer.";
+
+ /* VTK point data output ? */
+ flag.fillnull = G_define_flag();
+ flag.fillnull->key = 'f';
+ flag.fillnull->description = "Fill through NULL value areas in DEMs";
+
+ /* export a 3D hull polygon for every layer ? */
+ flag.hull = G_define_flag();
+ flag.hull->key = 'h';
+ flag.hull->description = "Export convex hull polygons for layers";
+
+ /* GRASS ASCII vector point data output ? */
+ flag.grasspts = G_define_flag();
+ flag.grasspts->key = 'g';
+ flag.grasspts->description = "Export voxel model as vector points";
+
+ /* VTK point data output ? */
+ flag.points = G_define_flag();
+ flag.points->key = 'p';
+ flag.points->description = "Export VTK point data instead of cell data";
+
+ /* quiet operation? */
+ flag.quiet = G_define_flag();
+ flag.quiet->key = 'q';
+ flag.quiet->description = "Disable on-screen progress display";
+
+ /* skip DEM topology checks? */
+ flag.skip = G_define_flag();
+ flag.skip->key = 's';
+ flag.skip->description = "Skip topology error checking";
+
+ /* generate a VTK file? */
+ flag.vtk = G_define_flag();
+ flag.vtk->key = 'v';
+ flag.vtk->description =
+ "Generate a .vtk file for visualisation with e.g. paraview";
+
+ /* adjust Z-region settings to fit DEM values? */
+ flag.zadjust = G_define_flag();
+ flag.zadjust->key = 'z';
+ flag.zadjust->description = "Fit active region's Z range to input DEMs";
+
+ int error;
+
+
+ /* INIT GLOBAL VARIABLES */
+ VERBOSE = 1;
+ NSLICES = 0;
+ PRECISION = 1000;
+
+ /* setup some basic GIS stuff */
+ G_gisinit(argv[0]);
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ if (flag.quiet->answer) {
+ VERBOSE = 0;
+ }
+
+ Rast_get_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* */
+ /* STEP 1: VALIDATE input data and topology */
+ /* */
+
+ fprintf(stdout, "Validating data: \n");
+
+ /* check if input DEMs are available */
+ /* allocate all necessary handles */
+ for (nfiles = 0; parm.maps->answers[nfiles]; nfiles++) ;
+
+ if (DEBUG)
+ fprintf(stdout, "nfiles = %i\n", nfiles);
+
+ fd = (int *)G_malloc(nfiles * sizeof(int));
+ dcell = (DCELL **) G_malloc(nfiles * sizeof(DCELL *));
+
+ if (nfiles < 2) {
+ G_fatal_error("Need at least two raster maps!");
+ }
+
+ for (i = 0; i < nfiles; i++) {
+ dcell[i] = Rast_allocate_d_buf();
+ name = parm.maps->answers[i];
+ mapset = G_find_raster2(name, "");
+ if (!mapset)
+ G_fatal_error(_("%s - raster map not found"), name);
+ fd[i] = Rast_open_old(name, mapset);
+ if (fd[i] < 0)
+ G_fatal_error(_("%s - can't open raster map"), name);
+ }
+
+ nvalues = 0;
+ ncolors = 0;
+ /* if values and/or colors have been specified for the model: check these parameters */
+ if (parm.values->answers != NULL)
+ for (nvalues = 0; parm.values->answers[nvalues]; nvalues++) ;
+ if (nvalues > 0) {
+ if (nvalues != nfiles) {
+ G_fatal_error(_("Number of values (%i) does not match number of DEMs (%i)!"),
+ nvalues, nfiles);
+ }
+ }
+ if (parm.colors->answers != NULL) {
+ for (ncolors = 0; parm.colors->answers[ncolors]; ncolors++) ;
+ /* KILL ME */
+ fprintf(stdout, "DEBUG: custom RGB colouring not implemented yet.\n");
+ /* KILL ME */
+ }
+ if (ncolors > 0) {
+ R = G_malloc(ncolors * sizeof(unsigned short int));
+ G = G_malloc(ncolors * sizeof(unsigned short int));
+ B = G_malloc(ncolors * sizeof(unsigned short int));
+ if (ncolors != nfiles) {
+ G_fatal_error(_("Number of RGB triplets (%i) does not match number of DEMs (%i)!"),
+ ncolors, nfiles);
+ }
+ for (i = 0; i < ncolors; i++) {
+ if (!isRGB(parm.colors->answers[i], &R[i], &G[i], &B[i])) {
+ G_fatal_error(_("%s is not a valid RRR:GGG:BBB code!"),
+ parm.colors->answers[i]);
+ }
+ }
+ }
+
+ check_zrange(fd, dcell, nfiles);
+
+ if (!flag.skip->answer) {
+ check_topologies(fd, dcell, nfiles);
+ }
+
+ /* let's see how many slices we will make */
+ dslices = ((window.top - window.bottom) / window.tb_res);
+ NSLICES = (int)rint(dslices);
+ if (NSLICES < 1) {
+ G_fatal_error
+ ("This would produce < 1 slice. Adjust active 3d region resolution.");
+ }
+ if (DEBUG) {
+ fprintf(stdout, "%i DEMS, NSLICES = %i (%f)\n", nfiles, NSLICES,
+ dslices);
+ fflush(stdout);
+ }
+
+
+ /* */
+ /* STEP 2: INTERPOLATE volume from DEMs */
+ /* */
+
+ /* allocate array for storing voxel data and initialise it with DNULLVALUE */
+ ascii = (double ***)G_malloc(NSLICES * sizeof(double **));
+ for (i = 0; i < NSLICES; i++) {
+ ascii[i] = (double **)G_malloc(nrows * sizeof(double *));
+ for (j = 0; j < nrows; j++) {
+ ascii[i][j] = G_malloc(ncols * sizeof(double));
+ for (k = 0; k < ncols; k++) {
+ ascii[i][j][k] = DNULLVALUE;
+ }
+ }
+ }
+ interpolate(fd, dcell, nfiles, nvalues, ascii);
+
+
+ /* */
+ /* STEP 3: OUTPUT volume data */
+ /* */
+
+ /* TEMP SECTION */
+ /* this will only stay in here temporarily until this
+ program will be able to write voxel maps directly !!! */
+ output_ascii(ascii);
+ /* now call r3.in.ascii to write a voxel model for us */
+
+ if (VERBOSE)
+ fprintf(stdout, "Running r3.in.ascii to create voxel model: \n");
+
+ sprintf(sys, "r3.in.ascii in=%s out=%s nv=%f", TMPFILE,
+ parm.output->answer, DNULLVALUE);
+ if (DEBUG > 0) {
+ fprintf(stdout, "%s\n", sys);
+ }
+ error = system(sys);
+ if (error != EXIT_SUCCESS) {
+ G_fatal_error("Error running command: %s", sys);
+ }
+ if (VERBOSE)
+ fprintf(stdout, " DONE \n");
+
+ /* END (TEMP SECTION) */
+
+
+ /* VTK output */
+ if (flag.vtk->answer) {
+ output_vtk();
+ }
+
+ /* ASCII vector points output */
+ if (flag.grasspts->answer) {
+ output_ascii_points(ascii);
+ }
+
+ /* ASCII convex hulls ouput */
+ if (flag.hull->answer) {
+ output_ascii_hulls(ascii, nvalues, nfiles);
+ }
+
+ /* show cell count stats? */
+ if (flag.countcells->answer) {
+ count_3d_cells(ascii, nfiles, nvalues);
+ }
+
+ fprintf(stdout, "\nJOB DONE.\n");
+ fflush(stdout);
+
+ return (EXIT_SUCCESS);
+}
Property changes on: grass-addons/grass7/raster/r.vol.dem/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list