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

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