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

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