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

brixplugin.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: brixplugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.20 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * BRIX format electron density maps.
00020  *
00021  * More info for format can be found at 
00022  * <http://zombie.imsb.au.dk/~mok/brix/brix-1.html#ss1.2>
00023  * 
00024  * TODO: Speed up inner loop of read_data, perhaps reading one full brick
00025  *       from the file at a time.
00026  */
00027 
00028 #include <stdlib.h>
00029 #include <stdio.h>
00030 #include <ctype.h>
00031 #include <math.h>
00032 #include <string.h>
00033 
00034 #if defined(_AIX)
00035 #include <strings.h>
00036 #endif
00037 
00038 #if defined(WIN32) || defined(WIN64)
00039 #define strcasecmp stricmp
00040 #endif
00041 
00042 #ifndef M_PI
00043 #define M_PI 3.14159265358979323846
00044 #endif
00045 
00046 #include "molfile_plugin.h"
00047 
00048 typedef struct {
00049   FILE *fd;
00050   int nsets;
00051   float prod, plus;
00052   molfile_volumetric_t *vol;
00053 } brix_t;
00054 
00055 
00056 static void *open_brix_read(const char *filepath, const char *filetype,
00057     int *natoms) {
00058   FILE *fd;
00059   brix_t *brix;
00060   char keyWord[81];
00061   // File header data:
00062   int org_x, org_y, org_z, ext_x, ext_y, ext_z;
00063   float grid_x, grid_y, grid_z, cell_x, cell_y, cell_z, 
00064         cell_alpha, cell_beta, cell_gamma, prod, plus, sigma, 
00065         xaxis[3], yaxis[3], zaxis[3], z1, z2, z3;
00066   
00067   fd = fopen(filepath, "rb");
00068   if (!fd) {
00069     fprintf(stderr, "brixplugin) Error opening file.\n");
00070     return NULL;
00071   }
00072 
00073   // read in BRIX file header information -- stored as a 512-element char array
00074   fscanf(fd, "%3s", keyWord);
00075   if (strcmp(keyWord, ":-)") != 0) {
00076     fprintf(stderr, "brixplugin) Error improperly formatted header.\n");
00077     return NULL;
00078   }
00079 
00080   // "Origin: the origin of the map in grid units, along X, Y, Z."
00081   fscanf(fd, " %s %d %d %d", keyWord, &org_x, &org_y, &org_z);
00082   if (strcasecmp(keyWord, "origin") != 0) {
00083     fprintf(stderr, "brixplugin) Error reading origin.\n");
00084     return NULL;
00085   }
00086 
00087   // "Extent: the extent (size) of the map, along X, Y, Z, in grid units"
00088   fscanf(fd, " %s %d %d %d", keyWord, &ext_x, &ext_y, &ext_z);
00089   if (strcasecmp(keyWord, "extent") != 0) {
00090     fprintf(stderr, "brixplugin) Error reading extent.\n");
00091     return NULL;
00092   }
00093 
00094   // "Grid: number of grid points along the whole unit cell, along X, Y, Z."
00095   fscanf(fd, " %s %f %f %f", keyWord, &grid_x, &grid_y, &grid_z);
00096   if (strcasecmp(keyWord, "grid") != 0) {
00097     fprintf(stderr, "brixplugin) Error reading grid.\n");
00098     return NULL;
00099   }
00100 
00101   // "Cell: Unit cell parameters"
00102   // cell x, y, and z are the length of the unit cell along each axis,
00103   // cell alpha, beta, and cell_gamma are the angles between axes.
00104   fscanf(fd, " %s %f %f %f %f %f %f", keyWord,
00105          &cell_x, &cell_y, &cell_z, &cell_alpha, &cell_beta, &cell_gamma);
00106   if (strcasecmp(keyWord, "cell") != 0) {
00107     fprintf(stderr, "brixplugin) Error reading cell.\n");
00108     return NULL;
00109   }
00110 
00111   // Convert angles to radians
00112   cell_alpha *= M_PI / 180.0;
00113   cell_beta *= M_PI / 180.0;
00114   cell_gamma *= M_PI / 180.0;
00115 
00116   // "Prod, plus: Constants that bring the electron density from byte to 
00117   // normal scale."
00118   fscanf(fd, " %s %f", keyWord, &prod);
00119   if (strcasecmp(keyWord, "prod") != 0) {
00120     fprintf(stderr, "brixplugin) Error reading prod.\n");
00121     return NULL;
00122   }
00123   fscanf(fd, " %s %f", keyWord, &plus);
00124   if (strcasecmp(keyWord, "plus") != 0) {
00125     fprintf(stderr, "brixplugin) Error reading plus.\n");
00126     return NULL;
00127   }
00128 
00129   // "Sigma: Rms deviation of electron density map."
00130   fscanf(fd, " %s %f", keyWord, &sigma);
00131   if (strcasecmp(keyWord, "sigma") != 0) {
00132     fprintf(stderr, "brixplugin) Error reading sigma.\n");
00133     return NULL;
00134   }
00135 
00136   // Allocate and initialize the brix structure
00137   brix = new brix_t;
00138   brix->fd = fd;
00139   brix->vol = NULL;
00140   *natoms = MOLFILE_NUMATOMS_NONE;
00141   brix->nsets = 1; // this file contains only one data set
00142 
00143   brix->prod = prod;
00144   brix->plus = plus;
00145 
00146   brix->vol = new molfile_volumetric_t[1];
00147   strcpy(brix->vol[0].dataname, "BRIX Electron Density Map");
00148 
00149   // Calculate non-orthogonal unit cell coordinates
00150   xaxis[0] = cell_x / grid_x;
00151   xaxis[1] = 0;
00152   xaxis[2] = 0;
00153 
00154   yaxis[0] = cos(cell_gamma) * cell_y / grid_y;
00155   yaxis[1] = sin(cell_gamma) * cell_y / grid_y;
00156   yaxis[2] = 0;
00157 
00158   z1 = cos(cell_beta);
00159   z2 = (cos(cell_alpha) - cos(cell_beta)*cos(cell_gamma)) / sin(cell_gamma);
00160   z3 = sqrt(1.0 - z1*z1 - z2*z2);
00161   zaxis[0] = z1 * cell_z / grid_z;
00162   zaxis[1] = z2 * cell_z / grid_z;
00163   zaxis[2] = z3 * cell_z / grid_z;
00164 
00165   brix->vol[0].origin[0] = xaxis[0] * org_x + yaxis[0] * org_y + 
00166                            zaxis[0] * org_z;
00167   brix->vol[0].origin[1] = yaxis[1] * org_y + zaxis[1] * org_z;
00168   brix->vol[0].origin[2] = zaxis[2] * org_z;
00169 
00170   brix->vol[0].xaxis[0] = xaxis[0] * (ext_x - 1);
00171   brix->vol[0].xaxis[1] = 0;
00172   brix->vol[0].xaxis[2] = 0;
00173 
00174   brix->vol[0].yaxis[0] = yaxis[0] * (ext_y - 1);
00175   brix->vol[0].yaxis[1] = yaxis[1] * (ext_y - 1);
00176   brix->vol[0].yaxis[2] = 0;
00177 
00178   brix->vol[0].zaxis[0] = zaxis[0] * (ext_z - 1);
00179   brix->vol[0].zaxis[1] = zaxis[1] * (ext_z - 1);
00180   brix->vol[0].zaxis[2] = zaxis[2] * (ext_z - 1);
00181 
00182   brix->vol[0].xsize = ext_x;
00183   brix->vol[0].ysize = ext_y;
00184   brix->vol[0].zsize = ext_z;
00185 
00186   brix->vol[0].has_color = 0;
00187 
00188   return brix;
00189 }
00190 
00191 static int read_brix_metadata(void *v, int *nsets, 
00192   molfile_volumetric_t **metadata) {
00193   brix_t *brix = (brix_t *)v;
00194   *nsets = brix->nsets; 
00195   *metadata = brix->vol;  
00196 
00197   return MOLFILE_SUCCESS;
00198 }
00199 
00200 static int read_brix_data(void *v, int set, float *datablock,
00201                          float *colorblock) {
00202   brix_t *brix = (brix_t *)v;
00203   float * cell = datablock;
00204   int xsize, ysize, zsize, xysize, xbrix, ybrix, zbrix, cellIndex;
00205   int x, y, z, xbrik, ybrik, zbrik;
00206   unsigned char brick[512];
00207   FILE * fd = brix->fd;
00208   float div, plus;
00209 
00210   // Read 512-byte "bricks" of data. Each brick contains data for 8*8*8 
00211   // gridpoints.
00212   fseek(fd, 512, SEEK_SET);
00213 
00214   div = 1.0 / brix->prod;
00215   plus = brix->plus;
00216 
00217   xsize = brix->vol[0].xsize;
00218   ysize = brix->vol[0].ysize;
00219   zsize = brix->vol[0].zsize;
00220   xysize = xsize * ysize;
00221 
00222   xbrix = (int) ceil((float) xsize / 8.0);
00223   ybrix = (int) ceil((float) ysize / 8.0);
00224   zbrix = (int) ceil((float) zsize / 8.0);
00225 
00226   // For every density value, 
00227   // cellIndex = (x + xbrik*8) + (y + ybrik*8)*xsize + (z + zbrik*8) * xysize
00228   cellIndex = 0;
00229   for (zbrik = 0; zbrik < zbrix; zbrik++) {
00230     for (ybrik = 0; ybrik < ybrix; ybrik++) {
00231       for (xbrik = 0; xbrik < xbrix; xbrik++) {
00232         if (feof(fd)) {
00233           fprintf(stderr, "brixplugin) Unexpected end-of-file.\n");
00234           return MOLFILE_ERROR;
00235         }
00236         if (ferror(fd)) {
00237           fprintf(stderr, "brixplugin) Error reading file.\n");
00238           return MOLFILE_ERROR;
00239         }
00240 
00241         // Read a data brick into the buffer.
00242         fread(brick, sizeof(char), 512, fd);
00243 
00244         // Copy the brick values into the cell
00245         for (z = 0; z < 8; z++) {
00246           for (y = 0; y < 8; y++) {
00247             for (x = 0; x < 8; x++) {
00248 
00249               if ( ((x + xbrik*8) < xsize) && ((y + ybrik*8) < ysize) &&
00250                    ((z + zbrik*8) < zsize) )
00251                 cell[cellIndex] = 
00252                   div * ((float) brick[x+y*8+z*64] - plus);
00253 
00254               cellIndex++;
00255             } // end for(x)
00256  
00257             cellIndex += xsize - 8;
00258           } // end for(y)
00259 
00260           cellIndex += xysize - 8*xsize;
00261         } // end for(z)
00262  
00263         cellIndex += 8 - 8*xysize; 
00264       } // end for(xbrik)
00265 
00266       cellIndex += 8 * (xsize - xbrix);
00267     } // end for(ybrik)
00268 
00269     cellIndex += 8 * (xysize - xsize*ybrik);
00270   } // end for(zbrik)
00271  
00272   return MOLFILE_SUCCESS;
00273 }
00274 
00275 static void close_brix_read(void *v) {
00276   brix_t *brix = (brix_t *)v;
00277 
00278   fclose(brix->fd);
00279   if (brix->vol != NULL)
00280     delete [] brix->vol; 
00281   delete brix;
00282 }
00283 
00284 /*
00285  * Initialization stuff here
00286  */
00287 static molfile_plugin_t plugin;
00288 
00289 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00290   memset(&plugin, 0, sizeof(molfile_plugin_t));
00291   plugin.abiversion = vmdplugin_ABIVERSION;
00292   plugin.type = MOLFILE_PLUGIN_TYPE;
00293   plugin.name = "brix";
00294   plugin.prettyname = "BRIX Density Map";
00295   plugin.author = "Eamon Caddigan";
00296   plugin.majorv = 0;
00297   plugin.minorv = 8;
00298   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00299   plugin.filename_extension = "brix,brx";
00300   plugin.open_file_read = open_brix_read;
00301   plugin.read_volumetric_metadata = read_brix_metadata;
00302   plugin.read_volumetric_data = read_brix_data;
00303   plugin.close_file_read = close_brix_read;
00304 
00305   return VMDPLUGIN_SUCCESS; 
00306 }
00307 
00308 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00309   (*cb)(v, (vmdplugin_t *)&plugin);
00310   return VMDPLUGIN_SUCCESS;
00311 }
00312 
00313 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00314 

Generated on Thu Nov 7 03:09:33 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002