Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

vtkplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 2007-2011 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: vtkplugin.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.5 $      $Date: 2015/10/26 22:14:45 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Plugin for reading uniform grids and vector fields written 
00019  *   in the VTK ASCII format
00020  ***************************************************************************/
00021 
00022 //
00023 // Prototype file reader developed using VTK refs and a couple samples:
00024 //   http://www.vtk.org/Wiki/VTK/Tutorials/3DDataTypes
00025 //   http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
00026 //
00027 
00028 // This parser currently expects the data to be formatted like this:
00029 // # vtk DataFile Version X.0
00030 // titletext
00031 // ASCII
00032 // DATASET STRUCTURED_POINTS
00033 // DIMENSIONS XXX YYY ZZZ
00034 // SPACING XX YY ZZ
00035 // ORIGIN XX YY XX
00036 // POINT_DATA n
00037 //
00038 // And either a field like this:
00039 // FIELD fieldname numarrays
00040 // arrayname0 numcomponents numtuples datatype
00041 // val0 val1 ...
00042 //
00043 // Or a list of vectors like this:
00044 // VECTORS dataname datatype
00045 // val0 val1 ...
00046 //
00047 
00048 #include <stdlib.h>
00049 #include <stdio.h>
00050 #include <ctype.h>
00051 #include <math.h>
00052 #include <string.h>
00053 
00054 #if defined(_AIX)
00055 #include <strings.h>
00056 #endif
00057 
00058 #if defined(WIN32) || defined(WIN64)
00059 #define strcasecmp  stricmp
00060 #define strncasecmp strnicmp
00061 #endif
00062 
00063 #include "molfile_plugin.h"
00064 #include "largefiles.h"
00065 
00066 #define THISPLUGIN plugin
00067 #include "vmdconio.h"
00068 
00069 #define LINESIZE 2040
00070 
00071 typedef struct {
00072   FILE *fd;
00073   char title[257]; 
00074   int nsets;
00075   molfile_volumetric_t *vol;
00076   int isBinary; 
00077 } vtk_t;
00078 
00079 
00080 // Get a string from a stream, printing any errors that occur
00081 static char *vtkgets(char *s, int n, FILE *stream) {
00082   char *returnVal;
00083 
00084   if (feof(stream)) {
00085     printf("vtkplugin) Unexpected end-of-file.\n");
00086     return NULL;
00087   } else if (ferror(stream)) {
00088     printf("vtkplugin) Error reading file.\n");
00089     return NULL;
00090   } else {
00091     returnVal = fgets(s, n, stream);
00092     if (returnVal == NULL) {
00093       printf("vtkplugin) Error reading line.\n");
00094     }
00095   }
00096 
00097   return returnVal;
00098 }
00099 
00100 
00101 static int vtkgetstrcmp(char *s, int n, FILE *stream, const char *cmpstr) {
00102   char *str = vtkgets(s, n, stream);
00103   int rc = strncmp(cmpstr, str, strlen(cmpstr));
00104   if (rc) {
00105     printf("vtkplugin) found '%s', expected '%s'\n", str, cmpstr);
00106   }
00107   return rc;
00108 }
00109    
00110 
00111 static void *open_vtk_read(const char *filepath, const char *filetype,
00112     int *natoms) {
00113   FILE *fd;
00114   vtk_t *vtk;
00115   char inbuf[LINESIZE];
00116   int xsize, ysize, zsize;
00117   float orig[3], xdelta[3], ydelta[3], zdelta[3];
00118  
00119   memset(orig, 0, sizeof(orig));
00120   memset(xdelta, 0, sizeof(xdelta));
00121   memset(ydelta, 0, sizeof(ydelta));
00122   memset(zdelta, 0, sizeof(zdelta));
00123  
00124   fd = fopen(filepath, "rb");
00125   if (!fd) {
00126     printf("vtkplugin) Error opening file.\n");
00127     return NULL;
00128   }
00129 
00130   // allocate and initialize the vtk structure
00131   vtk = new vtk_t;
00132   memset(vtk, 0, sizeof(vtk_t));
00133   vtk->fd = fd;
00134   vtk->vol = NULL;
00135   vtk->isBinary = 0;
00136   *natoms = MOLFILE_NUMATOMS_NONE;
00137   vtk->nsets = 1; /* this file contains only one data set */
00138 
00139   /* skip comments */
00140   do {
00141     if (vtkgets(inbuf, LINESIZE, fd) == NULL) 
00142       return NULL;
00143   } while (inbuf[0] == '#');
00144 
00145   // read VTK title line 
00146   printf("vtkplugin) Dataset title: '%s'\n", inbuf);
00147   strncpy(vtk->title, inbuf, sizeof(vtk->title) - 1);
00148   vtk->title[256]='\0'; // force-terminate potentially truncated title string
00149 
00150   if (vtkgetstrcmp(inbuf, LINESIZE, fd, "ASCII")) return NULL;
00151   if (vtkgetstrcmp(inbuf, LINESIZE, fd, "DATASET STRUCTURED_POINTS")) return NULL;
00152 
00153   // get the grid dimensions
00154   if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00155     delete vtk;
00156     return NULL;
00157   }
00158   if (sscanf(inbuf, "DIMENSIONS %d %d %d", &xsize, &ysize, &zsize) != 3) {
00159     printf("vtkplugin) Error reading grid dimensions!\n");
00160     delete vtk;
00161     return NULL;
00162   }
00163 
00164   // get the grid spacing
00165   if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00166     delete vtk;
00167     return NULL;
00168   }
00169   if (sscanf(inbuf, "SPACING %e %e %e", xdelta, ydelta+1, zdelta+2) != 3) {
00170     printf("vtkplugin) Error reading cell dimensions!\n");
00171     delete vtk;
00172     return NULL;
00173   }
00174 
00175   // get the grid origin
00176   if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00177     delete vtk;
00178     return NULL;
00179   }
00180   if (sscanf(inbuf, "ORIGIN %e %e %e", orig, orig+1, orig+2) != 3) {
00181     printf("vtkplugin) Error reading grid origin!\n");
00182     delete vtk;
00183     return NULL;
00184   }
00185 
00186   // get number of grid points
00187   if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00188     delete vtk;
00189     return NULL;
00190   }
00191   int numgridpoints = 0;
00192   if (sscanf(inbuf, "POINT_DATA %d", &numgridpoints) != 1) {
00193     printf("vtkplugin) Error reading grid point counts!\n");
00194     delete vtk;
00195     return NULL;
00196   }
00197 
00198   // get field or vector list depending on file format variant we have
00199   if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00200     delete vtk;
00201     return NULL;
00202   }
00203 
00204   char tmp[256];
00205   sscanf(inbuf, "%s", tmp);
00206   if (!strcmp(tmp, "FIELD")) {
00207     char fieldname[256];
00208     int numarrays=0;
00209     sscanf(inbuf, "FIELD %s %d", fieldname, &numarrays);
00210     printf("vtkplugin) FIELD: name '%s', %d arrays\n", fieldname, numarrays);
00211 
00212     // eat the array name for field zero
00213     if (vtkgets(inbuf, LINESIZE, fd) == NULL) {
00214       delete vtk;
00215       return NULL;
00216     }
00217   } else if (!strcmp(tmp, "VECTORS")) {
00218     // prepare to eat vectors
00219     char fieldname[256];
00220     int numvecs=0;
00221     sscanf(inbuf, "VECTORS %s %d", fieldname, &numvecs);
00222     printf("vtkplugin) VECTORS: name '%s', %d arrays\n", fieldname, numvecs);
00223   } else {
00224     printf("vtkplugin) Unrecognized file structure, aborting!:\n");
00225     printf("vtkplugin) line contents: '%s'\n", inbuf);
00226     delete vtk;
00227     return NULL;
00228   } 
00229 
00230   vtk->vol = new molfile_volumetric_t[1];
00231   memset(vtk->vol, 0, sizeof(molfile_volumetric_t));
00232   strcpy(vtk->vol[0].dataname, "VTK volumetric map");
00233 
00234   /* Set the unit cell origin and basis vectors */
00235   for (int i=0; i<3; i++) {
00236     vtk->vol[0].origin[i] = orig[i];
00237     vtk->vol[0].xaxis[i] = xdelta[i] * ((xsize-1 > 0) ? (xsize-1) : 1);
00238     vtk->vol[0].yaxis[i] = ydelta[i] * ((ysize-1 > 0) ? (ysize-1) : 1);
00239     vtk->vol[0].zaxis[i] = zdelta[i] * ((zsize-1 > 0) ? (zsize-1) : 1);
00240   }
00241 
00242   vtk->vol[0].xsize = xsize;
00243   vtk->vol[0].ysize = ysize;
00244   vtk->vol[0].zsize = zsize;
00245 
00246 #if vmdplugin_ABIVERSION > 16
00247   vtk->vol[0].has_scalar = 1;
00248   vtk->vol[0].has_gradient = 1;
00249   vtk->vol[0].has_variance = 0;
00250 #endif
00251   vtk->vol[0].has_color = 0; // no color data
00252 
00253   return vtk;
00254 }
00255 
00256 
00257 static int read_vtk_metadata(void *v, int *nsets, 
00258   molfile_volumetric_t **metadata) {
00259   vtk_t *vtk = (vtk_t *)v;
00260   *nsets = vtk->nsets; 
00261   *metadata = vtk->vol;  
00262 
00263   return MOLFILE_SUCCESS;
00264 }
00265 
00266 
00267 static int read_vtk_data(void *v, int set, float *datablock,
00268                          float *colorblock) {
00269   vtk_t *vtk = (vtk_t *)v;
00270   FILE *fd = vtk->fd;
00271   int x, y, z, xsize, ysize, zsize, xysize, total;
00272 
00273   if (vtk->isBinary)
00274     return MOLFILE_ERROR;
00275 
00276   xsize = vtk->vol[0].xsize;
00277   ysize = vtk->vol[0].ysize;
00278   zsize = vtk->vol[0].zsize;
00279   xysize = xsize * ysize;
00280   total = xysize * zsize;
00281 
00282   double scalemag = 1.0;
00283   const char *userscalefactor=getenv("VMDVTKPLUGINSCALEVOXELMAG");
00284   if (userscalefactor) {
00285     scalemag = atof(userscalefactor);
00286     if (scalemag != 0.0) {
00287       printf("vtkplugin) Applying user scaling factor to voxel scalar/gradient values: %g\n", scalemag); 
00288     } else {
00289       printf("vtkplugin) Warning: ignoring user scaling factor due to parse error or zero-value\n"); 
00290     }
00291   } else {
00292     printf("vtkplugin) No user scaling factor set, using scale factor 1.0.\n");
00293   }
00294 
00295   float maxmag = 0.0f;
00296   strcpy(vtk->vol[0].dataname, "volgradient");
00297   for (z=0; z<zsize; z++) {
00298     for (y=0; y<ysize; y++) {
00299       for (x=0; x<xsize; x++) {
00300         double vx, vy, vz;
00301         fscanf(fd, "%lf %lf %lf", &vx, &vy, &vz);
00302 
00303 #if 1
00304         // XXX hack to allow user override of voxel magnitude
00305         // during loading...
00306         vx *= scalemag;
00307         vy *= scalemag;
00308         vz *= scalemag;
00309 #endif
00310 
00311         // compute scalar magnitude from vector field
00312         double mag = sqrt(vx*vx + vy*vy + vz*vz);
00313 
00314         int addr = z*xsize*ysize + y*xsize + x;
00315         datablock[addr] = mag;
00316 
00317         if (mag > maxmag)
00318           maxmag = mag;
00319       }
00320     }
00321   }
00322   printf("vtkplugin) maxmag: %g\n", maxmag);
00323 
00324   return MOLFILE_SUCCESS;
00325 }
00326 
00327 
00328 #if vmdplugin_ABIVERSION > 16
00329 
00330 static int read_vtk_data_ex(void *v, molfile_volumetric_readwrite_t *p) {
00331   vtk_t *vtk = (vtk_t *)v;
00332   FILE *fd = vtk->fd;
00333   int x, y, z, xsize, ysize, zsize, xysize, total;
00334 
00335   if (vtk->isBinary)
00336     return MOLFILE_ERROR;
00337 
00338   if (!p->scalar || !p->gradient) 
00339     return MOLFILE_ERROR;
00340 
00341   xsize = vtk->vol[0].xsize;
00342   ysize = vtk->vol[0].ysize;
00343   zsize = vtk->vol[0].zsize;
00344   xysize = xsize * ysize;
00345   total = xysize * zsize;
00346 
00347   double scalemag = 1.0;
00348   const char *userscalefactor=getenv("VMDVTKPLUGINSCALEVOXELMAG");
00349   if (userscalefactor) {
00350     scalemag = atof(userscalefactor);
00351     if (scalemag != 0.0) {
00352       printf("vtkplugin) Applying user scaling factor to voxel scalar/gradient values: %g\n", scalemag); 
00353     } else {
00354       printf("vtkplugin) Warning: ignoring user scaling factor due to parse error or zero-value\n"); 
00355     }
00356   } else {
00357     printf("vtkplugin) No user scaling factor set, using scale factor 1.0.\n");
00358   }
00359 
00360   float maxmag = 0.0f;
00361   strcpy(vtk->vol[0].dataname, "volgradient");
00362   for (z=0; z<zsize; z++) {
00363     for (y=0; y<ysize; y++) {
00364       for (x=0; x<xsize; x++) {
00365         double vx, vy, vz;
00366         fscanf(fd, "%lf %lf %lf", &vx, &vy, &vz);
00367 
00368 #if 1
00369         // XXX hack to allow user override of voxel magnitude
00370         // during loading...
00371         vx *= scalemag;
00372         vy *= scalemag;
00373         vz *= scalemag;
00374 #endif
00375 
00376         // compute scalar magnitude from vector field
00377         double mag = sqrt(vx*vx + vy*vy + vz*vz);
00378 
00379         int addr = z*xsize*ysize + y*xsize + x;
00380         p->scalar[addr] = mag;
00381 
00382         if (mag > maxmag)
00383           maxmag = mag;
00384 
00385         // assign vector field to gradient map
00386         // index into vector field of 3-component vectors
00387         int addr3 = addr *= 3;
00388         p->gradient[addr3    ] = vx;
00389         p->gradient[addr3 + 1] = vy;
00390         p->gradient[addr3 + 2] = vz;
00391       }
00392     }
00393   }
00394   printf("vtkplugin) maxmag: %g\n", maxmag);
00395 
00396   return MOLFILE_SUCCESS;
00397 }
00398 
00399 #endif
00400 
00401 
00402 static void close_vtk_read(void *v) {
00403   vtk_t *vtk = (vtk_t *)v;
00404   
00405   fclose(vtk->fd);
00406   if (vtk->vol != NULL)
00407     delete [] vtk->vol; 
00408   delete vtk;
00409 }
00410 
00411 
00412 //
00413 // Initialization stuff here
00414 //
00415 static molfile_plugin_t vtkplugin;
00416 
00417 VMDPLUGIN_API int VMDPLUGIN_init(void) {
00418   memset(&vtkplugin, 0, sizeof(molfile_plugin_t));
00419   vtkplugin.abiversion = vmdplugin_ABIVERSION;
00420   vtkplugin.type = MOLFILE_PLUGIN_TYPE;
00421   vtkplugin.name = "vtk";
00422   vtkplugin.prettyname = "VTK grid reader";
00423   vtkplugin.author = "John Stone";
00424   vtkplugin.majorv = 0;
00425   vtkplugin.minorv = 2;
00426   vtkplugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00427   vtkplugin.filename_extension = "vtk";
00428   vtkplugin.open_file_read = open_vtk_read;
00429   vtkplugin.read_volumetric_metadata = read_vtk_metadata;
00430   vtkplugin.read_volumetric_data = read_vtk_data;
00431 #if vmdplugin_ABIVERSION > 16
00432   vtkplugin.read_volumetric_data_ex = read_vtk_data_ex;
00433 #endif
00434   vtkplugin.close_file_read = close_vtk_read;
00435   return VMDPLUGIN_SUCCESS;
00436 }
00437 
00438 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00439   (*cb)(v, (vmdplugin_t *)&vtkplugin);
00440   return VMDPLUGIN_SUCCESS;
00441 }
00442 
00443 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00444 
00445 

Generated on Fri Jan 17 02:57:54 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002