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

dsn6plugin.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: dsn6plugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.25 $       $Date: 2016/11/28 05:01:53 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* 
00019  * DSN6 format electron density maps.
00020  *
00021  * More info for format can be found at 
00022  * <http://www.uoxray.uoregon.edu/tnt/manual/node104.html>
00023  * TODO: Check byte-swapping and alignment issues in read_data()
00024  *
00025  * apparently there are some gotchas there, and mapman does some things to 
00026  * fix up some variants of DSN6:
00027  * Gerard "DVD" Kleywegt  gerard@xray.bmc.uu.se: 
00028  *   "dale's description is largely correct, but in elements 13 to 15
00029  *    the actual cell angles are stored (multiplied by a scale factor)
00030  *    rather than their cosines. also, turbo-frodo-style dsn6 maps
00031  *    differ slightly from o-style dsn6 maps. shown below is a chunk
00032  *    of code (yes, fortran) from mapman that fills the header record"
00033  * See original email for mapman code snippet:
00034  *      http://o-info.bioxray.dk/pipermail/o-info/2002-June/005993.html
00035  * 
00036  *
00037  */
00038 
00039 #include <stdlib.h>
00040 #include <stdio.h>
00041 #include <ctype.h>
00042 #include <math.h>
00043 #include <string.h>
00044 
00045 #if defined(_AIX)
00046 #include <strings.h>
00047 #endif
00048 
00049 #if defined(WIN32) || defined(WIN64)
00050 #define strcasecmp stricmp
00051 #endif
00052 
00053 #ifndef M_PI
00054 #define M_PI 3.14159265358979323846
00055 #endif
00056 
00057 #include "molfile_plugin.h"
00058 #include "endianswap.h"
00059 
00060 typedef struct {
00061   FILE *fd;
00062   int nsets;
00063   float prod, plus;
00064   molfile_volumetric_t *vol;
00065 } dsn6_t;
00066 
00067 
00068 static void *open_dsn6_read(const char *filepath, const char *filetype,
00069     int *natoms) {
00070   FILE *fd;
00071   dsn6_t *dsn6;
00072   short fileHeader[19];
00073   int start_x, start_y, start_z, extent_x, extent_y, extent_z;
00074   float scale, A, B, C, alpha, beta, gamma, 
00075         xaxis[3], yaxis[3], zaxis[3], z1, z2, z3;
00076   
00077   fd = fopen(filepath, "rb");
00078   if (!fd) {
00079     fprintf(stderr, "Error opening file.\n");
00080     return NULL;
00081   }
00082 
00083   // Read the header into a 19-element int array. The integers are stored
00084   // according to the endianness of the machine used to write the file, so
00085   // swapping may be necessary.
00086   fread(fileHeader, sizeof(short), 19, fd);
00087 
00088   // Check byte-order, swapping if necessary
00089   if (fileHeader[18] == 25600)
00090     swap2_aligned(fileHeader, 19);
00091   else if (fileHeader[18] != 100) {
00092     fprintf(stderr, "Error reading file header.\n");
00093     return NULL;
00094   }
00095   // else fileHeader[18] is 100, byte-order is fine 
00096   // (this value is hard-coded into the file format)
00097 
00098   // Unit cell origin, in grid coordinates
00099   start_x = fileHeader[0];
00100   start_y = fileHeader[1];
00101   start_z = fileHeader[2];
00102 
00103   // Unit cell size
00104   extent_x = fileHeader[3];
00105   extent_y = fileHeader[4];
00106   extent_z = fileHeader[5];
00107 
00108   // Constant cell scaling factor
00109   scale = 1.0 / fileHeader[17];
00110 
00111   // Unit cell edge, in angstroms / sample
00112   A = scale * fileHeader[9] / fileHeader[6];
00113   B = scale * fileHeader[10] / fileHeader[7];
00114   C = scale * fileHeader[11] / fileHeader[8];
00115 
00116   // Cell angles, in radians
00117   alpha = scale * fileHeader[12] * M_PI / 180.0;
00118   beta = scale * fileHeader[13] * M_PI / 180.0;
00119   gamma = scale * fileHeader[14] * M_PI / 180.0;
00120 
00121   // Allocate and initialize the dsn6 structure
00122   dsn6 = new dsn6_t;
00123   dsn6->fd = fd;
00124   dsn6->vol = NULL;
00125   *natoms = MOLFILE_NUMATOMS_NONE;
00126   dsn6->nsets = 1; // this file contains only one data set
00127 
00128   // Grid data scaling constants
00129   dsn6->prod = (float) fileHeader[15] / fileHeader[18];
00130   dsn6->plus = fileHeader[16];
00131 
00132   dsn6->vol = new molfile_volumetric_t[1];
00133   strcpy(dsn6->vol[0].dataname, "DSN6 Electron Density Map");
00134 
00135   // Calculate non-orthogonal unit cell coordinates
00136   xaxis[0] = A;
00137   xaxis[1] = 0;
00138   xaxis[2] = 0;
00139 
00140   yaxis[0] = cos(gamma) * B;
00141   yaxis[1] = sin(gamma) * B;
00142   yaxis[2] = 0;
00143 
00144   z1 = cos(beta);
00145   z2 = (cos(alpha) - cos(beta)*cos(gamma)) / sin(gamma);
00146   z3 = sqrt(1.0 - z1*z1 - z2*z2);
00147   zaxis[0] = z1 * C;
00148   zaxis[1] = z2 * C;
00149   zaxis[2] = z3 * C;
00150 
00151   // Convert the origin from grid space to cartesian coordinates
00152   dsn6->vol[0].origin[0] = xaxis[0] * start_x + yaxis[0] * start_y + 
00153                            zaxis[0] * start_z;
00154   dsn6->vol[0].origin[1] = yaxis[1] * start_y + zaxis[1] * start_z;
00155   dsn6->vol[0].origin[2] = zaxis[2] * start_z;
00156 
00157   dsn6->vol[0].xaxis[0] = xaxis[0] * (extent_x-1);
00158   dsn6->vol[0].xaxis[1] = 0;
00159   dsn6->vol[0].xaxis[2] = 0;
00160 
00161   dsn6->vol[0].yaxis[0] = yaxis[0] * (extent_y-1);
00162   dsn6->vol[0].yaxis[1] = yaxis[1] * (extent_y-1);
00163   dsn6->vol[0].yaxis[2] = 0;
00164 
00165   dsn6->vol[0].zaxis[0] = zaxis[0] * (extent_z-1);
00166   dsn6->vol[0].zaxis[1] = zaxis[1] * (extent_z-1);
00167   dsn6->vol[0].zaxis[2] = zaxis[2] * (extent_z-1);
00168 
00169   dsn6->vol[0].xsize = extent_x;
00170   dsn6->vol[0].ysize = extent_y;
00171   dsn6->vol[0].zsize = extent_z;
00172 
00173   dsn6->vol[0].has_color = 0;
00174 
00175   return dsn6;
00176 }
00177 
00178 static int read_dsn6_metadata(void *v, int *nsets, 
00179   molfile_volumetric_t **metadata) {
00180   dsn6_t *dsn6 = (dsn6_t *)v;
00181   *nsets = dsn6->nsets; 
00182   *metadata = dsn6->vol;  
00183 
00184   return MOLFILE_SUCCESS;
00185 }
00186 
00187 static int read_dsn6_data(void *v, int set, float *datablock,
00188                          float *colorblock) {
00189   dsn6_t *dsn6 = (dsn6_t *)v;
00190   float * cell = datablock;
00191   unsigned char brick[512];
00192   unsigned char* brickPtr = NULL;
00193   int xsize, ysize, zsize, xysize, xbrix, ybrix, zbrix, cellIndex;
00194   int x, y, z, xbrik, ybrik, zbrik;
00195   FILE * fd = dsn6->fd;
00196   float div, plus;
00197 
00198   // Read 512-byte "bricks" of data. Each brick contains data for 8*8*8 
00199   // gridpoints.
00200   fseek(fd, 512, SEEK_SET);
00201 
00202   div = 1.0 / dsn6->prod;
00203   plus = dsn6->plus;
00204 
00205   xsize = dsn6->vol[0].xsize;
00206   ysize = dsn6->vol[0].ysize;
00207   zsize = dsn6->vol[0].zsize;
00208   xysize = xsize * ysize;
00209 
00210   xbrix = (int) ceil((float) xsize / 8.0);
00211   ybrix = (int) ceil((float) ysize / 8.0);
00212   zbrix = (int) ceil((float) zsize / 8.0);
00213 
00214   cellIndex = 0;
00215   for (zbrik = 0; zbrik < zbrix; zbrik++) {
00216     for (ybrik = 0; ybrik < ybrix; ybrik++) {
00217       for (xbrik = 0; xbrik < xbrix; xbrik++) {
00218         // Read the next brick into the buffer and swap its bytes.
00219         if (feof(fd)) {
00220           fprintf(stderr, "Unexpected end-of-file.\n");
00221           return MOLFILE_ERROR;
00222         }
00223         if (ferror(fd)) {
00224           fprintf(stderr, "Error reading file.\n");
00225           return MOLFILE_ERROR;
00226         }
00227 
00228         fread(brick, sizeof(char), 512, fd);
00229         swap2_unaligned(brick, 512*sizeof(char));
00230         brickPtr = brick;
00231 
00232         for (z = 0; z < 8; z++) {
00233           if ((z + zbrik*8) >= zsize) {
00234             cellIndex += (8 - z) * xysize;
00235             break;
00236           }
00237 
00238           for (y = 0; y < 8; y++) {
00239             if ((y + ybrik*8) >= ysize) {
00240               cellIndex += (8 - y) * xsize;
00241               brickPtr += (8 - y) * 8;
00242               break;
00243             }
00244 
00245             for (x = 0; x < 8; x++) {
00246               if ((x + xbrik*8) >= xsize) {
00247                 cellIndex += 8 - x;
00248                 brickPtr += 8 - x;
00249                 break;
00250               }
00251 
00252               // cell[(x+xbrik*8) + (y+ybrik*8)*xsize + (z+zbrik*8)*xysize] =
00253               // div * ((float) brick[x + 8*y + 64*z] - plus)
00254               cell[cellIndex] = div * ((float) *brickPtr - plus);
00255 
00256               brickPtr++;
00257               cellIndex++;
00258             } // end for(x)
00259            
00260             cellIndex += xsize - 8;
00261           } // end for(y)
00262          
00263           cellIndex += xysize - 8*xsize;
00264         } // end for(z)
00265       
00266         cellIndex += 8 - 8*xysize; 
00267       } // end for(xbrik)
00268 
00269       cellIndex += 8 * (xsize - xbrix);
00270     } // end for(ybrik)
00271 
00272     cellIndex += 8 * (xysize - xsize*ybrik);
00273   } // end for(zbrik)
00274  
00275   return MOLFILE_SUCCESS;
00276 }
00277 
00278 static void close_dsn6_read(void *v) {
00279   dsn6_t *dsn6 = (dsn6_t *)v;
00280 
00281   fclose(dsn6->fd);
00282   if (dsn6->vol != NULL)
00283     delete [] dsn6->vol; 
00284   delete dsn6;
00285 }
00286 
00287 /*
00288  * Initialization stuff here
00289  */
00290 static molfile_plugin_t plugin;
00291 
00292 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00293   memset(&plugin, 0, sizeof(molfile_plugin_t));
00294   plugin.abiversion = vmdplugin_ABIVERSION;
00295   plugin.type = MOLFILE_PLUGIN_TYPE;
00296   plugin.name = "DSN6";
00297   plugin.prettyname = "DSN6";
00298   plugin.author = "Eamon Caddigan";
00299   plugin.majorv = 0;
00300   plugin.minorv = 6;
00301   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00302   plugin.filename_extension = "ds6,dsn6,omap";
00303   plugin.open_file_read = open_dsn6_read;
00304   plugin.read_volumetric_metadata = read_dsn6_metadata;
00305   plugin.read_volumetric_data = read_dsn6_data;
00306   plugin.close_file_read = close_dsn6_read;
00307   return VMDPLUGIN_SUCCESS; 
00308 }
00309 
00310 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00311   (*cb)(v, (vmdplugin_t *)&plugin);
00312   return VMDPLUGIN_SUCCESS;
00313 }
00314 
00315 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
00316 

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