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

situsplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 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: situsplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.16 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * Situs EM map file reader
00020  *   map is a cubic lattice
00021  *   all data is stored in plain ASCII for convenience
00022  *
00023  * Format of the file is:
00024  *   voxel size in Angstroms
00025  *   coordinates of first voxel
00026  *   integer X/Y/Z counts
00027  *   voxels follow in X fastest, Y next fastest, Z slowest
00028  */
00029 
00030 #include <stdlib.h>
00031 #include <stdio.h>
00032 #include <math.h>
00033 #include <string.h>
00034 
00035 #if defined(_AIX)
00036 #include <strings.h>
00037 #endif
00038 
00039 #define MIN(X,Y) (((X)<(Y))? (X) : (Y))
00040 #define MAX(X,Y) (((X)>(Y))? (X) : (Y))
00041 
00042 #define TOLERANCE 1e-4
00043 
00044 #ifndef NAN //not a number
00045   const float NAN = sqrtf(-1.f); //need some kind of portable NAN definition
00046 #endif
00047 
00048 #include "molfile_plugin.h"
00049 
00050 typedef struct {
00051   FILE *fd;
00052   int nsets;
00053   molfile_volumetric_t *vol;
00054 } situs_t;
00055 
00056 
00057 static void *open_situs_read(const char *filepath, const char *filetype,
00058     int *natoms) {
00059   FILE *fd;
00060   situs_t *situs;
00061   float scale;
00062   int xsize, ysize, zsize;
00063   float orig[3]; 
00064   
00065   fd = fopen(filepath, "r");
00066   if (!fd) {
00067     printf("situsplugin) Error opening file.\n");
00068     return NULL;
00069   }
00070 
00071   /* get the scale */
00072   if (fscanf(fd, "%f", &scale) != 1) {;
00073     printf("situsplugin) Error reading voxel scale.\n");
00074     return NULL;
00075   }
00076 
00077   if (fscanf(fd, "%f %f %f", orig, orig+1, orig+2) != 3) {
00078     printf("situsplugin) Error reading grid origin.\n");
00079     return NULL;
00080   }
00081 
00082   /* get the number of grid points */
00083   if (fscanf(fd, "%d %d %d", &xsize, &ysize, &zsize) != 3) {
00084     printf("situsplugin) Error reading grid dimensions.\n");
00085     return NULL;
00086   }
00087 
00088   /* allocate and initialize the situs structure */
00089   situs = new situs_t;
00090   situs->fd = fd;
00091   situs->vol = NULL;
00092   *natoms = MOLFILE_NUMATOMS_NONE;
00093   situs->nsets = 1; /* this file contains only one data set */
00094 
00095   situs->vol = new molfile_volumetric_t[1];
00096   strcpy(situs->vol[0].dataname, "Situs map");
00097 
00098   /* Set the unit cell origin and basis vectors */
00099   for (int i=0; i<3; i++) {
00100     situs->vol[0].origin[i] = orig[i];
00101     situs->vol[0].xaxis[i] = 0.0;
00102     situs->vol[0].yaxis[i] = 0.0;
00103     situs->vol[0].zaxis[i] = 0.0;
00104   }
00105   situs->vol[0].xaxis[0] = scale * (xsize-1);
00106   situs->vol[0].yaxis[1] = scale * (ysize-1);
00107   situs->vol[0].zaxis[2] = scale * (zsize-1);
00108 
00109   situs->vol[0].xsize = xsize;
00110   situs->vol[0].ysize = ysize;
00111   situs->vol[0].zsize = zsize;
00112 
00113   situs->vol[0].has_color = 0; /* Situs maps contain no color information. */
00114 
00115   return situs;
00116 }
00117 
00118 static int read_situs_metadata(void *v, int *nsets, 
00119   molfile_volumetric_t **metadata) {
00120   situs_t *situs = (situs_t *)v;
00121   *nsets = situs->nsets; 
00122   *metadata = situs->vol;  
00123 
00124   return MOLFILE_SUCCESS;
00125 }
00126 
00127 static int read_situs_data(void *v, int set, float *datablock,
00128                          float *colorblock) {
00129   situs_t *situs = (situs_t *)v;
00130   FILE *fd = situs->fd;
00131   int xsize, ysize, zsize, xysize, count, total;
00132 
00133   xsize = situs->vol[0].xsize;
00134   ysize = situs->vol[0].ysize;
00135   zsize = situs->vol[0].zsize;
00136   xysize = xsize * ysize;
00137   total = xysize * zsize;
00138 
00139   /* Read the values from the file */
00140   for (count=0; count<total; count++) {
00141     if (fscanf(fd, "%f", datablock + count) != 1) {
00142       printf("situsplugin) Failed reading situs map data\n");
00143       return MOLFILE_ERROR;
00144     }
00145   }   
00146 
00147   return MOLFILE_SUCCESS;
00148 }
00149 
00150 static void close_situs_read(void *v) {
00151   situs_t *situs = (situs_t *)v;
00152   
00153   fclose(situs->fd);
00154   if (situs->vol != NULL)
00155     delete [] situs->vol; 
00156   delete situs;
00157 }
00158 
00159 static void *open_situs_write(const char *filepath, const char *filetype, int natoms) {
00160 
00161   FILE *fd = fopen(filepath, "w");
00162   if (!fd) {
00163     fprintf(stderr, "situsplugin) Could not open path '%s' for writing.\n",
00164       filepath);
00165   }
00166   return fd;
00167 }
00168 
00169 static void close_situs_write(void *v) {
00170   if (v) {
00171     fclose((FILE *)v);
00172   }
00173 }
00174 
00175 /* The Situs format requires the same grid spacing in all directions.
00176  * The following code allows us to re-sample the map so that we can 
00177  * write a Situs file when the uniform grid spacing requirement is not met.
00178  * */
00179 
00181 float situs_voxel_value_safe(int x, int y, int z, const int xsize, const int ysize, const int zsize, const float *data) {
00182   int xx, yy, zz; 
00183   xx = (x > 0) ? ((x < xsize) ? x : xsize-1) : 0;
00184   yy = (y > 0) ? ((y < ysize) ? y : ysize-1) : 0;
00185   zz = (z > 0) ? ((z < zsize) ? z : zsize-1) : 0;
00186   int index = zz*xsize*ysize + yy*xsize + xx;
00187   return data[index];
00188 }
00189 
00191 float situs_voxel_value_interpolate(float xv, float yv, float zv, const int xsize, const int ysize, const int zsize, const float *data) {
00192   int x = (int) xv;
00193   int y = (int) yv;
00194   int z = (int) zv;
00195   // fractional offset
00196   float xf = xv - x;
00197   float yf = yv - y;
00198   float zf = zv - z;
00199   float xlerps[4];
00200   float ylerps[2];
00201   float tmp;
00202 
00203   tmp = situs_voxel_value_safe(x, y, z, xsize, ysize, zsize, data);
00204   xlerps[0] = tmp + xf*(situs_voxel_value_safe(x+1, y, z, xsize, ysize, zsize, data) - tmp);
00205 
00206   tmp = situs_voxel_value_safe(x, y+1, z, xsize, ysize, zsize, data);
00207   xlerps[1] = tmp + xf*(situs_voxel_value_safe(x+1, y+1, z, xsize, ysize, zsize, data) - tmp);
00208 
00209   tmp = situs_voxel_value_safe(x, y, z+1, xsize, ysize, zsize, data);
00210   xlerps[2] = tmp + xf*(situs_voxel_value_safe(x+1, y, z+1, xsize, ysize, zsize, data) - tmp);
00211 
00212   tmp = situs_voxel_value_safe(x, y+1, z+1, xsize, ysize, zsize, data);
00213   xlerps[3] = tmp + xf*(situs_voxel_value_safe(x+1, y+1, z+1, xsize, ysize, zsize, data) - tmp);
00214 
00215   ylerps[0] = xlerps[0] + yf*(xlerps[1] - xlerps[0]);
00216   ylerps[1] = xlerps[2] + yf*(xlerps[3] - xlerps[2]);
00217 
00218   return ylerps[0] + zf*(ylerps[1] - ylerps[0]);
00219 }
00220 
00223 float situs_voxel_value_interpolate_from_coord(float xpos, float ypos, float zpos, const float *origin, const float *xdelta, const float *ydelta, const float *zdelta, const int xsize, const int ysize, const int zsize, float *data) {
00224   xpos = (xpos-origin[0])/xdelta[0];
00225   ypos = (ypos-origin[1])/ydelta[1];
00226   zpos = (zpos-origin[2])/zdelta[2];
00227   int gx = (int) xpos; // XXX this is wrong for non-orthog cells.
00228   int gy = (int) ypos;
00229   int gz = (int) zpos;
00230   if (gx < 0 || gx >= xsize) return NAN;
00231   if (gy < 0 || gy >= ysize) return NAN;
00232   if (gz < 0 || gz >= zsize) return NAN;
00233     
00234   return situs_voxel_value_interpolate(xpos, ypos, zpos, xsize, ysize, zsize, data);
00235 
00236 }
00237 
00238 static int write_situs_data(void *v, molfile_volumetric_t *metadata, float *datablock, float *colorblock) {
00239 
00240   FILE *fd = (FILE *)v;
00241   const int xsize = metadata->xsize;
00242   const int ysize = metadata->ysize;
00243   const int zsize = metadata->zsize;
00244   const int xysize = xsize * ysize;
00245 
00246   float xaxis[3], yaxis[3], zaxis[3];
00247   float xdelta[3], ydelta[3], zdelta[3];
00248   float origin[3];
00249 
00250   int i, j, k;
00251 
00252   for (i=0; i<3; i++) {
00253     origin[i] = metadata->origin[i];
00254     xaxis[i] = metadata->xaxis[i];
00255     yaxis[i] = metadata->yaxis[i];
00256     zaxis[i] = metadata->zaxis[i];
00257     xdelta[i] = xaxis[i]/(xsize - 1);
00258     ydelta[i] = yaxis[i]/(ysize - 1);
00259     zdelta[i] = zaxis[i]/(zsize - 1);
00260   }
00261 
00262   /* Situs format requires an orthogonal cell */
00263   if (fabs(xaxis[1]) > TOLERANCE || fabs(xaxis[2]) > TOLERANCE ||
00264       fabs(yaxis[0]) > TOLERANCE || fabs(yaxis[2]) > TOLERANCE ||
00265       fabs(zaxis[0]) > TOLERANCE || fabs(zaxis[1]) > TOLERANCE) {
00266     fprintf(stderr, "situsplugin) Could not write situs file: this format requires an orthogonal cell.\n");
00267     return MOLFILE_ERROR;
00268   }
00269 
00270   /* Situs format requires the same grid spacing in all dimensions */
00271   float xres = xdelta[0]*xdelta[0] + xdelta[1]*xdelta[1] + xdelta[2]*xdelta[2];
00272   float yres = ydelta[0]*ydelta[0] + ydelta[1]*ydelta[1] + ydelta[2]*ydelta[2];
00273   float zres = zdelta[0]*zdelta[0] + zdelta[1]*zdelta[1] + zdelta[2]*zdelta[2];
00274   if ( (fabs(xres-yres) > TOLERANCE) || (fabs(xres-zres) > TOLERANCE) ) {
00275 
00276     //fprintf(stderr, "situsplugin) Could not write situs file: this format requires the same grid spacing in all dimensions.\n");
00277     //return MOLFILE_ERROR;
00278     
00279     fprintf(stderr, "situsplugin) Warning: This format requires the same grid spacing in all dimensions. The map will be re-sampled to meet this requirement. The resulting cell may be slightly smaller than the original one.\n");
00280 
00281     // Use the smallest grid spacing for the re-sampled map
00282     float new_delta = MIN(xdelta[0], ydelta[1]);
00283     new_delta = MIN(new_delta, zdelta[2]);
00284 
00285     // Calculate the map dimensions for the new grid spacing 
00286     int new_xsize = int(xaxis[0]/new_delta);
00287     int new_ysize = int(yaxis[1]/new_delta);
00288     int new_zsize = int(zaxis[2]/new_delta);
00289     int new_xysize = new_xsize * new_ysize;
00290     int new_size = new_xsize * new_ysize * new_zsize;
00291 
00292     // Resample map
00293     float *new_data = (float *)malloc(3*new_size*sizeof(float));
00294     for (i=0; i<new_xsize; i++) {
00295       float xpos = origin[0] + i*new_delta;
00296       for (j=0; j<new_ysize; j++) {
00297         float ypos = origin[1] + j*new_delta;
00298         for (k=0; k<new_zsize; k++) {
00299           float zpos = origin[2] + k*new_delta;
00300           new_data[i + j*new_xsize + k*new_xysize] = situs_voxel_value_interpolate_from_coord(xpos, ypos, zpos, origin, xdelta, ydelta, zdelta, xsize, ysize, zsize, datablock);
00301         }
00302       }
00303     }
00304 
00305     /* Write a situs header */
00306     fprintf(fd, "%g %g %g %g %d %d %d\n\n", new_delta, origin[0], origin[1], origin[2], new_xsize, new_ysize, new_zsize);
00307   
00308     /* Write the main data array */
00309     int numentries = 1;
00310     for (k=0; k<new_zsize; k++) {
00311       for (j=0; j<new_ysize; j++) {
00312         for (i=0; i<new_xsize; i++) {
00313           fprintf(fd, "%g ", new_data[i + j*new_xsize + k*new_xysize]);
00314           if (numentries % 10 == 0) fprintf(fd, "\n");
00315           numentries++;
00316         }
00317       }
00318     }
00319 
00320     free(new_data);
00321 
00322   } else {
00323   
00324     /* Write a situs header */
00325     fprintf(fd, "%g %g %g %g %d %d %d\n\n", xdelta[0], origin[0], origin[1], origin[2], xsize, ysize, zsize);
00326   
00327     /* Write the main data array */
00328     int numentries = 1;
00329     for (k=0; k<zsize; k++) {
00330       for (j=0; j<ysize; j++) {
00331         for (i=0; i<xsize; i++) {
00332           fprintf(fd, "%g ", datablock[i + j*xsize + k*xysize]);
00333           if (numentries % 10 == 0) fprintf(fd, "\n");
00334           numentries++;
00335         }
00336       }
00337     }
00338 
00339   }
00340 
00341   fflush(fd);
00342   return MOLFILE_SUCCESS;
00343 }
00344 
00345 /*
00346  * Initialization stuff here
00347  */
00348 static molfile_plugin_t plugin;
00349 
00350 VMDPLUGIN_API int VMDPLUGIN_init(void) {
00351   memset(&plugin, 0, sizeof(molfile_plugin_t));
00352   plugin.abiversion = vmdplugin_ABIVERSION;
00353   plugin.type = MOLFILE_PLUGIN_TYPE;
00354   plugin.name = "situs";
00355   plugin.prettyname = "Situs Density Map";
00356   plugin.author = "John Stone, Leonardo Trabuco";
00357   plugin.majorv = 1;
00358   plugin.minorv = 5;
00359   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00360   plugin.filename_extension = "sit,situs";
00361   plugin.open_file_read = open_situs_read;
00362   plugin.read_volumetric_metadata = read_situs_metadata;
00363   plugin.read_volumetric_data = read_situs_data;
00364   plugin.close_file_read = close_situs_read;
00365 #if vmdplugin_ABIVERSION > 9
00366   plugin.open_file_write = open_situs_write;
00367   plugin.write_volumetric_data = write_situs_data;
00368   plugin.close_file_write = close_situs_write;
00369 #endif
00370   return VMDPLUGIN_SUCCESS;
00371 }
00372 
00373 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00374   (*cb)(v, (vmdplugin_t *)&plugin);
00375   return VMDPLUGIN_SUCCESS;
00376 }
00377 
00378 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00379 

Generated on Thu Mar 28 03:08:19 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002