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

ccp4plugin.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: ccp4plugin.C,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.44 $       $Date: 2020/07/06 05:37:19 $
00015  *
00016  ***************************************************************************/
00017 
00018 // EMDB map format specs:
00019 //   ftp://ftp.wwpdb.org/pub/emdb/doc/Map-format/current/EMDB_map_format.pdf
00020 //
00021 // IMOD variants of MRC/CCP4 map format:
00022 //   http://bio3d.colorado.edu/imod/doc/mrc_format.txt
00023 //
00024 // CCP4 electron density map file format description:
00025 //   http://www2.mrc-lmb.cam.ac.uk/image2000.html
00026 //   http://www.ccp4.ac.uk/html/maplib.html
00027 //   http://iims.ebi.ac.uk/3dem-mrc-maps/distribution/mrc_maps.txt
00028 //
00029 // Recent MRC/CCP4 format extension proposals for cryo-EM/ET
00030 //   http://www.sciencedirect.com/science/article/pii/S104784771500074X
00031 //
00032 // TODO: Fix translation/scaling problems found when using 
00033 //       non-orthogonal unit cells.
00034 //
00035 #include <stdio.h>
00036 #include <stdlib.h>
00037 #include <math.h>
00038 #include <string.h>
00039 
00040 #ifndef M_PI
00041 #define M_PI 3.14159265358979323846
00042 #endif
00043 
00044 #include "molfile_plugin.h"
00045 #include "endianswap.h"
00046 
00047 #define CCP4HDSIZE 1024
00048 
00049 // CCP4 and IMOD MRC image types
00050 // based on notes in:
00051 //   http://bio3d.colorado.edu/imod/doc/mrc_format.txt
00052 #define MRC_TYPE_BYTE   0
00053 #define MRC_TYPE_SHORT  1
00054 #define MRC_TYPE_FLOAT  2
00055 #define MRC_TYPE_SHORT2 3
00056 #define MRC_TYPE_FLOAT2 4
00057 #define MRC_TYPE_USHORT 6  /* non-standard */
00058 #define MRC_TYPE_UCHAR3 16 /* non-standard */
00059 
00060 // Generic data flags (plugin-internal)
00061 #define DATA_FLAG_SIGNED               0x1 
00062 #define DATA_FLAG_IMOD_FORMAT          0x2
00063 
00064 // IMOD flags
00065 #define IMOD_FILE_ID            0x494d4f44
00066 #define IMOD_MAGIC_STAMP        1146047817
00067 #define IMOD_FLAG_SIGNED               0x1
00068 #define IMOD_FLAG_HEADER_SPACING       0x2
00069 #define IMOD_FLAG_ORIGIN_INVERTED_SIGN 0x4
00070 
00071 typedef struct {
00072   FILE *fd;
00073   int voxtype; 
00074   long voxcount;
00075   int dataflags;
00076   int imodstamp;
00077   int imodflags;
00078   int nversion;    // MRC/IMOD version conformance tag
00079   char exttype[4]; // extended header type
00080   int nsets;
00081   int swap;
00082   int xyz2crs[3];
00083   float amin; 
00084   float amax; 
00085   float amean;
00086   long dataOffset;
00087   molfile_volumetric_t *vol;
00088 } ccp4_t;
00089 
00090 typedef struct {
00091   unsigned char red;
00092   unsigned char blue;
00093   unsigned char green;
00094 } uchar3;
00095 
00096 static void *open_ccp4_read(const char *filepath, const char *filetype,
00097                             int *natoms) {
00098   FILE *fd;
00099   ccp4_t *ccp4;
00100   char mapString[4], symData[81];
00101   int nxyzstart[3], extent[3], grid[3], crs2xyz[3], voxtype, symBytes;
00102   float origin2k[3];
00103   int swap, i, xIndex, yIndex, zIndex;
00104   long dataOffset, filesize;
00105   float cellDimensions[3], cellAngles[3], xaxis[3], yaxis[3], zaxis[3];
00106   float alpha, beta, gamma, xScale, yScale, zScale, z1, z2, z3;
00107   float amin, amax, amean;
00108   int dataflags=0,imodstamp=0,imodflags=0,nversion=0;
00109   char exttype[4] = {' ', ' ', ' ', ' '};
00110  
00111   fd = fopen(filepath, "rb");
00112   if (!fd) {
00113     printf("ccp4plugin) Error opening file %s\n", filepath);
00114     return NULL;
00115   }
00116 
00117   if ((fread(extent, sizeof(int), 3, fd) != 3) ||
00118       (fread(&voxtype, sizeof(int), 1, fd) != 1) ||
00119       (fread(nxyzstart, sizeof(int), 3, fd) != 3) ||
00120       (fread(grid, sizeof(int), 3, fd) != 3) ||
00121       (fread(cellDimensions, sizeof(float), 3, fd) != 3) ||
00122       (fread(cellAngles, sizeof(float), 3, fd) != 3) ||
00123       (fread(crs2xyz, sizeof(int), 3, fd) != 3) ||
00124       (fread(&amin, sizeof(float), 1, fd) != 1) || 
00125       (fread(&amax, sizeof(float), 1, fd) != 1) || 
00126       (fread(&amean, sizeof(float), 1, fd) != 1)) {
00127     printf("ccp4plugin) Error: Improperly formatted line.\n");
00128     return NULL;
00129   }
00130 
00131   // Check the number of bytes used for storing symmetry operators
00132   // (word 23, byte 92)
00133   fseek(fd, 23 * 4, SEEK_SET);
00134   if ((fread(&symBytes, sizeof(int), 1, fd) != 1) ) {
00135     printf("ccp4plugin) Error: Failed reading symmetry bytes record.\n");
00136     return NULL;
00137   }
00138 
00139   // read MRC2000 Origin record at word 49, byte 196, and use it if necessary
00140   // http://www2.mrc-lmb.cam.ac.uk/image2000.html
00141   fseek(fd, 49 * 4, SEEK_SET);
00142   if (fread(origin2k, sizeof(float), 3, fd) != 3) {
00143     printf("ccp4plugin) Error: unable to read ORIGIN records at offset 196.\n");
00144   }
00145 
00146   // Check for IMOD stamp at offset 152, indicating an IMOD file 
00147   fseek(fd, 152, SEEK_SET);
00148   if (fread(&imodstamp, sizeof(int), 1, fd) != 1) {
00149     printf("ccp4plugin) Error: failed to read IMOD stamp from MRC file.\n");
00150     return NULL;
00151   }
00152   if (fread(&imodflags, sizeof(int), 1, fd) != 1) {
00153     printf("ccp4plugin) Error: failed to read IMOD flags from MRC file.\n");
00154     return NULL;
00155   }
00156 
00157   fseek(fd, 150, SEEK_SET);
00158   if (fread(&exttype, 4*sizeof(char), 1, fd) != 1) {
00159     printf("ccp4plugin) Error: failed to read exttype from MRC file.\n");
00160     return NULL;
00161   }
00162 
00163   // Check for CCP4/IMOD version tag at offset 154
00164   fseek(fd, 154, SEEK_SET);
00165   if (fread(&nversion, sizeof(int), 1, fd) != 1) {
00166     printf("ccp4plugin) Error: failed to read CCP4/IMOD VERSION from MRC file.\n");
00167     return NULL;
00168   }
00169 
00170   // Check file endianism using some heuristics
00171   swap = 0;
00172   int tmp[3];
00173   memcpy(tmp, extent, 3*sizeof(int));
00174   if (tmp[0] > 65536 || tmp[1] > 65536 || tmp[2] > 65536) {
00175     swap4_aligned(tmp, 3);
00176     if (tmp[0] > 65536 || tmp[1] > 65536 || tmp[2] > 65536) {
00177       swap = 0;
00178       printf("ccp4plugin) Guessing file endianism: native\n");
00179     } else {
00180       swap = 1;
00181       printf("ccp4plugin) Guessing file endianism: swapped\n");
00182     }
00183   }
00184   if (voxtype > 16 && swap == 0) {
00185     tmp[0] = voxtype;
00186     swap4_aligned(tmp, 1);
00187     if (tmp[0] <= 16) {
00188       swap = 1;
00189       printf("ccp4plugin) Guessing file endianism: swapped\n");
00190     }
00191   }
00192 
00193   // Byte-swap header information if needed
00194   if (swap == 1) {
00195     swap4_aligned(extent, 3);
00196     swap4_aligned(&voxtype, 1);
00197     swap4_aligned(nxyzstart, 3);
00198     swap4_aligned(origin2k, 3);
00199     swap4_aligned(grid, 3);
00200     swap4_aligned(cellDimensions, 3);
00201     swap4_aligned(cellAngles, 3);
00202     swap4_aligned(crs2xyz, 3);
00203     swap4_aligned(&amin, 1);
00204     swap4_aligned(&amax, 1);
00205     swap4_aligned(&amean, 1);
00206     swap4_aligned(&symBytes, 1);
00207     swap4_aligned(&imodstamp, 1);
00208     swap4_aligned(&imodflags, 1);
00209     swap4_aligned(&nversion, 1);
00210   }
00211 
00212   // Check for the string "MAP" at word 52 byte 208, indicating a CCP4 file.
00213   fseek(fd, 52 * 4, SEEK_SET);
00214   if (fgets(mapString, 4, fd) == NULL) {
00215     printf("ccp4plugin) Error: unable to read 'MAP' string, not a valid CCP4/IMOD MRC file.\n");
00216     return NULL;
00217   }
00218 
00219   if ((strcmp(mapString, "MAP") != 0) && (imodstamp != IMOD_MAGIC_STAMP)) {
00220     // Older versions of IMOD (2.6.19 or below) do not have the "MAP " string.
00221     // If the IMOD stamp is there its probably a valid MRC file
00222     printf("ccp4plugin) Warning: 'MAP' string missing!\n");
00223     printf("ccp4plugin) This usually indicates an invalid IMOD file,\n");
00224     printf("ccp4plugin) but some old IMOD versions did not specify 'MAP'.\n");
00225     printf("ccp4plugin) File loading will proceed but may ultimately fail.\n");
00226   }
00227 
00228   // Check if we found an IMOD file or not
00229   if (imodstamp == IMOD_MAGIC_STAMP) {
00230     printf("ccp4plugin) MRC file generated by IMOD-compatible program.\n");
00231 
00232     dataflags |= DATA_FLAG_IMOD_FORMAT; // mark that this file came from IMOD
00233 
00234     if (imodflags & IMOD_FLAG_SIGNED) {
00235       dataflags |= DATA_FLAG_SIGNED;
00236       printf("ccp4plugin)  IMOD flag: data uses signed-bytes\n");
00237     } else {
00238       printf("ccp4plugin)  IMOD flag: data uses unsigned-bytes\n");
00239     }
00240 
00241     if (imodflags & IMOD_FLAG_HEADER_SPACING)
00242       printf("ccp4plugin)  IMOD flag: pixel spacing set in extended header\n");
00243 
00244     if (imodflags & IMOD_FLAG_ORIGIN_INVERTED_SIGN)
00245       printf("ccp4plugin)  IMOD flag: origin sign is inverted.\n");
00246   } else {
00247     imodflags = 0;
00248     printf("ccp4plugin) No IMOD stamp found.\n");
00249 
00250     if (amin < 0) {
00251       printf("ccp4plugin) Note: amin < 0, signed voxels are implied\n");
00252       dataflags |= DATA_FLAG_SIGNED;
00253     } else {
00254       printf("ccp4plugin) Note: amin >= 0, unsigned voxels are implied\n");
00255     }
00256   }
00257 
00258   // Check the data type of the file.
00259   switch (voxtype) {
00260     case MRC_TYPE_BYTE:
00261       printf("ccp4plugin) voxel type: byte\n");
00262       break;
00263 
00264     case MRC_TYPE_SHORT:
00265       printf("ccp4plugin) voxel type: short (16-bit signed int)\n");
00266       break;
00267 
00268     case MRC_TYPE_FLOAT:
00269       printf("ccp4plugin) voxel type: float (32-bit real)\n");
00270       break;
00271 
00272     case MRC_TYPE_SHORT2:
00273       printf("ccp4plugin) voxel type: short2 (2x 16-bit signed int)\n");
00274       printf("ccp4plugin) Error: unimplemented voxel format\n");
00275       return NULL;
00276 
00277     case MRC_TYPE_FLOAT2:
00278       printf("ccp4plugin) voxel type: float2 (2x 32-bit real)\n");
00279       printf("ccp4plugin) Error: unimplemented voxel format\n");
00280       return NULL;
00281 
00282     case MRC_TYPE_USHORT:
00283       printf("ccp4plugin) voxel type: ushort (16-bit unsigned int)\n");
00284       break;
00285 
00286     case MRC_TYPE_UCHAR3:
00287       printf("ccp4plugin) voxel type: uchar3 (3x unsigned char)\n");
00288       break;
00289 
00290     default:
00291       printf("ccp4plugin) Error: Only byte, short (16-bit integer) or float (32-bit real) data types are supported.\n");
00292       return NULL;
00293   }
00294 
00295 #if 1
00296   printf("ccp4plugin)   nversion: %d\n", nversion);
00297   printf("ccp4plugin)     extent: %d x %d x %d\n",
00298          extent[0], extent[1], extent[2]);
00299   printf("ccp4plugin)  nxyzstart: %d, %d, %d\n", 
00300          nxyzstart[0], nxyzstart[1], nxyzstart[2]);
00301   printf("ccp4plugin)   origin2k: %f, %f, %f\n", 
00302          origin2k[0], origin2k[1], origin2k[2]);
00303   printf("ccp4plugin)       grid: %d x %d x %d\n", grid[0], grid[1], grid[2]);
00304   printf("ccp4plugin)    celldim: %f x %f x %f\n", 
00305          cellDimensions[0], cellDimensions[1], cellDimensions[2]);
00306   printf("cpp4plugin) cellangles: %f, %f, %f\n", 
00307          cellAngles[0], cellAngles[1], cellAngles[2]);
00308   printf("ccp4plugin)    crs2xyz: %d %d %d\n", 
00309          crs2xyz[0], crs2xyz[1], crs2xyz[2]);
00310   printf("ccp4plugin)       amin: %g  amax: %g  amean: %g\n", 
00311          amin, amax, amean);
00312   printf("ccp4plugin)   symBytes: %d\n", symBytes);
00313   printf("ccp4plugin)    extType: %c%c%c%c\n", 
00314          exttype[0], exttype[1], exttype[2], exttype[3]);
00315 #endif
00316 
00317   // Check the dataOffset: this fixes the problem caused by files claiming
00318   // to have symmetry records when they do not.
00319   fseek(fd, 0, SEEK_END);
00320   filesize = ftell(fd);
00321 
00322   // compute data offset using file size and voxel type info
00323   long voxcount = long(extent[0]) * long(extent[1]) * long(extent[2]);
00324   if (voxtype == MRC_TYPE_BYTE) {
00325     dataOffset = filesize - sizeof(char)*voxcount;
00326   } else if (voxtype == MRC_TYPE_FLOAT) {
00327     dataOffset = filesize - sizeof(float)*voxcount;
00328   } else if (voxtype == MRC_TYPE_SHORT || voxtype == MRC_TYPE_USHORT) {
00329     dataOffset = filesize - sizeof(short)*voxcount;
00330   } else if (voxtype == MRC_TYPE_UCHAR3) {
00331     dataOffset = filesize - sizeof(uchar3)*voxcount;
00332   } else {
00333     printf("ccp4plugin) unimplemented voxel type!\n");
00334   }
00335 
00336   if (dataOffset != (CCP4HDSIZE + symBytes)) {
00337     if (dataOffset == CCP4HDSIZE) {
00338       // Bogus symmetry record information
00339       printf("ccp4plugin) Warning: file contains bogus symmetry record.\n");
00340       symBytes = 0;
00341     } else if (dataOffset < CCP4HDSIZE) {
00342       printf("ccp4plugin) Error: File appears truncated and doesn't match header.\n");
00343       return NULL;
00344     } else if ((dataOffset > CCP4HDSIZE) && (dataOffset < (1024*1024))) {
00345       // Fix for loading SPIDER files which are larger than usual
00346       // In this specific case, we must absolutely trust the symBytes record
00347       dataOffset = CCP4HDSIZE + symBytes; 
00348       printf("ccp4plugin) Warning: File is larger than expected and doesn't match header.\n");
00349       printf("ccp4plugin) Warning: Continuing file load, good luck!\n");
00350     } else {
00351       printf("ccp4plugin) Error: File is MUCH larger than expected and doesn't match header.\n");
00352       return NULL;
00353     }
00354   }
00355 
00356   // Read symmetry records -- organized as 80-byte lines of text.
00357   if (symBytes != 0) {
00358     if (dataflags & DATA_FLAG_IMOD_FORMAT) {
00359       printf("ccp4plugin) IMOD format, not printing symmetry string section\n");
00360     } else {
00361       printf("ccp4plugin) Symmetry records found:\n");
00362       fseek(fd, CCP4HDSIZE, SEEK_SET);
00363       for (i = 0; i < symBytes/80; i++) {
00364         fgets(symData, 81, fd);
00365         printf("ccp4plugin) %s\n", symData);
00366       }
00367     }
00368   }
00369 
00370   // check extent and grid interval counts
00371   if (grid[0] == 0 && extent[0] > 0) {
00372     grid[0] = extent[0] - 1;
00373     printf("ccp4plugin) Warning: Fixed X interval count\n");
00374   }
00375   if (grid[1] == 0 && extent[1] > 0) {
00376     grid[1] = extent[1] - 1;
00377     printf("ccp4plugin) Warning: Fixed Y interval count\n");
00378   }
00379   if (grid[2] == 0 && extent[2] > 0) {
00380     grid[2] = extent[2] - 1;
00381     printf("ccp4plugin) Warning: Fixed Z interval count\n");
00382   }
00383 
00384   // Allocate and initialize the ccp4 structure
00385   ccp4 = new ccp4_t;
00386   memset(ccp4, 0, sizeof(ccp4_t));
00387   ccp4->fd = fd;
00388   ccp4->vol = NULL;
00389   *natoms = MOLFILE_NUMATOMS_NONE;
00390   ccp4->nsets = 1; // this EDM file contains only one data set
00391   ccp4->voxtype = voxtype;
00392   ccp4->voxcount = voxcount;
00393   ccp4->dataflags = dataflags;
00394   ccp4->imodstamp = imodstamp;
00395   ccp4->imodflags = imodflags;
00396   ccp4->nversion = nversion;
00397   memcpy(ccp4->exttype, exttype, sizeof(exttype));
00398   ccp4->swap = swap;
00399   ccp4->dataOffset = dataOffset;
00400   ccp4->amin = amin;
00401   ccp4->amax = amax;
00402   ccp4->amean = amean;
00403 
00404   ccp4->vol = new molfile_volumetric_t[1];
00405   memset(ccp4->vol, 0, sizeof(molfile_volumetric_t));
00406   strcpy(ccp4->vol[0].dataname, "CCP4 Electron Density Map");
00407 
00408   // Mapping between CCP4 column, row, section and VMD x, y, z.
00409   if (crs2xyz[0] == 0 && crs2xyz[1] == 0 && crs2xyz[2] == 0) {
00410     printf("ccp4plugin) Warning: All crs2xyz records are zero.\n");
00411     printf("ccp4plugin) Warning: Setting crs2xyz to 1, 2, 3\n");
00412     crs2xyz[0] = 1;
00413     crs2xyz[0] = 2;
00414     crs2xyz[0] = 3;
00415   }
00416 
00417   ccp4->xyz2crs[crs2xyz[0]-1] = 0;
00418   ccp4->xyz2crs[crs2xyz[1]-1] = 1;
00419   ccp4->xyz2crs[crs2xyz[2]-1] = 2;
00420   xIndex = ccp4->xyz2crs[0];
00421   yIndex = ccp4->xyz2crs[1];
00422   zIndex = ccp4->xyz2crs[2];
00423 
00424   // calculate non-orthogonal unit cell coordinates
00425   alpha = (M_PI / 180.0) * cellAngles[0];
00426   beta = (M_PI / 180.0) * cellAngles[1];
00427   gamma = (M_PI / 180.0) * cellAngles[2];
00428 
00429   if (cellDimensions[0] == 0.0 && 
00430       cellDimensions[1] == 0.0 &&
00431       cellDimensions[2] == 0.0) {
00432     printf("ccp4plugin) Warning: Cell dimensions are all zero.\n");
00433     printf("ccp4plugin) Warning: Setting to 1.0, 1.0, 1.0 for viewing.\n");
00434     printf("ccp4plugin) Warning: Map file will not align with other structures.\n");
00435     cellDimensions[0] = 1.0;
00436     cellDimensions[1] = 1.0;
00437     cellDimensions[2] = 1.0;
00438   } 
00439 
00440   xScale = cellDimensions[0] / grid[0];
00441   yScale = cellDimensions[1] / grid[1];
00442   zScale = cellDimensions[2] / grid[2];
00443 
00444   if ((alpha == 0) || (beta == 0) || (gamma == 0)) {
00445     printf("ccp4plugin) orthorhombic unit cell axes\n");
00446     xaxis[0] = xScale;
00447     xaxis[1] = 0;
00448     xaxis[2] = 0;
00449 
00450     yaxis[0] = 0;
00451     yaxis[1] = yScale;
00452     yaxis[2] = 0;
00453 
00454     zaxis[0] = 0;
00455     zaxis[1] = 0;
00456     zaxis[2] = zScale;
00457   } else {
00458     // calculate non-orthogonal unit cell coordinates
00459     printf("ccp4plugin) computing non-orthorhombic unit cell axes\n");
00460     xaxis[0] = xScale;
00461     xaxis[1] = 0;
00462     xaxis[2] = 0;
00463 
00464     yaxis[0] = cos(gamma) * yScale;
00465     yaxis[1] = sin(gamma) * yScale;
00466     yaxis[2] = 0;
00467 
00468     z1 = cos(beta);
00469     z2 = (cos(alpha) - cos(beta)*cos(gamma)) / sin(gamma);
00470     z3 = sqrt(1.0 - z1*z1 - z2*z2);
00471     zaxis[0] = z1 * zScale;
00472     zaxis[1] = z2 * zScale;
00473     zaxis[2] = z3 * zScale;
00474   }
00475 
00476 
00477 #if 0
00478   printf("ccp4plugin)  scale: %g %g %g\n", xScale, yScale, zScale);
00479   printf("ccp4plugin)    xax: %g %g %g\n", xaxis[0], xaxis[1], xaxis[2]);
00480   printf("ccp4plugin)    yax: %g %g %g\n", yaxis[0], yaxis[1], yaxis[2]);
00481   printf("ccp4plugin)    zax: %g %g %g\n", zaxis[0], zaxis[1], zaxis[2]);
00482 #endif
00483 
00484 #if 1
00485   // Handle both MRC-2000 and older format maps
00486   if (origin2k[0] == 0.0f && origin2k[1] == 0.0f && origin2k[2] == 0.0f) {
00487     printf("ccp4plugin) using CCP4 n[xyz]start origin\n");
00488     ccp4->vol[0].origin[0] = xaxis[0] * nxyzstart[xIndex] + 
00489                              yaxis[0] * nxyzstart[yIndex] +
00490                              zaxis[0] * nxyzstart[zIndex];
00491     ccp4->vol[0].origin[1] = yaxis[1] * nxyzstart[yIndex] +
00492                              zaxis[1] * nxyzstart[zIndex];
00493     ccp4->vol[0].origin[2] = zaxis[2] * nxyzstart[zIndex];
00494   } else {
00495     // Use ORIGIN records rather than old n[xyz]start records
00496     //   http://www2.mrc-lmb.cam.ac.uk/image2000.html
00497     // XXX the ORIGIN field is only used by the EM community, and
00498     //     has undefined meaning for non-orthogonal maps and/or
00499     //     non-cubic voxels, etc.
00500     printf("ccp4plugin) using MRC2000 origin\n");
00501     ccp4->vol[0].origin[0] = origin2k[xIndex];
00502     ccp4->vol[0].origin[1] = origin2k[yIndex];
00503     ccp4->vol[0].origin[2] = origin2k[zIndex];
00504   }
00505 #else
00506   // old code that only pays attention to old MRC nxstart/nystart/nzstart
00507   ccp4->vol[0].origin[0] = xaxis[0] * nxyzstart[xIndex] + 
00508                            yaxis[0] * nxyzstart[yIndex] +
00509                            zaxis[0] * nxyzstart[zIndex];
00510   ccp4->vol[0].origin[1] = yaxis[1] * nxyzstart[yIndex] +
00511                            zaxis[1] * nxyzstart[zIndex];
00512   ccp4->vol[0].origin[2] = zaxis[2] * nxyzstart[zIndex];
00513 #endif
00514 
00515   ccp4->vol[0].xaxis[0] = xaxis[0] * (extent[xIndex]-1);
00516   ccp4->vol[0].xaxis[1] = 0;
00517   ccp4->vol[0].xaxis[2] = 0;
00518 
00519   ccp4->vol[0].yaxis[0] = yaxis[0] * (extent[yIndex]-1);
00520   ccp4->vol[0].yaxis[1] = yaxis[1] * (extent[yIndex]-1);
00521   ccp4->vol[0].yaxis[2] = 0;
00522 
00523   ccp4->vol[0].zaxis[0] = zaxis[0] * (extent[zIndex]-1);
00524   ccp4->vol[0].zaxis[1] = zaxis[1] * (extent[zIndex]-1);
00525   ccp4->vol[0].zaxis[2] = zaxis[2] * (extent[zIndex]-1);
00526 
00527   ccp4->vol[0].xsize = extent[xIndex];
00528   ccp4->vol[0].ysize = extent[yIndex];
00529   ccp4->vol[0].zsize = extent[zIndex];
00530 
00531   ccp4->vol[0].has_color = 0;
00532 
00533 #if 0
00534   printf("ccp4plugin) origin: %.3f %.3f %.3f\n",
00535          ccp4->vol[0].origin[0],
00536          ccp4->vol[0].origin[1],
00537          ccp4->vol[0].origin[2]);
00538 
00539   printf("ccp4plugin)   xax': %g %g %g\n",
00540          ccp4->vol[0].xaxis[0], ccp4->vol[0].xaxis[1], ccp4->vol[0].xaxis[2]);
00541   printf("ccp4plugin)   yax': %g %g %g\n",
00542          ccp4->vol[0].yaxis[0], ccp4->vol[0].yaxis[1], ccp4->vol[0].yaxis[2]);
00543   printf("ccp4plugin)   zax': %g %g %g\n",
00544          ccp4->vol[0].zaxis[0], ccp4->vol[0].zaxis[1], ccp4->vol[0].zaxis[2]);
00545 
00546   printf("ccp4plugin)     sz: %d %d %d\n",
00547          ccp4->vol[0].xsize, ccp4->vol[0].ysize, ccp4->vol[0].zsize);
00548 #endif
00549 
00550 
00551   return ccp4;
00552 }
00553 
00554 static int read_ccp4_metadata(void *v, int *nsets, 
00555   molfile_volumetric_t **metadata) {
00556   ccp4_t *ccp4 = (ccp4_t *)v;
00557   *nsets = ccp4->nsets; 
00558   *metadata = ccp4->vol;  
00559 
00560   return MOLFILE_SUCCESS;
00561 }
00562 
00563 static int read_ccp4_data(void *v, int set, float *datablock,
00564                          float *colorblock) {
00565   ccp4_t *ccp4 = (ccp4_t *)v;
00566   int x, y, z, xSize, ySize, zSize, extent[3], coord[3];
00567   long xySize;
00568   FILE *fd = ccp4->fd;
00569 
00570   xSize = ccp4->vol[0].xsize;
00571   ySize = ccp4->vol[0].ysize;
00572   zSize = ccp4->vol[0].zsize;
00573   xySize = xSize * ySize;
00574 
00575   // coord = <col, row, sec>
00576   // extent = <colSize, rowSize, secSize>
00577   extent[ccp4->xyz2crs[0]] = xSize;
00578   extent[ccp4->xyz2crs[1]] = ySize;
00579   extent[ccp4->xyz2crs[2]] = zSize;
00580 
00581   fseek(fd, ccp4->dataOffset, SEEK_SET);
00582 
00583   // Read entire rows of data from the file, then write into the
00584   // datablock with the correct slice ordering.
00585   if ((ccp4->voxtype == MRC_TYPE_BYTE) && (ccp4->dataflags & DATA_FLAG_SIGNED)) {
00586     printf("ccp4plugin) reading signed-byte voxel data\n");
00587     signed char *rowdata = new signed char[extent[0]];
00588     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00589       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00590         if (feof(fd)) {
00591           printf("ccp4plugin) Unexpected end-of-file.\n");
00592           delete [] rowdata;
00593           return MOLFILE_ERROR;
00594         }
00595         if (ferror(fd)) {
00596           printf("ccp4plugin) Problem reading the file.\n");
00597           delete [] rowdata;
00598           return MOLFILE_ERROR;
00599         }
00600         if ( fread(rowdata, sizeof(char), extent[0], fd) != extent[0] ) {
00601           printf("ccp4plugin) Error reading data row.\n");
00602           delete [] rowdata;
00603           return MOLFILE_ERROR;
00604         }
00605 
00606         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00607           x = coord[ccp4->xyz2crs[0]];
00608           y = coord[ccp4->xyz2crs[1]];
00609           z = coord[ccp4->xyz2crs[2]];
00610           datablock[x + long(y*xSize) + long(z*xySize)] = rowdata[coord[0]];
00611         }
00612       }
00613     }
00614     delete [] rowdata;
00615   } else if ((ccp4->voxtype == MRC_TYPE_BYTE) && !(ccp4->dataflags & DATA_FLAG_SIGNED)) {
00616     printf("ccp4plugin) reading unsigned-byte voxel data\n");
00617     unsigned char *rowdata = new unsigned char[extent[0]];
00618     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00619       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00620         if (feof(fd)) {
00621           printf("ccp4plugin) Unexpected end-of-file.\n");
00622           delete [] rowdata;
00623           return MOLFILE_ERROR;
00624         }
00625         if (ferror(fd)) {
00626           printf("ccp4plugin) Problem reading the file.\n");
00627           delete [] rowdata;
00628           return MOLFILE_ERROR;
00629         }
00630         if ( fread(rowdata, sizeof(unsigned char), extent[0], fd) != extent[0] ) {
00631           printf("ccp4plugin) Error reading data row.\n");
00632           delete [] rowdata;
00633           return MOLFILE_ERROR;
00634         }
00635 
00636         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00637           x = coord[ccp4->xyz2crs[0]];
00638           y = coord[ccp4->xyz2crs[1]];
00639           z = coord[ccp4->xyz2crs[2]];
00640           datablock[x + long(y*xSize) + long(z*xySize)] = rowdata[coord[0]];
00641         }
00642       }
00643     }
00644     delete [] rowdata;
00645   } else if (ccp4->voxtype == MRC_TYPE_FLOAT) {
00646     printf("ccp4plugin) reading float (32-bit real) voxel data\n");
00647     float *rowdata = new float[extent[0]];
00648     int x, y, z;
00649     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00650       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00651         if (feof(fd)) {
00652           printf("ccp4plugin) Unexpected end-of-file.\n");
00653           delete [] rowdata;
00654           return MOLFILE_ERROR;
00655         }
00656         if (ferror(fd)) {
00657           printf("ccp4plugin) Problem reading the file.\n");
00658           delete [] rowdata;
00659           return MOLFILE_ERROR;
00660         }
00661         if (fread(rowdata, sizeof(float), extent[0], fd) != extent[0] ) {
00662           printf("ccp4plugin) Error reading data row.\n");
00663           delete [] rowdata;
00664           return MOLFILE_ERROR;
00665         }
00666 
00667         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00668           x = coord[ccp4->xyz2crs[0]];
00669           y = coord[ccp4->xyz2crs[1]];
00670           z = coord[ccp4->xyz2crs[2]];
00671           datablock[x + long(y*xSize) + long(z*xySize)] = rowdata[coord[0]];
00672         }
00673       }
00674     }
00675     delete [] rowdata;
00676     if (ccp4->swap == 1)
00677       swap4_aligned(datablock, xySize * zSize);
00678   } else if (ccp4->voxtype == MRC_TYPE_SHORT) {
00679     printf("ccp4plugin) reading short (16-bit int) voxel data\n");
00680     short *rowdata = new short[extent[0]];
00681     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00682       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00683         if (feof(fd)) {
00684           printf("ccp4plugin) Unexpected end-of-file.\n");
00685           delete [] rowdata;
00686           return MOLFILE_ERROR;
00687         }
00688         if (ferror(fd)) {
00689           printf("ccp4plugin) Problem reading the file.\n");
00690           delete [] rowdata;
00691           return MOLFILE_ERROR;
00692         }
00693         if (fread(rowdata, sizeof(short), extent[0], fd) != extent[0] ) {
00694           printf("ccp4plugin) Error reading data row.\n");
00695           delete [] rowdata;
00696           return MOLFILE_ERROR;
00697         }
00698         if (ccp4->swap == 1)
00699           swap2_aligned(rowdata, extent[0]);
00700         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00701           x = coord[ccp4->xyz2crs[0]];
00702           y = coord[ccp4->xyz2crs[1]];
00703           z = coord[ccp4->xyz2crs[2]];
00704           datablock[x + long(y*xSize) + long(z*xySize)] = rowdata[coord[0]];
00705         }
00706       }
00707     }
00708     delete [] rowdata;
00709   } else if (ccp4->voxtype == MRC_TYPE_SHORT2) {
00710     /* IMOD developers said that this is not used anymore and not worth our time to implement */
00711     printf("TYPE_SHORT2 not implemented yet...\n");
00712     return MOLFILE_ERROR;
00713   } else if (ccp4->voxtype == MRC_TYPE_USHORT) {
00714     printf("ccp4plugin) reading unsigned short (16-bit int) voxel data\n");
00715     unsigned short *rowdata = new unsigned short[extent[0]];
00716     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00717       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00718         if (feof(fd)) {
00719           printf("ccp4plugin) Unexpected end-of-file.\n");
00720           delete [] rowdata;
00721           return MOLFILE_ERROR;
00722         }
00723         if (ferror(fd)) {
00724           printf("ccp4plugin) Problem reading the file.\n");
00725           delete [] rowdata;
00726           return MOLFILE_ERROR;
00727         }
00728         if (fread(rowdata, sizeof(unsigned short), extent[0], fd) != extent[0] ) {
00729           printf("ccp4plugin) Error reading data row.\n");
00730           delete [] rowdata;
00731           return MOLFILE_ERROR;
00732         }
00733         if (ccp4->swap == 1)
00734           swap2_aligned(rowdata, extent[0]);
00735         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00736           x = coord[ccp4->xyz2crs[0]];
00737           y = coord[ccp4->xyz2crs[1]];
00738           z = coord[ccp4->xyz2crs[2]];
00739           datablock[x + long(y*xSize) + long(z*xySize)] = rowdata[coord[0]];
00740         }
00741       }
00742     }
00743     delete [] rowdata;
00744   } else if (ccp4->voxtype == MRC_TYPE_UCHAR3) {
00745     printf("ccp4plugin) reading unsigned char * 3 (8-bit uchar * 3) voxel data\n");
00746     uchar3 *rowdata = new uchar3[extent[0]];
00747     float grayscale;
00748     for (coord[2] = 0; coord[2] < extent[2]; coord[2]++) {
00749       for (coord[1] = 0; coord[1] < extent[1]; coord[1]++) {
00750         if (feof(fd)) {
00751           printf("ccp4plugin) Unexpected end-of-file.\n");
00752           delete [] rowdata;
00753           return MOLFILE_ERROR;
00754         }
00755         if (ferror(fd)) {
00756           printf("ccp4plugin) Problem reading the file.\n");
00757           delete [] rowdata;
00758           return MOLFILE_ERROR;
00759         }
00760         if ( fread(rowdata, sizeof(uchar3), extent[0], fd) != extent[0] ) {
00761           printf("ccp4plugin) Error reading data row.\n");
00762           delete [] rowdata;
00763           return MOLFILE_ERROR;
00764         }
00765         for (coord[0] = 0; coord[0] < extent[0]; coord[0]++) {
00766           x = coord[ccp4->xyz2crs[0]];
00767           y = coord[ccp4->xyz2crs[1]];
00768           z = coord[ccp4->xyz2crs[2]];
00769           grayscale = rowdata[coord[0]].red + rowdata[coord[0]].blue + rowdata[coord[0]].green;
00770           datablock[x + long(y*xSize) + long(z*xySize)] = grayscale/3.0;
00771         }
00772       }
00773     }
00774     delete [] rowdata;
00775   }
00776   return MOLFILE_SUCCESS;
00777 }
00778 
00779 static void close_ccp4_read(void *v) {
00780   ccp4_t *ccp4 = (ccp4_t *)v;
00781 
00782   fclose(ccp4->fd);
00783   delete [] ccp4->vol; 
00784   delete ccp4;
00785 }
00786 
00787 /*
00788  * Initialization stuff here
00789  */
00790 static molfile_plugin_t plugin;
00791 
00792 VMDPLUGIN_API int VMDPLUGIN_init(void) { 
00793   memset(&plugin, 0, sizeof(molfile_plugin_t));
00794   plugin.abiversion = vmdplugin_ABIVERSION;
00795   plugin.type = MOLFILE_PLUGIN_TYPE;
00796   plugin.name = "ccp4";
00797   plugin.prettyname = "CCP4, MRC Density Map";
00798   plugin.author = "Eamon Caddigan, Brendan McMorrow, John Stone";
00799   plugin.majorv = 2;
00800   plugin.minorv = 2;
00801   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
00802   plugin.filename_extension = "ccp4,mrc,map";
00803   plugin.open_file_read = open_ccp4_read;
00804   plugin.read_volumetric_metadata = read_ccp4_metadata;
00805   plugin.read_volumetric_data = read_ccp4_data;
00806   plugin.close_file_read = close_ccp4_read;
00807   return VMDPLUGIN_SUCCESS; 
00808 }
00809 
00810 
00811 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00812   (*cb)(v, (vmdplugin_t *)&plugin);
00813   return VMDPLUGIN_SUCCESS;
00814 }
00815 
00816 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }

Generated on Tue Aug 11 03:06:28 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002