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

edmplugin.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: edmplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.34 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <math.h>
00021 #include <string.h>
00022 
00023 #include "molfile_plugin.h"
00024 
00025 #ifndef NAN //not a number
00026   const float NAN = sqrtf(-1.f); //need some kind of portable NAN definition
00027 #endif
00028 
00029 #define TOLERANCE 1e-4
00030 
00031 typedef struct {
00032   FILE *fd;
00033   int nsets;
00034   molfile_volumetric_t *vol;
00035 } edm_t;
00036 
00037 static void eatline(FILE * fd) {
00038   char readbuf[1025];
00039   fgets(readbuf, 1024, fd);    // go on to next line
00040 }  
00041 
00042 static void *open_edm_read(const char *filepath, const char *filetype,
00043     int *natoms) {
00044   FILE *fd;
00045   edm_t *edm;
00046   int ntitle, na, nb, nc, xsize, ysize, zsize;
00047   int amin, amax, bmin, bmax, cmin, cmax;
00048   float a, b, c, alpha, beta, gamma;
00049   float xdelta, ydelta, zdelta;
00050   float alpha1, beta1, gamma1;
00051   float xaxis[3], yaxis[3], zaxis[3], z1, z2, z3;
00052   int i, convcnt;
00053   
00054   fd = fopen(filepath, "rb");
00055   if (!fd) 
00056     return NULL;
00057 
00058   edm = new edm_t;
00059   edm->fd = fd;
00060   edm->vol = NULL;
00061   *natoms = MOLFILE_NUMATOMS_NONE;
00062 
00063   edm->vol = new molfile_volumetric_t[1];
00064  
00065   edm->nsets = 1; // this EDM file contains only one data set
00066 
00067   // read in EDM file header information
00068   eatline(edm->fd);               // skip first header line 
00069 
00070   convcnt = fscanf(edm->fd, "%d", &ntitle); // read number of title lines
00071   if (convcnt != 1) {
00072     printf("edmplugin) failed to read in title line count\n");
00073     fclose(edm->fd);
00074     delete [] edm->vol;
00075     delete edm;
00076     return NULL;
00077   }
00078     
00079   eatline(edm->fd);               // go on to next line
00080 
00081   // skip past title and comment lines in header
00082   for (i=0; i<ntitle; i++) {
00083     eatline(edm->fd);             // skip a line
00084   }
00085 
00086   // read in the box dimensions and grid spacing deltas
00087   convcnt = fscanf(edm->fd, "%d %d %d %d %d %d %d %d %d",
00088          &na, &amin, &amax, &nb, &bmin, &bmax, &nc, &cmin, &cmax);
00089   if (convcnt != 9) {
00090     printf("edmplugin) failed to read in box dimensions\n");
00091     fclose(edm->fd);
00092     delete [] edm->vol;
00093     delete edm;
00094     return NULL;
00095   }
00096 
00097   eatline(edm->fd);               // go on to next line
00098   
00099   // calculate number of samples in each dimension
00100   xsize = amax - amin + 1;    
00101   ysize = bmax - bmin + 1;    
00102   zsize = cmax - cmin + 1;    
00103   edm->vol[0].xsize = xsize;
00104   edm->vol[0].ysize = ysize;
00105   edm->vol[0].zsize = zsize;
00106   edm->vol[0].has_color = 0;
00107 
00108   // read in 6 values for unit cell box orientation 
00109   convcnt = fscanf(edm->fd, "%f %f %f %f %f %f", 
00110                    &a, &b, &c, &alpha, &beta, &gamma);
00111   if (convcnt != 6) {
00112     printf("edmplugin) failed to read in box lengths and angles\n");
00113     fclose(edm->fd);
00114     delete [] edm->vol;
00115     delete edm;
00116     return NULL;
00117   }
00118   eatline(edm->fd);            // go on to next line
00119 
00120   // find box coordinates 
00121   xdelta = a / (float) na;
00122   ydelta = b / (float) nb;
00123   zdelta = c / (float) nc;
00124 
00125   strcpy(edm->vol[0].dataname, "X-PLOR Electron Density Map");
00126 
00127   // convert degrees to radians
00128   alpha1 = 3.14159265358979323846 * alpha / 180.0;
00129   beta1  = 3.14159265358979323846 *  beta / 180.0;
00130   gamma1 = 3.14159265358979323846 * gamma / 180.0;
00131 
00132   // calculate non-orthogonal unit cell coordinates
00133   xaxis[0] = xdelta;
00134   xaxis[1] = 0;
00135   xaxis[2] = 0;
00136 
00137   yaxis[0] = cos(gamma1) * ydelta;
00138   yaxis[1] = sin(gamma1) * ydelta;
00139   yaxis[2] = 0;
00140 
00141   z1 = cos(beta1);
00142   z2 = (cos(alpha1) - cos(beta1)*cos(gamma1)) / sin(gamma1);
00143   z3 = sqrt(1.0 - z1*z1 - z2*z2);
00144   zaxis[0] = z1 * zdelta;
00145   zaxis[1] = z2 * zdelta;
00146   zaxis[2] = z3 * zdelta;
00147 
00148   edm->vol[0].origin[0] = xaxis[0] * amin + yaxis[0] * bmin + zaxis[0] * cmin;
00149   edm->vol[0].origin[1] = yaxis[1] * bmin + zaxis[1] * cmin;
00150   edm->vol[0].origin[2] = zaxis[2] * cmin;
00151 
00152   edm->vol[0].xaxis[0] = xaxis[0] * (xsize-1);
00153   edm->vol[0].xaxis[1] = 0;
00154   edm->vol[0].xaxis[2] = 0;
00155 
00156   edm->vol[0].yaxis[0] = yaxis[0] * (ysize-1);
00157   edm->vol[0].yaxis[1] = yaxis[1] * (ysize-1);
00158   edm->vol[0].yaxis[2] = 0;
00159 
00160   edm->vol[0].zaxis[0] = zaxis[0] * (zsize-1);
00161   edm->vol[0].zaxis[1] = zaxis[1] * (zsize-1);
00162   edm->vol[0].zaxis[2] = zaxis[2] * (zsize-1);
00163 
00164   // Check that the EDM file is stored in the "ZYX" format we expect,
00165   // and return NULL if it is not a supported file type.
00166   char planeorder[4];
00167   memset(planeorder, 0, sizeof(planeorder));
00168   convcnt = fscanf(edm->fd, "%3s", planeorder);
00169   if (convcnt != 1) {
00170     printf("edmplugin) failed to read in plane order\n");
00171     fclose(edm->fd);
00172     delete [] edm->vol;
00173     delete edm;
00174     return NULL;
00175   }
00176 
00177   if (strcmp(planeorder, "ZYX")) { 
00178     printf("edmplugin) unsupported plane ordering %s\n", planeorder);
00179     fclose(edm->fd);
00180     delete [] edm->vol;
00181     delete edm;
00182     return NULL;
00183   }
00184   eatline(edm->fd);               // go on to next line
00185 
00186   return edm;
00187 }
00188 
00189 static int read_edm_metadata(void *v, int *nsets, 
00190   molfile_volumetric_t **metadata) {
00191   edm_t *edm = (edm_t *)v;
00192   *nsets = edm->nsets; 
00193   *metadata = edm->vol;  
00194 
00195   return MOLFILE_SUCCESS;
00196 }
00197 
00198 static int read_edm_data(void *v, int set, float *datablock,
00199                          float *colorblock) {
00200   edm_t *edm = (edm_t *)v;
00201   float * cell = datablock;
00202   int z, sentinel, convcnt;
00203   char readbuf[16];
00204  
00205   int xsize = edm->vol[0].xsize;
00206   int ysize = edm->vol[0].ysize;
00207   int zsize = edm->vol[0].zsize;
00208 
00209   // number of lines of text per slice
00210   int nperslice = xsize * ysize;
00211   int nlines = (int)(nperslice / 6.0);
00212   if ((nlines * 6) < nperslice) {
00213     nlines++;
00214   }
00215 
00216   for (z=0; z<zsize; z++) {
00217     int c;
00218     eatline(edm->fd);                // Eat the Z-plane index and throw away
00219 
00220     // read one plane of data, one cell at a time
00221     for (c=0; c<nperslice; c++) {
00222       convcnt = fscanf(edm->fd, "%f", cell);
00223       if (convcnt != 1) {
00224         printf("edmplugin) failed reading cell data\n");
00225         printf("edmplugin) cell %d of %d, slice %d\n", c, nperslice, z);
00226         return MOLFILE_ERROR; // bad file format encountered
00227       }
00228       cell++;
00229     }
00230     eatline(edm->fd);                // go on to next line
00231   }
00232 
00233   // read the -9999 end-of-file sentinel record
00234   sentinel = 0;
00235   fgets(readbuf, 13, edm->fd);
00236   sscanf(readbuf, "%d", &sentinel);  
00237   if (sentinel != -9999) {
00238     printf("edmplugin) EOF sentinel != -9999\n");
00239     // return MOLFILE_ERROR; // bad file format encountered, no sentinel record
00240   }
00241  
00242   return MOLFILE_SUCCESS;
00243 }
00244 
00245 static void close_edm_read(void *v) {
00246   edm_t *edm = (edm_t *)v;
00247 
00248   fclose(edm->fd);
00249   delete [] edm->vol; 
00250   delete edm;
00251 }
00252 
00253 static void *open_edm_write(const char *filepath, const char *filetype, int natoms) {
00254 
00255   FILE *fd = fopen(filepath, "w");
00256   if (!fd) {
00257     fprintf(stderr, "edmplugin) Could not open path '%s' for writing.\n",
00258       filepath);
00259   }
00260   return fd;
00261 }
00262 
00263 static void close_edm_write(void *v) {
00264   if (v) {
00265     fclose((FILE *)v);
00266   }
00267 }
00268 
00269 // XXX - The following interpolation code is duplicated from situsplugin.C
00270 //       (which is in turn duplicated from VMD). The only difference is that
00271 //       we are returning zeroes instead of NANs for out-of-range queries
00272 //       in the code below. These routines should be made available from 
00273 //       a centralized place to all molfile plugins, eliminating code 
00274 //       duplication.
00275 
00277 float edm_voxel_value_safe(int x, int y, int z, const int xsize, const int ysize, const int zsize, const float *data) {
00278   int xx, yy, zz; 
00279   xx = (x > 0) ? ((x < xsize) ? x : xsize-1) : 0;
00280   yy = (y > 0) ? ((y < ysize) ? y : ysize-1) : 0;
00281   zz = (z > 0) ? ((z < zsize) ? z : zsize-1) : 0;
00282   int index = zz*xsize*ysize + yy*xsize + xx;
00283   return data[index];
00284 }
00285 
00287 float edm_voxel_value_interpolate(float xv, float yv, float zv, const int xsize, const int ysize, const int zsize, const float *data) {
00288   int x = (int) xv;
00289   int y = (int) yv;
00290   int z = (int) zv;
00291   // fractional offset
00292   float xf = xv - x;
00293   float yf = yv - y;
00294   float zf = zv - z;
00295   float xlerps[4];
00296   float ylerps[2];
00297   float tmp;
00298 
00299   tmp = edm_voxel_value_safe(x, y, z, xsize, ysize, zsize, data);
00300   xlerps[0] = tmp + xf*(edm_voxel_value_safe(x+1, y, z, xsize, ysize, zsize, data) - tmp);
00301 
00302   tmp = edm_voxel_value_safe(x, y+1, z, xsize, ysize, zsize, data);
00303   xlerps[1] = tmp + xf*(edm_voxel_value_safe(x+1, y+1, z, xsize, ysize, zsize, data) - tmp);
00304 
00305   tmp = edm_voxel_value_safe(x, y, z+1, xsize, ysize, zsize, data);
00306   xlerps[2] = tmp + xf*(edm_voxel_value_safe(x+1, y, z+1, xsize, ysize, zsize, data) - tmp);
00307 
00308   tmp = edm_voxel_value_safe(x, y+1, z+1, xsize, ysize, zsize, data);
00309   xlerps[3] = tmp + xf*(edm_voxel_value_safe(x+1, y+1, z+1, xsize, ysize, zsize, data) - tmp);
00310 
00311   ylerps[0] = xlerps[0] + yf*(xlerps[1] - xlerps[0]);
00312   ylerps[1] = xlerps[2] + yf*(xlerps[3] - xlerps[2]);
00313 
00314   return ylerps[0] + zf*(ylerps[1] - ylerps[0]);
00315 }
00316 
00319 float edm_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) {
00320   xpos = (xpos-origin[0])/xdelta[0];
00321   ypos = (ypos-origin[1])/ydelta[1];
00322   zpos = (zpos-origin[2])/zdelta[2];
00323   int gx = (int) xpos; // XXX this is wrong for non-orthog cells.
00324   int gy = (int) ypos;
00325   int gz = (int) zpos;
00326 
00327 //  if (gx < 0 || gx >= xsize) return NAN;
00328 //  if (gy < 0 || gy >= ysize) return NAN;
00329 //  if (gz < 0 || gz >= zsize) return NAN;
00330 
00331   // Pad with zeroes
00332   if (gx < 0 || gx >= xsize) return 0;
00333   if (gy < 0 || gy >= ysize) return 0;
00334   if (gz < 0 || gz >= zsize) return 0;
00335     
00336   return edm_voxel_value_interpolate(xpos, ypos, zpos, xsize, ysize, zsize, data);
00337 
00338 }
00339 
00340 static int write_edm_data(void *v, molfile_volumetric_t *metadata, float *datablock, float *colorblock) {
00341 
00342   FILE *fd = (FILE *)v;
00343   const int xsize = metadata->xsize;
00344   const int ysize = metadata->ysize;
00345   const int zsize = metadata->zsize;
00346 
00347   float xaxis[3], yaxis[3], zaxis[3];
00348   float xdelta[3], ydelta[3], zdelta[3];
00349   float origin[3], porigin[3];
00350 
00351   int i, j, k;
00352 
00353   int na, amin, amax, nb, bmin, bmax, nc, cmin, cmax;
00354   float a, b, c, alpha, beta, gamma;
00355 
00356   for (i=0; i<3; i++) {
00357     origin[i] = metadata->origin[i];
00358     xaxis[i] = metadata->xaxis[i];
00359     yaxis[i] = metadata->yaxis[i];
00360     zaxis[i] = metadata->zaxis[i];
00361     xdelta[i] = xaxis[i]/(xsize - 1);
00362     ydelta[i] = yaxis[i]/(ysize - 1);
00363     zdelta[i] = zaxis[i]/(zsize - 1);
00364   }
00365 
00366   // The origin and the box length must be multiples of the grid spacing,
00367   // so we pad the box accordingly, keeping the same grid spacing and 
00368   // resampling the map.
00369 
00370   // For now, we will only write orthonogal cells
00371   if (fabs(xaxis[1]) > TOLERANCE || fabs(xaxis[2]) > TOLERANCE ||
00372       fabs(yaxis[0]) > TOLERANCE || fabs(yaxis[2]) > TOLERANCE ||
00373       fabs(zaxis[0]) > TOLERANCE || fabs(zaxis[1]) > TOLERANCE) {
00374     fprintf(stderr, "edmplugin) Could not write X-PLOR file: only orthogonal cells are currently supported.\n");
00375     return MOLFILE_ERROR;
00376   }
00377 
00378   amin = (int) floorf(origin[0]/xdelta[0]);
00379   bmin = (int) floorf(origin[1]/ydelta[1]);
00380   cmin = (int) floorf(origin[2]/zdelta[2]);
00381 
00382   porigin[0] = amin * xdelta[0];
00383   porigin[1] = bmin * ydelta[1];
00384   porigin[2] = cmin * zdelta[2];
00385 
00386   amax = (int) ceilf((xaxis[0]+origin[0])/xdelta[0]);
00387   bmax = (int) ceilf((yaxis[1]+origin[1])/ydelta[1]);
00388   cmax = (int) ceilf((zaxis[2]+origin[2])/zdelta[2]);
00389 
00390   na = amax - amin + 1;
00391   nb = bmax - bmin + 1;
00392   nc = cmax - cmin + 1;
00393 
00394   // The new cell axes may be slightly larger than the original ones
00395   a = na*xdelta[0];
00396   b = nb*ydelta[1];
00397   c = nc*zdelta[2];
00398 
00399   // Assuming the cell is orthognal
00400   alpha = beta = gamma = 90;
00401 
00402   // Write header
00403   fprintf(fd,"\n 2 !NTITLE\n"); // number of title lines
00404   fprintf(fd,"REMARKS FILENAME=\"\"\n");
00405   fprintf(fd,"REMARKS created by VMD\n");
00406 
00407   // Write the box dimensions and grid spacing deltas
00408   fprintf(fd,"%d %d %d %d %d %d %d %d %d\n", na, amin, amax, nb, bmin, bmax,
00409                                              nc, cmin, cmax);
00410   fprintf(fd,"%g %g %g %g %g %g\n", a, b, c, alpha, beta, gamma);
00411 
00412   // Write plane order
00413   fprintf(fd,"ZYX\n");
00414 
00415   // Copy voldata to a padded array
00416   int psize = na*nb*nc;
00417   int pxysize = na*nb;
00418 
00419   // Resample map
00420   float *pdata = (float*) malloc(psize*sizeof(float));
00421   for (i=0; i<na; i++) {
00422     float xpos = porigin[0] + i*xdelta[0];
00423     for (j=0; j<nb; j++) {
00424       float ypos = porigin[1] + j*ydelta[1];
00425       for (k=0; k<nc; k++) {
00426         float zpos = porigin[2] + k*zdelta[2];
00427         pdata[i + j*na + k*pxysize] = edm_voxel_value_interpolate_from_coord(xpos, ypos, zpos, origin, xdelta, ydelta, zdelta, xsize, ysize, zsize, datablock);
00428       }
00429     }
00430   }
00431 
00432   // Write each xy slice separately
00433   int count = 0;
00434   for (k=0; k<nc; k++) {
00435     if (count % 6) fprintf(fd, "\n");
00436     fprintf(fd, "%8d\n", k);
00437     count=0;
00438     for (j=0; j<nb; j++) {
00439       for (i=0; i<na; i++) {
00440         fprintf(fd, "%12.5e", pdata[k*pxysize + j*na + i]);
00441         if (++count % 6 == 0)
00442           fprintf(fd, "\n");
00443       }
00444     }
00445   }
00446   if (count % 6) fprintf(fd, "\n");
00447 
00448   // Write footer
00449   int sentinel = -9999;
00450   fprintf(fd, "%8d\n", sentinel);
00451 
00452   // Calculate average and standard deviation
00453   double avg = 0;
00454   double stddev = 0;
00455   double sum = 0;
00456   double sum2 = 0;
00457   for (i=0; i<psize; i++) {
00458     sum += pdata[i];
00459     sum2 += pdata[i]*pdata[i];
00460   }
00461   avg = sum/(double)psize;
00462   stddev = psize/(psize-1) * sqrt(sum2/psize - avg*avg);
00463   fprintf(fd, "%g %g\n", avg, stddev);
00464 
00465   free(pdata);
00466 
00467   fflush(fd);
00468   return MOLFILE_SUCCESS;
00469 }
00470 
00471 /*
00472  * Initialization stuff here
00473  */
00474 static molfile_plugin_t plugin;
00475 
00476 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00477   memset(&plugin, 0, sizeof(molfile_plugin_t));
00478   plugin.abiversion = vmdplugin_ABIVERSION;
00479   plugin.type = MOLFILE_PLUGIN_TYPE;
00480   plugin.name = "edm";
00481   plugin.prettyname = "XPLOR Electron Density Map";
00482   plugin.author = "John Stone, Leonardo Trabuco";
00483   plugin.majorv = 0;
00484   plugin.minorv = 9;
00485   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00486   plugin.filename_extension = "cns,edm,xplor";
00487   plugin.open_file_read = open_edm_read;
00488   plugin.read_volumetric_metadata = read_edm_metadata;
00489   plugin.read_volumetric_data = read_edm_data;
00490   plugin.close_file_read = close_edm_read;
00491 #if vmdplugin_ABIVERSION > 9
00492   plugin.open_file_write = open_edm_write;
00493   plugin.write_volumetric_data = write_edm_data;
00494   plugin.close_file_write = close_edm_write;
00495 #endif
00496   return VMDPLUGIN_SUCCESS; 
00497 }
00498 
00499 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00500   (*cb)(v, (vmdplugin_t *)&plugin);
00501   return VMDPLUGIN_SUCCESS;
00502 }
00503 
00504 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00505 
00506 
00507 
00508 #ifdef TEST_EDMPLUGIN
00509 
00510 int main(int argc, char *argv[]) {
00511   int natoms;
00512   void *v;
00513   int i, nsets, set;
00514   molfile_volumetric_t * meta;
00515 
00516   while (--argc) {
00517     ++argv;
00518     v = open_edm_read(*argv, "edm", &natoms);
00519     if (!v) {
00520       fprintf(stderr, "open_edm_read failed for file %s\n", *argv);
00521       return 1;
00522     }
00523 
00524     // try loading the EDM metadata now
00525     if (read_edm_metadata(v, &nsets, &meta)) {
00526       return 1; // failed to load edm file
00527     }
00528 
00529     for (set=0; set<nsets; set++) {
00530       printf("Loading volume set: %d\n", set);   
00531       
00532       int elements = meta[set].xsize * meta[set].ysize * meta[set].zsize;
00533       printf("   Grid Elements: %d\n", elements);
00534       printf(" Grid dimensions: X: %d Y: %d Z: %d\n", 
00535              meta[set].xsize, meta[set].ysize, meta[set].zsize);
00536 
00537       float * voldata = (float *) malloc(sizeof(float) * elements);
00538       float * coldata = NULL;
00539 
00540       if (meta[set].has_color) {
00541         coldata = (float *) malloc(sizeof(float) * elements * 3);
00542       }
00543 
00544       // try loading the EDM data sets now
00545       if (read_edm_data(v, set, voldata, coldata)) {
00546         return 1; // failed to load edm file
00547       }
00548 
00549       printf("First 6 elements:\n   ");
00550       for (i=0; i<6; i++) {
00551         printf("%g, ", voldata[i]);
00552       }
00553       printf("\n"); 
00554 
00555       printf("Last 6 elements:\n   ");
00556       for (i=elements - 6; i<elements; i++) {
00557         printf("%g, ", voldata[i]);
00558       }
00559       printf("\n"); 
00560     }
00561 
00562     close_edm_read(v);
00563   }
00564   return 0;
00565 }
00566 
00567 #endif
00568 
00569 
00570 
00571 

Generated on Fri Mar 29 03:10:05 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002