GRASS Programmer's Manual  6.4.2(2012)
g3dvolume.c
Go to the documentation of this file.
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, &region);
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines