GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
g3dvolume.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <sys/types.h>
4 #include <unistd.h>
5 #include <grass/G3d.h>
6 
7 /*---------------------------------------------------------------------------*/
8 
9 static int verifyVolumeVertices(map, v)
10 
11  void *map;
12  double v[2][2][2][3];
13 
14 {
15  if (!(G3d_isValidLocation(map, v[0][0][0][0], v[0][0][0][1],
16  v[0][0][0][2]) &&
17  G3d_isValidLocation(map, v[0][0][1][0], v[0][0][1][1],
18  v[0][0][1][2]) &&
19  G3d_isValidLocation(map, v[0][1][0][0], v[0][1][0][1],
20  v[0][1][0][2]) &&
21  G3d_isValidLocation(map, v[0][1][1][0], v[0][1][1][1],
22  v[0][1][1][2]) &&
23  G3d_isValidLocation(map, v[1][0][0][0], v[1][0][0][1],
24  v[1][0][0][2]) &&
25  G3d_isValidLocation(map, v[1][0][1][0], v[1][0][1][1],
26  v[1][0][1][2]) &&
27  G3d_isValidLocation(map, v[1][1][0][0], v[1][1][0][1],
28  v[1][1][0][2]) &&
29  G3d_isValidLocation(map, v[1][1][1][0], v[1][1][1][1],
30  v[1][1][1][2])))
31  G3d_fatalError("verifyCubeVertices: volume vertex out of range");
32  return 0;
33 }
34 
35 
36 /*---------------------------------------------------------------------------*/
37 
38 static int verifyVolumeEdges(nx, ny, nz)
39 
40  int nx, ny, nz;
41 
42 {
43  if ((nx <= 0) || (ny <= 0) || (nz <= 0))
44  G3d_fatalError("verifyCubeEdges: Volume edge out of range");
45  return 0;
46 }
47 
48 /*---------------------------------------------------------------------------*/
49 
50 void
51 G3d_getVolumeA(void *map, double u[2][2][2][3], int nx, int ny, int nz,
52  void *volumeBuf, int type)
53 {
54  typedef double doubleArray[3];
55 
56  doubleArray *u000, *u001, *u010, *u011;
57  doubleArray *u100, *u101, *u110, *u111;
58  double v00[3], v01[3], v10[3], v11[3];
59  double v0[3], v1[3];
60  double v[3];
61  double r, rp, s, sp, t, tp;
62  double dx, dy, dz;
63  int x, y, z, nxp, nyp, nzp;
64  double *doubleBuf;
65  float *floatBuf;
66 
67  doubleBuf = (double *)volumeBuf;
68  floatBuf = (float *)volumeBuf;
69 
70  verifyVolumeVertices(map, u);
71  verifyVolumeEdges(nx, ny, nz);
72 
73  nxp = nx * 2 + 1;
74  nyp = ny * 2 + 1;
75  nzp = nz * 2 + 1;
76 
77  u000 = (doubleArray *) u[0][0][0];
78  u001 = (doubleArray *) u[0][0][1];
79  u010 = (doubleArray *) u[0][1][0];
80  u011 = (doubleArray *) u[0][1][1];
81  u100 = (doubleArray *) u[1][0][0];
82  u101 = (doubleArray *) u[1][0][1];
83  u110 = (doubleArray *) u[1][1][0];
84  u111 = (doubleArray *) u[1][1][1];
85 
86  for (dz = 1; dz < nzp; dz += 2) {
87  r = 1. - (rp = dz / nz / 2.);
88  v00[0] = r * (*u000)[0] + rp * (*u100)[0];
89  v00[1] = r * (*u000)[1] + rp * (*u100)[1];
90  v00[2] = r * (*u000)[2] + rp * (*u100)[2];
91 
92  v01[0] = r * (*u001)[0] + rp * (*u101)[0];
93  v01[1] = r * (*u001)[1] + rp * (*u101)[1];
94  v01[2] = r * (*u001)[2] + rp * (*u101)[2];
95 
96  v10[0] = r * (*u010)[0] + rp * (*u110)[0];
97  v10[1] = r * (*u010)[1] + rp * (*u110)[1];
98  v10[2] = r * (*u010)[2] + rp * (*u110)[2];
99 
100  v11[0] = r * (*u011)[0] + rp * (*u111)[0];
101  v11[1] = r * (*u011)[1] + rp * (*u111)[1];
102  v11[2] = r * (*u011)[2] + rp * (*u111)[2];
103 
104  for (dy = 1; dy < nyp; dy += 2) {
105  s = 1. - (sp = dy / ny / 2.);
106  v0[0] = s * v00[0] + sp * v10[0];
107  v0[1] = s * v00[1] + sp * v10[1];
108  v0[2] = s * v00[2] + sp * v10[2];
109 
110  v1[0] = s * v01[0] + sp * v11[0];
111  v1[1] = s * v01[1] + sp * v11[1];
112  v1[2] = s * v01[2] + sp * v11[2];
113 
114  for (dx = 1; dx < nxp; dx += 2) {
115  t = 1. - (tp = dx / nx / 2.);
116  v[0] = t * v0[0] + tp * v1[0];
117  v[1] = t * v0[1] + tp * v1[1];
118  v[2] = t * v0[2] + tp * v1[2];
119 
120  G3d_location2coord(map, v[0], v[1], v[2], &x, &y, &z);
121  /* DEBUG
122  printf ("(%d %d %d) (%lf %lf %lf) (%d %d %d) %lf\n",
123  (int) dx / 2, (int) dy / 2, (int) dz / 2,
124  v[0], v[1], v[2],
125  x, y, z,
126  G3d_getDoubleRegion (map, x, y, z));
127  */
128  if (type == DCELL_TYPE)
129  *(doubleBuf + ((int)dz / 2) * nx * ny +
130  ((int)dy / 2) * nx + (int)dx / 2) =
131 G3d_getDoubleRegion(map, x, y, z);
132  else
133  *(floatBuf + ((int)dz / 2) * nx * ny +
134  ((int)dy / 2) * nx + (int)dx / 2) =
135 G3d_getFloatRegion(map, x, y, z);
136  }
137  }
138  }
139 }
140 
141 /*---------------------------------------------------------------------------*/
142 
143 void
144 G3d_getVolume(void *map,
145  double originNorth, double originWest, double originBottom,
146  double vxNorth, double vxWest, double vxBottom,
147  double vyNorth, double vyWest, double vyBottom,
148  double vzNorth, double vzWest, double vzBottom,
149  int nx, int ny, int nz, void *volumeBuf, int type)
150 {
151  double u[2][2][2][3];
152 
153  u[0][0][0][0] = originNorth;
154  u[0][0][0][1] = originWest;
155  u[0][0][0][2] = originBottom;
156 
157  u[0][0][1][0] = vxNorth;
158  u[0][0][1][1] = vxWest;
159  u[0][0][1][2] = vxBottom;
160 
161  u[1][0][0][0] = vzNorth;
162  u[1][0][0][1] = vzWest;
163  u[1][0][0][2] = vzBottom;
164 
165  u[1][0][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[1][0][0][0];
166  u[1][0][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[1][0][0][1];
167  u[1][0][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[1][0][0][2];
168 
169  u[0][1][0][0] = vyNorth;
170  u[0][1][0][1] = vyWest;
171  u[0][1][0][2] = vyBottom;
172 
173  u[0][1][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[0][1][0][0];
174  u[0][1][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[0][1][0][1];
175  u[0][1][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[0][1][0][2];
176 
177  u[1][1][0][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][0][0];
178  u[1][1][0][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][0][1];
179  u[1][1][0][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][0][2];
180 
181  u[1][1][1][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][1][0];
182  u[1][1][1][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][1][1];
183  u[1][1][1][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][1][2];
184 
185  G3d_getVolumeA(map, u, nx, ny, nz, volumeBuf, type);
186 }
187 
188 /*---------------------------------------------------------------------------*/
189 
190 void
191 G3d_getAlignedVolume(void *map,
192  double originNorth, double originWest,
193  double originBottom, double lengthNorth,
194  double lengthWest, double lengthBottom, int nx, int ny,
195  int nz, void *volumeBuf, int type)
196 {
197  G3d_getVolume(map,
198  originNorth, originWest, originBottom,
199  originNorth + lengthNorth, originWest, originBottom,
200  originNorth, originWest + lengthWest, originBottom,
201  originNorth, originWest, originBottom + lengthBottom,
202  nx, ny, nz, volumeBuf, type);
203 }
204 
205 /*---------------------------------------------------------------------------*/
206 
207 void
208 G3d_makeAlignedVolumeFile(void *map, const char *fileName,
209  double originNorth, double originWest,
210  double originBottom, double lengthNorth,
211  double lengthWest, double lengthBottom, int nx,
212  int ny, int nz)
213 {
214  void *volumeBuf;
215  void *mapVolume;
216  int x, y, z, eltLength;
217  G3D_Region region;
218 
219  volumeBuf = G3d_malloc(nx * ny * nz * sizeof(G3d_getFileType()));
220  if (volumeBuf == NULL)
221  G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_malloc");
222 
223  G3d_getAlignedVolume(map,
224  originNorth, originWest, originBottom,
225  lengthNorth, lengthWest, lengthBottom,
226  nx, ny, nz, volumeBuf, G3d_getFileType());
227 
228  region.north = originNorth;
229  region.south = originNorth + lengthNorth;
230  region.east = originWest;
231  region.west = originWest + lengthWest;
232  region.top = originBottom;
233  region.bottom = originBottom + lengthBottom;
234 
235  region.rows = ny;
236  region.cols = nx;
237  region.depths = nz;
238 
239  mapVolume = G3d_openCellNew(fileName, G3d_getFileType(),
240  G3D_USE_CACHE_DEFAULT, &region);
241  if (mapVolume == NULL)
242  G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_openCellNew");
243 
244  eltLength = G3d_length(G3d_getFileType());
245 
246  for (z = 0; z < nz; z++) {
247  for (y = 0; y < ny; y++) {
248  for (x = 0; x < nx; x++) {
249  /* G3d_putValueRegion? */
250  if (!G3d_putValue(mapVolume, x, y, z,
251  G_incr_void_ptr(volumeBuf,
252  (z * ny * nx + y * nx +
253  x) * eltLength),
254  G3d_fileTypeMap(mapVolume)))
255  G3d_fatalError
256  ("G3d_makeAlignedVolumeFile: error in G3d_putValue");
257  }
258  }
259  }
260 
261  if (!G3d_closeCell(mapVolume))
262  G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_closeCell");
263 
264  G3d_free(volumeBuf);
265 }
G3d_isValidLocation
int G3d_isValidLocation(G3D_Map *map, double north, double east, double top)
Returns 1 if region-coordinates (north, west, bottom) are inside the region of map. Returns 0 otherwise.
Definition: tilemath.c:327
G3d_free
void G3d_free(void *buf)
Same as free (ptr).
Definition: g3dalloc.c:71
r
float r
Definition: named_colr.c:8
G3d_getVolume
void G3d_getVolume(void *map, double originNorth, double originWest, double originBottom, double vxNorth, double vxWest, double vxBottom, double vyNorth, double vyWest, double vyBottom, double vzNorth, double vzWest, double vzBottom, int nx, int ny, int nz, void *volumeBuf, int type)
Definition: g3dvolume.c:144
y
int y
Definition: plot.c:34
G_incr_void_ptr
void * G_incr_void_ptr(const void *ptr, const size_t size)
Advance void pointer.
Definition: gis/raster.c:33
G3d_getAlignedVolume
void G3d_getAlignedVolume(void *map, double originNorth, double originWest, double originBottom, double lengthNorth, double lengthWest, double lengthBottom, int nx, int ny, int nz, void *volumeBuf, int type)
Definition: g3dvolume.c:191
G3d_location2coord
void G3d_location2coord(G3D_Map *map, double north, double east, double top, int *x, int *y, int *z)
Converts region-coordinates (north, west, bottom) into cell-coordinates (x, y, z).
Definition: tilemath.c:355
int
int
Definition: g3dcolor.c:48
G3d_closeCell
int G3d_closeCell(G3D_Map *map)
Closes g3d-file. If map is new and cache-mode is used for map then every tile which is not flushed be...
Definition: g3dclose.c:144
G3d_makeAlignedVolumeFile
void G3d_makeAlignedVolumeFile(void *map, const char *fileName, double originNorth, double originWest, double originBottom, double lengthNorth, double lengthWest, double lengthBottom, int nx, int ny, int nz)
Definition: g3dvolume.c:208
G3d_getDoubleRegion
double G3d_getDoubleRegion()
Is equivalent to G3d_getValueRegion (map, x, y, z, &amp;value, DCELL_TYPE); return value.
NULL
return NULL
Definition: dbfopen.c:1394
G3d_malloc
void * G3d_malloc(int nBytes)
Same as malloc (nBytes), except that in case of error G3d_error() is invoked.
Definition: g3dalloc.c:24
G3d_length
int G3d_length(int t)
Definition: g3dmisc.c:77
G3d_fileTypeMap
int G3d_fileTypeMap(G3D_Map *map)
Returns the type with which tiles of map are stored on file.
Definition: headerinfo.c:178
G3d_getVolumeA
void G3d_getVolumeA(void *map, double u[2][2][2][3], int nx, int ny, int nz, void *volumeBuf, int type)
Definition: g3dvolume.c:51
G3d_fatalError
void G3d_fatalError(const char *,...)
This function prints the error message msg, and terminates the program with an error status...
Definition: g3derror.c:58
G3d_getFloatRegion
float G3d_getFloatRegion(G3D_Map *map, int x, int y, int z)
Is equivalent to G3d_getValueRegion (map, x, y, z, &amp;value, FCELL_TYPE); return value.
Definition: tileio.c:172
G3d_putValue
int G3d_putValue(G3D_Map *map, int x, int y, int z, const void *value, int type)
After converting *value of type into the type specified at the initialization time (i...
Definition: tilewrite.c:540
G3d_openCellNew
void * G3d_openCellNew(const char *name, int typeIntern, int cache, G3D_Region *region)
Opens new g3d-file with name in the current mapset. Tiles are stored in memory with type which must b...
Definition: g3dopen.c:213
G3d_getFileType
int G3d_getFileType()
get G3d file type
Definition: g3ddefaults.c:243