GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <sys/types.h> 00004 #include <unistd.h> 00005 #include <grass/G3d.h> 00006 00007 /*---------------------------------------------------------------------------*/ 00008 00009 static int verifyVolumeVertices(map, v) 00010 00011 void *map; 00012 double v[2][2][2][3]; 00013 00014 { 00015 if (!(G3d_isValidLocation(map, v[0][0][0][0], v[0][0][0][1], 00016 v[0][0][0][2]) && 00017 G3d_isValidLocation(map, v[0][0][1][0], v[0][0][1][1], 00018 v[0][0][1][2]) && 00019 G3d_isValidLocation(map, v[0][1][0][0], v[0][1][0][1], 00020 v[0][1][0][2]) && 00021 G3d_isValidLocation(map, v[0][1][1][0], v[0][1][1][1], 00022 v[0][1][1][2]) && 00023 G3d_isValidLocation(map, v[1][0][0][0], v[1][0][0][1], 00024 v[1][0][0][2]) && 00025 G3d_isValidLocation(map, v[1][0][1][0], v[1][0][1][1], 00026 v[1][0][1][2]) && 00027 G3d_isValidLocation(map, v[1][1][0][0], v[1][1][0][1], 00028 v[1][1][0][2]) && 00029 G3d_isValidLocation(map, v[1][1][1][0], v[1][1][1][1], 00030 v[1][1][1][2]))) 00031 G3d_fatalError("verifyCubeVertices: volume vertex out of range"); 00032 return 0; 00033 } 00034 00035 00036 /*---------------------------------------------------------------------------*/ 00037 00038 static int verifyVolumeEdges(nx, ny, nz) 00039 00040 int nx, ny, nz; 00041 00042 { 00043 if ((nx <= 0) || (ny <= 0) || (nz <= 0)) 00044 G3d_fatalError("verifyCubeEdges: Volume edge out of range"); 00045 return 0; 00046 } 00047 00048 /*---------------------------------------------------------------------------*/ 00049 00050 void 00051 G3d_getVolumeA(void *map, double u[2][2][2][3], int nx, int ny, int nz, 00052 void *volumeBuf, int type) 00053 { 00054 typedef double doubleArray[3]; 00055 00056 doubleArray *u000, *u001, *u010, *u011; 00057 doubleArray *u100, *u101, *u110, *u111; 00058 double v00[3], v01[3], v10[3], v11[3]; 00059 double v0[3], v1[3]; 00060 double v[3]; 00061 double r, rp, s, sp, t, tp; 00062 double dx, dy, dz; 00063 int x, y, z, nxp, nyp, nzp; 00064 double *doubleBuf; 00065 float *floatBuf; 00066 00067 doubleBuf = (double *)volumeBuf; 00068 floatBuf = (float *)volumeBuf; 00069 00070 verifyVolumeVertices(map, u); 00071 verifyVolumeEdges(nx, ny, nz); 00072 00073 nxp = nx * 2 + 1; 00074 nyp = ny * 2 + 1; 00075 nzp = nz * 2 + 1; 00076 00077 u000 = (doubleArray *) u[0][0][0]; 00078 u001 = (doubleArray *) u[0][0][1]; 00079 u010 = (doubleArray *) u[0][1][0]; 00080 u011 = (doubleArray *) u[0][1][1]; 00081 u100 = (doubleArray *) u[1][0][0]; 00082 u101 = (doubleArray *) u[1][0][1]; 00083 u110 = (doubleArray *) u[1][1][0]; 00084 u111 = (doubleArray *) u[1][1][1]; 00085 00086 for (dz = 1; dz < nzp; dz += 2) { 00087 r = 1. - (rp = dz / nz / 2.); 00088 v00[0] = r * (*u000)[0] + rp * (*u100)[0]; 00089 v00[1] = r * (*u000)[1] + rp * (*u100)[1]; 00090 v00[2] = r * (*u000)[2] + rp * (*u100)[2]; 00091 00092 v01[0] = r * (*u001)[0] + rp * (*u101)[0]; 00093 v01[1] = r * (*u001)[1] + rp * (*u101)[1]; 00094 v01[2] = r * (*u001)[2] + rp * (*u101)[2]; 00095 00096 v10[0] = r * (*u010)[0] + rp * (*u110)[0]; 00097 v10[1] = r * (*u010)[1] + rp * (*u110)[1]; 00098 v10[2] = r * (*u010)[2] + rp * (*u110)[2]; 00099 00100 v11[0] = r * (*u011)[0] + rp * (*u111)[0]; 00101 v11[1] = r * (*u011)[1] + rp * (*u111)[1]; 00102 v11[2] = r * (*u011)[2] + rp * (*u111)[2]; 00103 00104 for (dy = 1; dy < nyp; dy += 2) { 00105 s = 1. - (sp = dy / ny / 2.); 00106 v0[0] = s * v00[0] + sp * v10[0]; 00107 v0[1] = s * v00[1] + sp * v10[1]; 00108 v0[2] = s * v00[2] + sp * v10[2]; 00109 00110 v1[0] = s * v01[0] + sp * v11[0]; 00111 v1[1] = s * v01[1] + sp * v11[1]; 00112 v1[2] = s * v01[2] + sp * v11[2]; 00113 00114 for (dx = 1; dx < nxp; dx += 2) { 00115 t = 1. - (tp = dx / nx / 2.); 00116 v[0] = t * v0[0] + tp * v1[0]; 00117 v[1] = t * v0[1] + tp * v1[1]; 00118 v[2] = t * v0[2] + tp * v1[2]; 00119 00120 G3d_location2coord(map, v[0], v[1], v[2], &x, &y, &z); 00121 /* DEBUG 00122 printf ("(%d %d %d) (%lf %lf %lf) (%d %d %d) %lf\n", 00123 (int) dx / 2, (int) dy / 2, (int) dz / 2, 00124 v[0], v[1], v[2], 00125 x, y, z, 00126 G3d_getDoubleRegion (map, x, y, z)); 00127 */ 00128 if (type == DCELL_TYPE) 00129 *(doubleBuf + ((int)dz / 2) * nx * ny + 00130 ((int)dy / 2) * nx + (int)dx / 2) = 00131 G3d_getDoubleRegion(map, x, y, z); 00132 else 00133 *(floatBuf + ((int)dz / 2) * nx * ny + 00134 ((int)dy / 2) * nx + (int)dx / 2) = 00135 G3d_getFloatRegion(map, x, y, z); 00136 } 00137 } 00138 } 00139 } 00140 00141 /*---------------------------------------------------------------------------*/ 00142 00143 void 00144 G3d_getVolume(void *map, 00145 double originNorth, double originWest, double originBottom, 00146 double vxNorth, double vxWest, double vxBottom, 00147 double vyNorth, double vyWest, double vyBottom, 00148 double vzNorth, double vzWest, double vzBottom, 00149 int nx, int ny, int nz, void *volumeBuf, int type) 00150 { 00151 double u[2][2][2][3]; 00152 00153 u[0][0][0][0] = originNorth; 00154 u[0][0][0][1] = originWest; 00155 u[0][0][0][2] = originBottom; 00156 00157 u[0][0][1][0] = vxNorth; 00158 u[0][0][1][1] = vxWest; 00159 u[0][0][1][2] = vxBottom; 00160 00161 u[1][0][0][0] = vzNorth; 00162 u[1][0][0][1] = vzWest; 00163 u[1][0][0][2] = vzBottom; 00164 00165 u[1][0][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[1][0][0][0]; 00166 u[1][0][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[1][0][0][1]; 00167 u[1][0][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[1][0][0][2]; 00168 00169 u[0][1][0][0] = vyNorth; 00170 u[0][1][0][1] = vyWest; 00171 u[0][1][0][2] = vyBottom; 00172 00173 u[0][1][1][0] = (u[0][0][1][0] - u[0][0][0][0]) + u[0][1][0][0]; 00174 u[0][1][1][1] = (u[0][0][1][1] - u[0][0][0][1]) + u[0][1][0][1]; 00175 u[0][1][1][2] = (u[0][0][1][2] - u[0][0][0][2]) + u[0][1][0][2]; 00176 00177 u[1][1][0][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][0][0]; 00178 u[1][1][0][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][0][1]; 00179 u[1][1][0][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][0][2]; 00180 00181 u[1][1][1][0] = (u[1][0][0][0] - u[0][0][0][0]) + u[0][1][1][0]; 00182 u[1][1][1][1] = (u[1][0][0][1] - u[0][0][0][1]) + u[0][1][1][1]; 00183 u[1][1][1][2] = (u[1][0][0][2] - u[0][0][0][2]) + u[0][1][1][2]; 00184 00185 G3d_getVolumeA(map, u, nx, ny, nz, volumeBuf, type); 00186 } 00187 00188 /*---------------------------------------------------------------------------*/ 00189 00190 void 00191 G3d_getAlignedVolume(void *map, 00192 double originNorth, double originWest, 00193 double originBottom, double lengthNorth, 00194 double lengthWest, double lengthBottom, int nx, int ny, 00195 int nz, void *volumeBuf, int type) 00196 { 00197 G3d_getVolume(map, 00198 originNorth, originWest, originBottom, 00199 originNorth + lengthNorth, originWest, originBottom, 00200 originNorth, originWest + lengthWest, originBottom, 00201 originNorth, originWest, originBottom + lengthBottom, 00202 nx, ny, nz, volumeBuf, type); 00203 } 00204 00205 /*---------------------------------------------------------------------------*/ 00206 00207 void 00208 G3d_makeAlignedVolumeFile(void *map, const char *fileName, 00209 double originNorth, double originWest, 00210 double originBottom, double lengthNorth, 00211 double lengthWest, double lengthBottom, int nx, 00212 int ny, int nz) 00213 { 00214 void *volumeBuf; 00215 void *mapVolume; 00216 int x, y, z, eltLength; 00217 G3D_Region region; 00218 00219 volumeBuf = G3d_malloc(nx * ny * nz * sizeof(G3d_getFileType())); 00220 if (volumeBuf == NULL) 00221 G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_malloc"); 00222 00223 G3d_getAlignedVolume(map, 00224 originNorth, originWest, originBottom, 00225 lengthNorth, lengthWest, lengthBottom, 00226 nx, ny, nz, volumeBuf, G3d_getFileType()); 00227 00228 region.north = originNorth; 00229 region.south = originNorth + lengthNorth; 00230 region.east = originWest; 00231 region.west = originWest + lengthWest; 00232 region.top = originBottom; 00233 region.bottom = originBottom + lengthBottom; 00234 00235 region.rows = ny; 00236 region.cols = nx; 00237 region.depths = nz; 00238 00239 mapVolume = G3d_openCellNew(fileName, G3d_getFileType(), 00240 G3D_USE_CACHE_DEFAULT, ®ion); 00241 if (mapVolume == NULL) 00242 G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_openCellNew"); 00243 00244 eltLength = G3d_length(G3d_getFileType()); 00245 00246 for (z = 0; z < nz; z++) { 00247 for (y = 0; y < ny; y++) { 00248 for (x = 0; x < nx; x++) { 00249 /* G3d_putValueRegion? */ 00250 if (!G3d_putValue(mapVolume, x, y, z, 00251 G_incr_void_ptr(volumeBuf, 00252 (z * ny * nx + y * nx + 00253 x) * eltLength), 00254 G3d_fileTypeMap(mapVolume))) 00255 G3d_fatalError 00256 ("G3d_makeAlignedVolumeFile: error in G3d_putValue"); 00257 } 00258 } 00259 } 00260 00261 if (!G3d_closeCell(mapVolume)) 00262 G3d_fatalError("G3d_makeAlignedVolumeFile: error in G3d_closeCell"); 00263 00264 G3d_free(volumeBuf); 00265 }