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