GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
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
11static 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
35static 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
44void 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
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
134void 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
180void 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{
189 ny, nz, volumeBuf, type);
190}
191
192/*---------------------------------------------------------------------------*/
193
194void 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
211 lengthNorth, lengthWest, lengthBottom, nx, ny, nz,
213
214 region.north = originNorth;
215 region.south = originNorth + lengthNorth;
216 region.east = originWest;
217 region.west = originWest + lengthWest;
218 region.top = originBottom;
220
221 region.rows = ny;
222 region.cols = nx;
223 region.depths = nz;
224
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,
240 (z * ny * nx + y * nx + x) * eltLength),
242 Rast3d_fatal_error("Rast3d_make_aligned_volume_file: error "
243 "in Rast3d_put_value");
244 }
245 }
246 }
247
250 "Rast3d_make_aligned_volume_file: error in Rast3d_close");
251
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:263
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:159
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:249
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...
int Rast3d_get_file_type(void)
get G3d file type
Definition defaults.c:212
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:89
int Rast3d_length(int)
void * Rast3d_malloc(int)
Same as malloc (nBytes), except that in case of error Rast3d_error() is invoked.
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:330
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:201
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 bottom
Definition raster3d.h:51
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