GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
volume.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 
6 #include <grass/raster.h>
7 #include <grass/raster3d.h>
8 
9 /*---------------------------------------------------------------------------*/
10 
11 static int verifyVolumeVertices(void *map, double v[2][2][2][3])
12 {
13  if (!(Rast3d_is_valid_location(map, v[0][0][0][0], v[0][0][0][1],
14  v[0][0][0][2]) &&
15  Rast3d_is_valid_location(map, v[0][0][1][0], v[0][0][1][1],
16  v[0][0][1][2]) &&
17  Rast3d_is_valid_location(map, v[0][1][0][0], v[0][1][0][1],
18  v[0][1][0][2]) &&
19  Rast3d_is_valid_location(map, v[0][1][1][0], v[0][1][1][1],
20  v[0][1][1][2]) &&
21  Rast3d_is_valid_location(map, v[1][0][0][0], v[1][0][0][1],
22  v[1][0][0][2]) &&
23  Rast3d_is_valid_location(map, v[1][0][1][0], v[1][0][1][1],
24  v[1][0][1][2]) &&
25  Rast3d_is_valid_location(map, v[1][1][0][0], v[1][1][0][1],
26  v[1][1][0][2]) &&
27  Rast3d_is_valid_location(map, v[1][1][1][0], v[1][1][1][1],
28  v[1][1][1][2])))
29  Rast3d_fatal_error("verifyCubeVertices: volume vertex out of range");
30  return 0;
31 }
32 
33 
34 /*---------------------------------------------------------------------------*/
35 
36 static int verifyVolumeEdges(int nx, int ny, int nz)
37 {
38  if ((nx <= 0) || (ny <= 0) || (nz <= 0))
39  Rast3d_fatal_error("verifyCubeEdges: Volume edge out of range");
40  return 0;
41 }
42 
43 /*---------------------------------------------------------------------------*/
44 
45 void
46 Rast3d_get_volume_a(void *map, double u[2][2][2][3], int nx, int ny, int nz,
47  void *volumeBuf, int type)
48 {
49  typedef double doubleArray[3];
50 
51  doubleArray *u000, *u001, *u010, *u011;
52  doubleArray *u100, *u101, *u110, *u111;
53  double v00[3], v01[3], v10[3], v11[3];
54  double v0[3], v1[3];
55  double v[3];
56  double r, rp, s, sp, t, tp;
57  double dx, dy, dz;
58  int x, y, z, nxp, nyp, nzp;
59  double *doubleBuf;
60  float *floatBuf;
61 
62  doubleBuf = (double *)volumeBuf;
63  floatBuf = (float *)volumeBuf;
64 
65  verifyVolumeVertices(map, u);
66  verifyVolumeEdges(nx, ny, nz);
67 
68  nxp = nx * 2 + 1;
69  nyp = ny * 2 + 1;
70  nzp = nz * 2 + 1;
71 
72  u000 = (doubleArray *) u[0][0][0];
73  u001 = (doubleArray *) u[0][0][1];
74  u010 = (doubleArray *) u[0][1][0];
75  u011 = (doubleArray *) u[0][1][1];
76  u100 = (doubleArray *) u[1][0][0];
77  u101 = (doubleArray *) u[1][0][1];
78  u110 = (doubleArray *) u[1][1][0];
79  u111 = (doubleArray *) u[1][1][1];
80 
81  for (dz = 1; dz < nzp; dz += 2) {
82  r = 1. - (rp = dz / nz / 2.);
83  v00[0] = r * (*u000)[0] + rp * (*u100)[0];
84  v00[1] = r * (*u000)[1] + rp * (*u100)[1];
85  v00[2] = r * (*u000)[2] + rp * (*u100)[2];
86 
87  v01[0] = r * (*u001)[0] + rp * (*u101)[0];
88  v01[1] = r * (*u001)[1] + rp * (*u101)[1];
89  v01[2] = r * (*u001)[2] + rp * (*u101)[2];
90 
91  v10[0] = r * (*u010)[0] + rp * (*u110)[0];
92  v10[1] = r * (*u010)[1] + rp * (*u110)[1];
93  v10[2] = r * (*u010)[2] + rp * (*u110)[2];
94 
95  v11[0] = r * (*u011)[0] + rp * (*u111)[0];
96  v11[1] = r * (*u011)[1] + rp * (*u111)[1];
97  v11[2] = r * (*u011)[2] + rp * (*u111)[2];
98 
99  for (dy = 1; dy < nyp; dy += 2) {
100  s = 1. - (sp = dy / ny / 2.);
101  v0[0] = s * v00[0] + sp * v10[0];
102  v0[1] = s * v00[1] + sp * v10[1];
103  v0[2] = s * v00[2] + sp * v10[2];
104 
105  v1[0] = s * v01[0] + sp * v11[0];
106  v1[1] = s * v01[1] + sp * v11[1];
107  v1[2] = s * v01[2] + sp * v11[2];
108 
109  for (dx = 1; dx < nxp; dx += 2) {
110  t = 1. - (tp = dx / nx / 2.);
111  v[0] = t * v0[0] + tp * v1[0];
112  v[1] = t * v0[1] + tp * v1[1];
113  v[2] = t * v0[2] + tp * v1[2];
114 
115  Rast3d_location2coord2(map, v[0], v[1], v[2], &x, &y, &z);
116  /* DEBUG
117  printf ("(%d %d %d) (%lf %lf %lf) (%d %d %d) %lf\n",
118  (int) dx / 2, (int) dy / 2, (int) dz / 2,
119  v[0], v[1], v[2],
120  x, y, z,
121  Rast3d_get_double_region (map, x, y, z));
122  */
123  if (type == DCELL_TYPE)
124  *(doubleBuf + ((int)dz / 2) * nx * ny +
125  ((int)dy / 2) * nx + (int)dx / 2) =
126 Rast3d_get_double_region(map, x, y, z);
127  else
128  *(floatBuf + ((int)dz / 2) * nx * ny +
129  ((int)dy / 2) * nx + (int)dx / 2) =
130 Rast3d_get_float_region(map, x, y, z);
131  }
132  }
133  }
134 }
135 
136 /*---------------------------------------------------------------------------*/
137 
138 void
140  double originNorth, double originWest, double originBottom,
141  double vxNorth, double vxWest, double vxBottom,
142  double vyNorth, double vyWest, double vyBottom,
143  double vzNorth, double vzWest, double vzBottom,
144  int nx, int ny, int nz, void *volumeBuf, int type)
145 {
146  double u[2][2][2][3];
147 
148  u[0][0][0][0] = originNorth;
149  u[0][0][0][1] = originWest;
150  u[0][0][0][2] = originBottom;
151 
152  u[0][0][1][0] = vxNorth;
153  u[0][0][1][1] = vxWest;
154  u[0][0][1][2] = vxBottom;
155 
156  u[1][0][0][0] = vzNorth;
157  u[1][0][0][1] = vzWest;
158  u[1][0][0][2] = vzBottom;
159 
160  u[1][0][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[1][0][0][0];
161  u[1][0][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[1][0][0][1];
162  u[1][0][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[1][0][0][2];
163 
164  u[0][1][0][0] = vyNorth;
165  u[0][1][0][1] = vyWest;
166  u[0][1][0][2] = vyBottom;
167 
168  u[0][1][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[0][1][0][0];
169  u[0][1][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[0][1][0][1];
170  u[0][1][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[0][1][0][2];
171 
172  u[1][1][0][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][0][0];
173  u[1][1][0][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][0][1];
174  u[1][1][0][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][0][2];
175 
176  u[1][1][1][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][1][0];
177  u[1][1][1][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][1][1];
178  u[1][1][1][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][1][2];
179 
180  Rast3d_get_volume_a(map, u, nx, ny, nz, volumeBuf, type);
181 }
182 
183 /*---------------------------------------------------------------------------*/
184 
185 void
187  double originNorth, double originWest,
188  double originBottom, double lengthNorth,
189  double lengthWest, double lengthBottom, int nx, int ny,
190  int nz, void *volumeBuf, int type)
191 {
192  Rast3d_get_volume(map,
193  originNorth, originWest, originBottom,
194  originNorth + lengthNorth, originWest, originBottom,
195  originNorth, originWest + lengthWest, originBottom,
196  originNorth, originWest, originBottom + lengthBottom,
197  nx, ny, nz, volumeBuf, type);
198 }
199 
200 /*---------------------------------------------------------------------------*/
201 
202 void
203 Rast3d_make_aligned_volume_file(void *map, const char *fileName,
204  double originNorth, double originWest,
205  double originBottom, double lengthNorth,
206  double lengthWest, double lengthBottom, int nx,
207  int ny, int nz)
208 {
209  void *volumeBuf;
210  void *mapVolume;
211  int x, y, z, eltLength;
212  RASTER3D_Region region;
213 
214  volumeBuf = Rast3d_malloc(nx * ny * nz * sizeof(Rast3d_get_file_type()));
215  if (volumeBuf == NULL)
216  Rast3d_fatal_error("Rast3d_make_aligned_volume_file: error in Rast3d_malloc");
217 
219  originNorth, originWest, originBottom,
220  lengthNorth, lengthWest, lengthBottom,
221  nx, ny, nz, volumeBuf, Rast3d_get_file_type());
222 
223  region.north = originNorth;
224  region.south = originNorth + lengthNorth;
225  region.east = originWest;
226  region.west = originWest + lengthWest;
227  region.top = originBottom;
228  region.bottom = originBottom + lengthBottom;
229 
230  region.rows = ny;
231  region.cols = nx;
232  region.depths = nz;
233 
234  mapVolume = Rast3d_open_cell_new(fileName, Rast3d_get_file_type(),
235  RASTER3D_USE_CACHE_DEFAULT, &region);
236  if (mapVolume == NULL)
237  Rast3d_fatal_error("Rast3d_make_aligned_volume_file: error in Rast3d_open_cell_new");
238 
239  eltLength = Rast3d_length(Rast3d_get_file_type());
240 
241  for (z = 0; z < nz; z++) {
242  for (y = 0; y < ny; y++) {
243  for (x = 0; x < nx; x++) {
244  /* Rast3d_putValueRegion? */
245  if (!Rast3d_put_value(mapVolume, x, y, z,
246  G_incr_void_ptr(volumeBuf,
247  (z * ny * nx + y * nx +
248  x) * eltLength),
249  Rast3d_file_type_map(mapVolume)))
251  ("Rast3d_make_aligned_volume_file: error in Rast3d_put_value");
252  }
253  }
254  }
255 
256  if (!Rast3d_close(mapVolume))
257  Rast3d_fatal_error("Rast3d_make_aligned_volume_file: error in Rast3d_close");
258 
259  Rast3d_free(volumeBuf);
260 }
#define RASTER3D_USE_CACHE_DEFAULT
Definition: raster3d.h:20
double east
Definition: raster3d.h:51
int Rast3d_put_value(RASTER3D_Map *, int, int, int, const void *, int)
After converting *value of type into the type specified at the initialization time (i...
Definition: putvalue.c:96
double top
Definition: raster3d.h:52
int Rast3d_file_type_map(RASTER3D_Map *)
Returns the type with which tiles of map are stored on file.
Definition: headerinfo.c:281
double bottom
Definition: raster3d.h:52
void * Rast3d_malloc(int)
Same as malloc (nBytes), except that in case of error Rast3d_error() is invoked.
int Rast3d_is_valid_location(RASTER3D_Region *, double, double, double)
Returns 1 if region-coordinates (north, east, top) are inside the region of map. Returns 0 otherwise...
Definition: region.c:263
int Rast3d_close(RASTER3D_Map *)
Close 3D raster map files.
#define G_incr_void_ptr(ptr, size)
Definition: defs/gis.h:100
void * Rast3d_open_cell_new(const char *, int, int, RASTER3D_Region *)
Opens new g3d-file with name in the current mapset. Tiles are stored in memory with type which must b...
void Rast3d_get_aligned_volume(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: volume.c:186
#define NULL
Definition: ccmath.h:32
#define x
float Rast3d_get_float_region(RASTER3D_Map *, int, int, int)
Is equivalent to Rast3d_get_value_region (map, x, y, z, &value, FCELL_TYPE); return value...
Definition: getvalue.c:165
void Rast3d_fatal_error(const char *,...) __attribute__((format(printf
double north
Definition: raster3d.h:50
void Rast3d_get_volume(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: volume.c:139
double t
Definition: r_raster.c:39
#define DCELL_TYPE
Definition: raster.h:13
double west
Definition: raster3d.h:51
void Rast3d_location2coord2(RASTER3D_Region *, double, double, double, int *, int *, int *)
Converts region-coordinates (north, east, top) into cell-coordinates (x, y, z). This function calls R...
Definition: region.c:347
double Rast3d_get_double_region(RASTER3D_Map *, int, int, int)
Is equivalent to Rast3d_get_value_region (map, x, y, z, &value, DCELL_TYPE); return value...
Definition: getvalue.c:207
int Rast3d_get_file_type(void)
get G3d file type
Definition: defaults.c:227
void Rast3d_get_volume_a(void *map, double u[2][2][2][3], int nx, int ny, int nz, void *volumeBuf, int type)
Definition: volume.c:46
double south
Definition: raster3d.h:50
void Rast3d_make_aligned_volume_file(void *map, const char *fileName, double originNorth, double originWest, double originBottom, double lengthNorth, double lengthWest, double lengthBottom, int nx, int ny, int nz)
Definition: volume.c:203
double r
Definition: r_raster.c:39
int Rast3d_length(int)
Definition: raster3d/misc.c:78
void Rast3d_free(void *)
Same as free (ptr).