00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
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
00051
00052
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
00059 #define MRC_TYPE_UCHAR3 16
00060
00061
00062 #define DATA_FLAG_SIGNED 0x1
00063 #define DATA_FLAG_IMOD_FORMAT 0x2
00064
00065
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;
00080 char exttype[4];
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
00133
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
00141
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
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
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
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
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
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
00222
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
00230 if (imodstamp == IMOD_MAGIC_STAMP) {
00231 printf("ccp4plugin) MRC file generated by IMOD-compatible program.\n");
00232
00233 dataflags |= DATA_FLAG_IMOD_FORMAT;
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
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
00319
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
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
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
00358
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
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
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
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;
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
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
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
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
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
00508
00509
00510
00511
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
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
00588
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
00596
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
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
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; }