25 void read_vol(
char *vol_file,
double *width,
double *origx,
double *origy,
26 double *origz,
unsigned *extx,
unsigned *exty,
unsigned *extz,
29 int n_range_viol0, n_range_viol1, n_range_viol2;
30 int ordermode = 7, cubic = 1, orom = 1;
33 int nxstart = 0, nystart = 0, nzstart = 0;
34 int ncstart = 0, nrstart = 0, nsstart = 0;
35 double widthx, widthy, widthz;
36 double xorigin, yorigin, zorigin;
40 n_range_viol0 =
test_mrc(vol_file, 0);
41 n_range_viol1 =
test_mrc(vol_file, 1);
44 if (n_range_viol2 < n_range_viol0 && n_range_viol2 < n_range_viol1)
45 read_situs(vol_file, width, origx, origy, origz, extx, exty, extz, phi);
47 read_mrc(vol_file, &orom, &cubic, &ordermode, &nc, &nr, &ns, &ncstart,
48 &nrstart, &nsstart, &widthx, &widthy, &widthz, &xorigin, &yorigin,
49 &zorigin, &alpha, &beta, &gamma, &phi_raw);
50 permute_map(ordermode, nc, nr, ns, &nx, &ny, &nz, ncstart, nrstart, nsstart,
51 &nxstart, &nystart, &nzstart, phi_raw, phi);
53 nsstart, nxstart, nystart, nzstart);
55 nx, ny, nz, nxstart, nystart, nzstart, xorigin, yorigin,
56 zorigin, width, origx, origy, origz, extx, exty, extz, phi);
57 printf(
"lib_vio> \n");
65 void write_vol(
char *vol_file,
double width,
double origx,
double origy,
66 double origz,
unsigned extx,
unsigned exty,
unsigned extz,
71 write_situs(vol_file, width, origx, origy, origz, extx, exty, extz, phi);
73 write_mrc(1, vol_file, width, origx, origy, origz, extx, exty, extz, phi);
80 void read_situs(
char *vol_file,
double *width,
double *origx,
double *origy,
81 double *origz,
unsigned *extx,
unsigned *exty,
unsigned *extz,
84 unsigned long nvox, count;
85 double dorigx, dorigy, dorigz, dwidth;
86 double phitest, dtemp;
87 const char *program =
"lib_vio";
90 fin = fopen(vol_file,
"r");
96 if (7 != fscanf(fin,
"%le %le %le %le %d %d %d", &dwidth, &dorigx, &dorigy, &dorigz, extx, exty, extz))
error_fscanf(program, vol_file);
101 printf(
"lib_vio> Situs formatted map file %s - Header information: \n", vol_file);
102 printf(
"lib_vio> Columns, rows, and sections: x=1-%d, y=1-%d, z=1-%d\n", *extx, *exty, *extz);
103 printf(
"lib_vio> 3D coordinates of first voxel: (%f,%f,%f)\n", *origx, *origy, *origz);
104 printf(
"lib_vio> Voxel size in Angstrom: %f \n", *width);
106 nvox = *extx * *exty * *extz;
109 printf(
"lib_vio> Reading density data... \n");
110 *phi = (
double *)
alloc_vect(nvox,
sizeof(
double));
112 for (count = 0; count < nvox; count++) {
113 if (fscanf(fin,
"%le", &dtemp) != 1) {
115 }
else *(*phi + count) = dtemp;
117 if (fscanf(fin,
"%le", &phitest) != EOF) {
121 printf(
"lib_vio> Volumetric data read from file %s\n", vol_file);
127 void write_situs(
char *vol_file,
double width,
double origx,
double origy,
128 double origz,
unsigned extx,
unsigned exty,
unsigned extz,
131 unsigned long nvox, count;
132 const char *program =
"lib_vio";
135 nvox = extx * exty * extz;
136 fout = fopen(vol_file,
"w");
141 printf(
"lib_vio> Writing density data... \n");
142 fprintf(fout,
"%f %f %f %f %d %d %d\n", width, origx, origy, origz, extx, exty, extz);
145 for (count = 0; count < nvox; count++) {
146 if ((count + 1) % 10 == 0)
fprintf(fout,
" %10.6f \n", *(phi + count));
147 else fprintf(fout,
" %10.6f ", *(phi + count));
150 printf(
"lib_vio> Volumetric data written in Situs format to file %s \n", vol_file);
153 printf(
"lib_vio> Situs formatted map file %s - Header information: \n", vol_file);
154 printf(
"lib_vio> Columns, rows, and sections: x=1-%d, y=1-%d, z=1-%d\n", extx, exty, extz);
155 printf(
"lib_vio> 3D coordinates of first voxel: (%f,%f,%f)\n", origx, origy, origz);
156 printf(
"lib_vio> Voxel size in Angstrom: %f \n", width);
163 void read_ascii(
char *vol_file,
unsigned long nvox,
double **fphi)
169 fin = fopen(vol_file,
"r");
173 printf(
"lib_vio> Reading ASCII data... \n");
175 for (count = 0; count < nvox; count++) {
176 if (fscanf(fin,
"%e", &currfloat) != 1) {
179 *(*fphi + count) = currfloat;
182 if (fscanf(fin,
"%e", &currfloat) != EOF) {
185 printf(
"lib_vio> Volumetric data read from file %s\n", vol_file);
192 double *widthx,
double *widthy,
double *widthz,
193 double *alpha,
double *
beta,
double *
gamma,
194 int *nxstart,
int *nystart,
int *nzstart,
195 unsigned *extx,
unsigned *exty,
unsigned *extz,
double **fphi)
202 float a_tmp, b_tmp, g_tmp, currfloat;
203 long mx, my, mz, mxend, myend, mzend;
204 long mxstart, mystart, mzstart;
205 int testa, testb, testg;
206 float xlen, ylen, zlen;
207 int indx, indy, indz;
211 fin = fopen(vol_file,
"r");
221 if (sscanf(nextline,
"%8ld%8ld%8ld%8ld%8ld%8ld%8ld%8ld%8ld", &mx, &mxstart, &mxend, &my, &mystart, &myend, &mz, &mzstart, &mzend) != 9) {
224 *extx = mxend - mxstart + 1;
225 *exty = myend - mystart + 1;
226 *extz = mzend - mzstart + 1;
227 nvox = *extx * *exty * *extz;
229 printf(
"lib_vio> X-PLOR map indexing (counting from 0): \n");
230 printf(
"lib_vio> NA = %8ld (# of X intervals in unit cell) \n", mx);
231 printf(
"lib_vio> AMIN = %8ld (start index X) \n", mxstart);
232 printf(
"lib_vio> AMAX = %8ld (end index X) \n", mxend);
233 printf(
"lib_vio> NB = %8ld (# of Y intervals in unit cell) \n", my);
234 printf(
"lib_vio> BMIN = %8ld (start index Y) \n", mystart);
235 printf(
"lib_vio> BMAX = %8ld (end index Y) \n", myend);
236 printf(
"lib_vio> NC = %8ld (# of Z intervals in unit cell) \n", mz);
237 printf(
"lib_vio> CMIN = %8ld (start index Z) \n", mzstart);
238 printf(
"lib_vio> CMAX = %8ld (end index Z) \n", mzend);
242 if (sscanf(nextline,
"%12f%12f%12f%12f%12f%12f", &xlen, &ylen, &zlen, &a_tmp, &b_tmp, &g_tmp) != 6) {
250 printf(
"lib_vio> X-PLOR unit cell info: \n");
251 printf(
"lib_vio> A = %8.3f (unit cell dimension) \n", xlen);
252 printf(
"lib_vio> B = %8.3f (unit cell dimension) \n", ylen);
253 printf(
"lib_vio> C = %8.3f (unit cell dimension) \n", zlen);
254 printf(
"lib_vio> ALPHA = %8.3f (unit cell angle) \n", *alpha);
255 printf(
"lib_vio> BETA = %8.3f (unit cell angle) \n", *beta);
256 printf(
"lib_vio> GAMMA = %8.3f (unit cell angle) \n", *gamma);
259 *widthx = xlen / (double) mx;
260 *widthy = ylen / (double) my;
261 *widthz = zlen / (double) mz;
267 testa = (int)
floor(100 * *alpha + 0.5);
268 testb = (int)
floor(100 * *beta + 0.5);
269 testg = (int)
floor(100 * *gamma + 0.5);
270 if (testa != 9000 || testb != 9000 || testg != 9000) *orom = 0;
271 if (*orom == 0 ||
floor((*widthx - *widthy) * 1000 + 0.5) != 0 ||
272 floor((*widthy - *widthz) * 1000 + 0.5) != 0 ||
273 floor((*widthx - *widthz) * 1000 + 0.5) != 0)
277 for (done = 0; done == 0;) {
278 if (fgets(nextline,
FLENGTH, fin) == NULL) {
281 if (*nextline ==
'Z' && *(nextline + 1) ==
'Y' && *(nextline + 2) ==
'X') {
285 if (*nextline ==
'z' && *(nextline + 1) ==
'y' && *(nextline + 2) ==
'x') {
293 for (indz = 0; indz < (int)*extz; indz++) {
296 if (sscanf(nextline,
"%8d", &idummy) != 1) {
299 if (idummy != indz) {
304 for (indy = 0; indy < (int)*exty; indy++)
for (indx = 0; indx < (int)*extx; indx++) {
305 if (fscanf(fin,
"%12f", &currfloat) != 1) {
308 *(*fphi +
gidz_general(indz, indy, indx, *exty, *extx)) = currfloat;
315 if (sscanf(nextline,
"%8d", &idummy) != 1 || idummy != -9999) {
321 printf(
"lib_vio> Volumetric data read from file %s\n", vol_file);
329 int i,
done, foundletter;
331 for (done = 0; done == 0;) {
332 if (fgets(*nextline,
FLENGTH, *fin) == NULL) {
336 if ((*(*nextline + i) ==
'!') || (*(*nextline + i) ==
'\n')) {
337 *(*nextline +
i) =
'\0';
341 for (i = 0; * (*nextline +
i) !=
'\0'; ++
i) {
342 if (*(*nextline + i) > 64 && *(*nextline +
i) < 68) {
346 if (*(*nextline + i) > 69 && *(*nextline +
i) < 91) {
350 if (*(*nextline + i) > 96 && *(*nextline +
i) < 100) {
354 if (*(*nextline + i) > 101 && *(*nextline +
i) < 123) {
359 if (foundletter == 0)
for (i = 0; * (*nextline +
i) !=
'\0'; ++
i)
360 if (*(*nextline + i) >=
'0' && *(*nextline +
i) <=
'9') {
369 void read_mrc(
char *vol_file,
int *orom,
int *cubic,
int *ordermode,
370 unsigned *nc,
unsigned *nr,
unsigned *ns,
371 int *ncstart,
int *nrstart,
int *nsstart,
372 double *widthx,
double *widthy,
double *widthz,
373 double *xorigin,
double *yorigin,
double *zorigin,
374 double *alpha,
double *
beta,
double *
gamma,
378 unsigned long count, nvox;
380 int nc_tmp, nr_tmp, ns_tmp, mx, my, mz;
382 float a_tmp, b_tmp, g_tmp;
383 float x_tmp, y_tmp, z_tmp;
384 int testa, testb, testg;
385 float xlen, ylen, zlen;
386 int i, swap, header_ok = 1;
387 int mapc, mapr, maps;
388 float dmin,
dmax, dmean,
dummy, currfloat, cfunsigned, prevfloat = 0.0f, pfunsigned = 0.0f,
rms;
389 double totdiffsigned = 0.0, totdiffunsigned = 0.0;
390 int n_range_viol0, n_range_viol1;
391 char mapchar1, mapchar2, mapchar3, mapchar4;
392 char machstchar1, machstchar2, machstchar3, machstchar4;
393 int ispg, nsymbt, lskflg;
394 float skwmat11, skwmat12, skwmat13, skwmat21, skwmat22, skwmat23, skwmat31, skwmat32, skwmat33, skwtrn1, skwtrn2, skwtrn3;
396 n_range_viol0 =
test_mrc(vol_file, 0);
397 n_range_viol1 =
test_mrc(vol_file, 1);
399 if (n_range_viol0 < n_range_viol1) {
401 if (n_range_viol0 > 0) {
402 printf(
"lib_vio> Warning: %i header field range violations detected in file %s \n", n_range_viol0, vol_file);
406 if (n_range_viol1 > 0) {
407 printf(
"lib_vio> Warning: %i header field range violations detected in file %s \n", n_range_viol1, vol_file);
412 fin = fopen(vol_file,
"rb");
416 printf(
"lib_vio> Reading header information from MRC or CCP4 file %s \n", vol_file);
417 header_ok *=
read_int(&nc_tmp, fin, swap);
418 header_ok *=
read_int(&nr_tmp, fin, swap);
419 header_ok *=
read_int(&ns_tmp, fin, swap);
423 header_ok *=
read_int(&mode, fin, swap);
424 header_ok *=
read_int(ncstart, fin, swap);
425 header_ok *=
read_int(nrstart, fin, swap);
426 header_ok *=
read_int(nsstart, fin, swap);
427 header_ok *=
read_int(&mx, fin, swap);
428 header_ok *=
read_int(&my, fin, swap);
429 header_ok *=
read_int(&mz, fin, swap);
439 header_ok *=
read_int(&mapc, fin, swap);
440 header_ok *=
read_int(&mapr, fin, swap);
441 header_ok *=
read_int(&maps, fin, swap);
445 header_ok *=
read_int(&ispg, fin, swap);
446 header_ok *=
read_int(&nsymbt, fin, swap);
447 header_ok *=
read_int(&lskflg, fin, swap);
448 header_ok *=
read_float(&skwmat11, fin, swap);
449 header_ok *=
read_float(&skwmat12, fin, swap);
450 header_ok *=
read_float(&skwmat13, fin, swap);
451 header_ok *=
read_float(&skwmat21, fin, swap);
452 header_ok *=
read_float(&skwmat22, fin, swap);
453 header_ok *=
read_float(&skwmat23, fin, swap);
454 header_ok *=
read_float(&skwmat31, fin, swap);
455 header_ok *=
read_float(&skwmat32, fin, swap);
456 header_ok *=
read_float(&skwmat33, fin, swap);
460 for (i = 38; i < 50; ++
i) header_ok *=
read_float(&dummy, fin, swap);
471 header_ok *=
read_char(&machstchar1, fin);
472 header_ok *=
read_char(&machstchar2, fin);
473 header_ok *=
read_char(&machstchar3, fin);
474 header_ok *=
read_char(&machstchar4, fin);
477 if (header_ok == 0) {
482 printf(
"lib_vio> NC = %8d (# columns)\n", *nc);
483 printf(
"lib_vio> NR = %8d (# rows)\n", *nr);
484 printf(
"lib_vio> NS = %8d (# sections)\n", *ns);
485 printf(
"lib_vio> MODE = %8d (data type: 0: 8-bit, 1: 16-bit; 2: 32-bit)\n", mode);
486 printf(
"lib_vio> NCSTART = %8d (index of first column, counting from 0)\n", *ncstart);
487 printf(
"lib_vio> NRSTART = %8d (index of first row, counting from 0)\n", *nrstart);
488 printf(
"lib_vio> NSSTART = %8d (index of first section, counting from 0)\n", *nsstart);
489 printf(
"lib_vio> MX = %8d (# of X intervals in unit cell)\n", mx);
490 printf(
"lib_vio> MY = %8d (# of Y intervals in unit cell)\n", my);
491 printf(
"lib_vio> MZ = %8d (# of Z intervals in unit cell)\n", mz);
492 printf(
"lib_vio> X length = %8.3f (unit cell dimension)\n", xlen);
493 printf(
"lib_vio> Y length = %8.3f (unit cell dimension)\n", ylen);
494 printf(
"lib_vio> Z length = %8.3f (unit cell dimension)\n", zlen);
495 printf(
"lib_vio> Alpha = %8.3f (unit cell angle)\n", *alpha);
496 printf(
"lib_vio> Beta = %8.3f (unit cell angle)\n", *beta);
497 printf(
"lib_vio> Gamma = %8.3f (unit cell angle)\n", *gamma);
498 printf(
"lib_vio> MAPC = %8d (columns axis: 1=X,2=Y,3=Z)\n", mapc);
499 printf(
"lib_vio> MAPR = %8d (rows axis: 1=X,2=Y,3=Z)\n", mapr);
500 printf(
"lib_vio> MAPS = %8d (sections axis: 1=X,2=Y,3=Z)\n", maps);
501 printf(
"lib_vio> DMIN = %8.3f (minimum density value - ignored)\n", dmin);
502 printf(
"lib_vio> DMAX = %8.3f (maximum density value - ignored)\n", dmax);
503 printf(
"lib_vio> DMEAN = %8.3f (mean density value - ignored)\n", dmean);
504 printf(
"lib_vio> ISPG = %8d (space group number - ignored)\n", ispg);
505 printf(
"lib_vio> NSYMBT = %8d (# bytes storing symmetry operators)\n", nsymbt);
506 printf(
"lib_vio> LSKFLG = %8d (skew matrix flag: 0:none, 1:follows)\n", lskflg);
508 printf(
"lib_vio> Warning: Will compute skew parameters from the above header info.\n");
509 printf(
"lib_vio> The following skew parameters in the header will be ignored:\n");
510 printf(
"lib_vio> S11 = %8.3f (skew matrix element 11)\n", skwmat11);
511 printf(
"lib_vio> S12 = %8.3f (skew matrix element 12)\n", skwmat12);
512 printf(
"lib_vio> S13 = %8.3f (skew matrix element 13)\n", skwmat13);
513 printf(
"lib_vio> S21 = %8.3f (skew matrix element 21)\n", skwmat21);
514 printf(
"lib_vio> S22 = %8.3f (skew matrix element 22)\n", skwmat22);
515 printf(
"lib_vio> S23 = %8.3f (skew matrix element 23)\n", skwmat23);
516 printf(
"lib_vio> S31 = %8.3f (skew matrix element 31)\n", skwmat31);
517 printf(
"lib_vio> S32 = %8.3f (skew matrix element 32)\n", skwmat32);
518 printf(
"lib_vio> S33 = %8.3f (skew matrix element 33)\n", skwmat33);
519 printf(
"lib_vio> T1 = %8.3f (skew translation element 1)\n", skwtrn1);
520 printf(
"lib_vio> T2 = %8.3f (skew translation element 2)\n", skwtrn2);
521 printf(
"lib_vio> T3 = %8.3f (skew translation element 3)\n", skwtrn3);
522 printf(
"lib_vio> End of skew parameters warning.\n");
524 printf(
"lib_vio> XORIGIN = %8.3f (X origin - MRC2000 only)\n", *xorigin);
525 printf(
"lib_vio> YORIGIN = %8.3f (Y origin - MRC2000 only)\n", *yorigin);
526 printf(
"lib_vio> ZORIGIN = %8.3f (Z origin - MRC2000 only)\n", *zorigin);
527 printf(
"lib_vio> MAP = %c%c%c%c (map string)\n", mapchar1, mapchar2, mapchar3, mapchar4);
528 printf(
"lib_vio> MACHST = %d %d %d %d (machine stamp - ignored)\n", machstchar1, machstchar2, machstchar3, machstchar4);
529 printf(
"lib_vio> RMS = %8.3f (density rms deviation -ignored)\n",
rms);
530 if (mapchar1 !=
'M' || mapchar2 !=
'A' || mapchar3 !=
'P') {
531 printf(
"lib_vio> Warning: MAP string not detected. It appears this is an older (pre-2000) or unconventional MRC format.\n");
532 printf(
"lib_vio> If in doubt use the map2map tool in manual mode and/or inspect the results.\n");
537 nvox = *nc * *nr * *ns;
544 for (count = 0; count < nvox; ++count) {
547 *(*fphi + count) = currfloat;
548 if (currfloat < 0) cfunsigned = currfloat + 256.0;
549 else cfunsigned = currfloat;
550 totdiffsigned += fabs(currfloat - prevfloat);
551 totdiffunsigned += fabs(cfunsigned - pfunsigned);
552 prevfloat = currfloat;
553 pfunsigned = cfunsigned;
557 if (totdiffsigned > totdiffunsigned) {
558 for (count = 0; count < nvox; count++)
if (*(*fphi + count) < 0) *(*fphi + count) += 256.0;
559 printf(
"lib_vio> Warning: It appears the MODE 0 data type uses the unsigned 8-bit convention, adding 256 to negative values.\n");
560 printf(
"lib_vio> If in doubt use the map2map tool in manual mode and/or inspect the results.\n");
569 for (count = 0; count < (
unsigned long)nsymbt; ++count)
if (
read_char_float(&currfloat, fin) == 0) {
573 for (count = 0; count < nvox; ++count) {
577 *(*fphi + count) = currfloat;
588 for (count = 0; count < (
unsigned long)nsymbt; ++count)
if (
read_char_float(&currfloat, fin) == 0) {
591 for (count = 0; count < nvox; ++count) {
595 *(*fphi + count) = currfloat;
605 *widthx = xlen / (double) mx;
606 *widthy = ylen / (double) my;
607 *widthz = zlen / (double) mz;
610 testa = (int)
floor(100 * *alpha + 0.5);
611 testb = (int)
floor(100 * *beta + 0.5);
612 testg = (int)
floor(100 * *gamma + 0.5);
613 if (testa != 9000 || testb != 9000 || testg != 9000) *orom = 0;
614 if (*orom == 0 ||
floor((*widthx - *widthy) * 1000 + 0.5) != 0 ||
floor((*widthy - *widthz) * 1000 + 0.5) != 0 ||
floor((*widthx - *widthz) * 1000 + 0.5) != 0) *cubic = 0;
618 if (mapc == 1 && mapr == 2 && maps == 3) {
621 if (mapc == 1 && mapr == 3 && maps == 2) {
624 if (mapc == 2 && mapr == 1 && maps == 3) {
627 if (mapc == 2 && mapr == 3 && maps == 1) {
630 if (mapc == 3 && mapr == 1 && maps == 2) {
633 if (mapc == 3 && mapr == 2 && maps == 1) {
636 if (*ordermode == 7) {
639 printf(
"lib_vio> Volumetric data read from file %s\n", vol_file);
644 void read_spider(
char *vol_file,
unsigned *extx,
unsigned *exty,
unsigned *extz,
double **fphi)
648 int i, swap, header_ok = 1, headlen;
650 float dummy, currfloat;
651 float nslice, nrow, iform, imami,
fmax, fmin, av, sig, nsam, headrec;
652 float iangle, phi,
theta,
gamma, xoff, yoff, zoff, scale, labbyt, lenbyt;
653 float istack, inuse, maxim, kangle, phi1, theta1, psi1, phi2, theta2, psi2;
654 int n_range_viol0, n_range_viol1;
659 if (n_range_viol0 < n_range_viol1) {
661 if (n_range_viol0 > 0) {
662 printf(
"lib_vio> Warning: %i header field range violations detected \n", n_range_viol0);
666 if (n_range_viol1 > 0) {
667 printf(
"lib_vio> Warning: %i header field range violations detected \n", n_range_viol1);
672 fin = fopen(vol_file,
"rb");
676 printf(
"lib_vio> Reading header information from SPIDER file %s \n", vol_file);
703 for (i = 0; i < 4; ++
i) header_ok *=
read_float(&dummy, fin, swap);
711 if (header_ok == 0) {
716 printf(
"lib_vio> NSLICE = %8.f (# sections)\n", nslice);
717 printf(
"lib_vio> NROW = %8.f (# rows)\n", nrow);
718 printf(
"lib_vio> IFORM = %8.f (file type specifier - ignored)\n", iform);
719 printf(
"lib_vio> IMAMI = %8.f (flag: 1 = the maximum and minimum are computed - ignored)\n", imami);
720 printf(
"lib_vio> FMAX = %8.3f (maximum density value - ignored)\n", fmax);
721 printf(
"lib_vio> FMIN = %8.3f (minimum density value - ignored)\n", fmin);
722 printf(
"lib_vio> AV = %8.3f (average density value - ignored)\n", av);
723 printf(
"lib_vio> SIG = %8.3f (density rms deviation - ignored)\n", sig);
724 printf(
"lib_vio> NSAM = %8.f (# columns)\n", nsam);
725 printf(
"lib_vio> HEADREC = %8.f (number of records in file header)\n", headrec);
726 printf(
"lib_vio> IANGLE = %8.f (flag: 1 = tilt angles filled - ignored)\n", iangle);
727 printf(
"lib_vio> PHI = %8.3f (tilt angle - ignored)\n", phi);
728 printf(
"lib_vio> THETA = %8.3f (tilt angle - ignored)\n", theta);
729 printf(
"lib_vio> GAMMA = %8.3f (tilt angle - ignored)\n", gamma);
730 printf(
"lib_vio> XOFF = %8.3f (X offset - ignored)\n", xoff);
731 printf(
"lib_vio> YOFF = %8.3f (Y offset - ignored)\n", yoff);
732 printf(
"lib_vio> ZOFF = %8.3f (Z offset - ignored)\n", zoff);
733 printf(
"lib_vio> SCALE = %8.3f (scale factor - ignored)\n", scale);
734 printf(
"lib_vio> LABBYT = %8.f (total number of bytes in header)\n", labbyt);
735 printf(
"lib_vio> LENBYT = %8.f (record length in bytes)\n", lenbyt);
736 printf(
"lib_vio> ISTACK = %8.f (flag; file contains a stack of images - ignored)\n", istack);
737 printf(
"lib_vio> INUSE = %8.f (flag; this image in stack is used - ignored)\n", inuse);
738 printf(
"lib_vio> MAXIM = %8.f (maximum image used in stack - ignored)\n", maxim);
739 printf(
"lib_vio> KANGLE = %8.f (flag: additional angles set - ignored)\n", kangle);
740 printf(
"lib_vio> PHI1 = %8.3f (additional rotation - ignored)\n", phi1);
741 printf(
"lib_vio> THETA1 = %8.3f (additional rotation - ignored)\n", theta1);
742 printf(
"lib_vio> PSI1 = %8.3f (additional rotation - ignored)\n", psi1);
743 printf(
"lib_vio> PHI2 = %8.3f (additional rotation - ignored)\n", phi2);
744 printf(
"lib_vio> THETA2 = %8.3f (additional rotation - ignored)\n", theta2);
745 printf(
"lib_vio> PSI2 = %8.3f (additional rotation - ignored)\n", psi2);
747 *extx = (
unsigned int)nsam;
748 *exty = (
unsigned int)nrow;
749 *extz = (
unsigned int)nslice;
750 nvox = *extx * *exty * *extz;
751 headlen = *extx * (int)ceil(256 / (*extx * 1.0));
754 for (count = 0; count < (
unsigned long)headlen; ++count)
if (
read_float_empty(fin) == 0) {
757 for (count = 0; count < nvox; ++count)
if (
read_float(&currfloat, fin, swap) == 0) {
760 *(*fphi + count) = currfloat;
765 if (((
int)
floor(100 * (
sizeof(
float)*headlen - labbyt) + 0.5)) != 0) {
768 printf(
"lib_vio> Volumetric data read from file %s\n", vol_file);
778 unsigned long nfloat;
781 fin = fopen(in_file,
"rb");
788 phi = (
float *)
alloc_vect(nfloat,
sizeof(
float));
790 for (count = 0; count < nfloat; count++)
793 printf(
"lib_vio> Binary data read as 32-bit 'float' type from file %s\n", in_file);
794 if (swap) printf(
"lib_vio> The byte order (endianism) has been swapped.\n");
797 fout = fopen(out_file,
"w");
802 printf(
"lib_vio> Writing ASCII data... \n");
803 for (count = 0; count < nfloat; count++) {
804 if ((count + 1) % 10 == 0)
fprintf(fout,
" %10.6f \n", *(phi + count));
805 else fprintf(fout,
" %10.6f ", *(phi + count));
809 printf(
"lib_vio> ASCII dump of binary file %s written to file %s. \n", in_file, out_file);
810 printf(
"lib_vio> Open / check ASCII file with text editor and extract map densities. \n");
811 printf(
"lib_vio> Then convert with map2map tool, option 1: ASCII.\n");
819 unsigned char *cptr, tmp;
821 if (fread(currfloat, 4, 1, fin) != 1)
return 0;
823 cptr = (
unsigned char *)currfloat;
838 unsigned char *cptr, tmp;
841 if (fread(&currshort, 2, 1, fin) != 1)
return 0;
843 cptr = (
unsigned char *)&currshort;
848 *currfloat = (float)currshort;
858 if (fread(&currfloat, 4, 1, fin) != 1)
return 0;
867 if (fread(&currchar, 1, 1, fin) != 1)
return 0;
868 *currfloat = (float)currchar;
877 if (fread(currchar, 1, 1, fin) != 1)
return 0;
885 unsigned char *cptr, tmp;
887 if (fread(currlong, 4, 1, fin) != 1)
return 0;
889 cptr = (
unsigned char *)currlong;
904 unsigned long finl = 0;
908 if (fgetc(*fin) == EOF)
break;
919 float xreg, xreg1, yreg, yreg1, zreg, zreg1;
921 xreg = fabs(fmod(origx + 0.00001 * width, width));
922 xreg1 = fabs(fmod(origx - 0.00001 * width, width));
923 yreg = fabs(fmod(origy + 0.00001 * width, width));
924 yreg1 = fabs(fmod(origy - 0.00001 * width, width));
925 zreg = fabs(fmod(origz + 0.00001 * width, width));
926 zreg1 = fabs(fmod(origz - 0.00001 * width, width));
927 if (xreg1 < xreg) xreg = xreg1;
928 if (yreg1 < yreg) yreg = yreg1;
929 if (zreg1 < zreg) zreg = zreg1;
930 if (xreg + yreg + zreg > 0.0001 * width)
return 0;
941 double width, origx, origy, origz;
942 unsigned extx, exty, extz;
943 int header_ok = 1, n_range_viols = 0;
945 fin = fopen(vol_file,
"r");
951 if (fscanf(fin,
"%le", &width) != 1) header_ok *= 0;
952 if (fscanf(fin,
"%le", &origx) != 1) header_ok *= 0;
953 if (fscanf(fin,
"%le", &origy) != 1) header_ok *= 0;
954 if (fscanf(fin,
"%le", &origz) != 1) header_ok *= 0;
955 if (fscanf(fin,
"%d", &extx) != 1) header_ok *= 0;
956 if (fscanf(fin,
"%d", &exty) != 1) header_ok *= 0;
957 if (fscanf(fin,
"%d", &extz) != 1) header_ok *= 0;
960 if (header_ok == 0)
return 999999;
962 n_range_viols += (extx > 5000);
963 n_range_viols += (extx < 0);
964 n_range_viols += (exty > 5000);
965 n_range_viols += (exty < 0);
966 n_range_viols += (extz > 5000);
967 n_range_viols += (extz < 0);
968 n_range_viols += (width > 2000);
969 n_range_viols += (width < 0);
970 n_range_viols += (origx > 1e6);
971 n_range_viols += (origx < -1e6);
972 n_range_viols += (origy > 1e6);
973 n_range_viols += (origy < -1e6);
974 n_range_viols += (origz > 1e6);
975 n_range_viols += (origz < -1e6);
976 return 2 * n_range_viols;
985 int situs_format = 0;
987 vl = strlen(vol_file);
989 if (strstr(vol_file + (vl - 4),
".sit") != NULL) situs_format = 1;
990 if (strstr(vol_file + (vl - 4),
".SIT") != NULL) situs_format = 1;
993 if (strstr(vol_file + (vl - 6),
".situs") != NULL) situs_format = 1;
994 if (strstr(vol_file + (vl - 6),
".SITUS") != NULL) situs_format = 1;
1013 int nc, nr, ns, mx, my, mz;
1014 int mode, ncstart, nrstart, nsstart;
1015 float xlen, ylen, zlen;
1016 int i, header_ok = 1, n_range_viols = 0;
1017 int mapc, mapr, maps;
1021 fin = fopen(vol_file,
"rb");
1027 header_ok *=
read_int(&nc, fin, swap);
1028 header_ok *=
read_int(&nr, fin, swap);
1029 header_ok *=
read_int(&ns, fin, swap);
1030 header_ok *=
read_int(&mode, fin, swap);
1031 header_ok *=
read_int(&ncstart, fin, swap);
1032 header_ok *=
read_int(&nrstart, fin, swap);
1033 header_ok *=
read_int(&nsstart, fin, swap);
1034 header_ok *=
read_int(&mx, fin, swap);
1035 header_ok *=
read_int(&my, fin, swap);
1036 header_ok *=
read_int(&mz, fin, swap);
1043 header_ok *=
read_int(&mapc, fin, swap);
1044 header_ok *=
read_int(&mapr, fin, swap);
1045 header_ok *=
read_int(&maps, fin, swap);
1049 for (i = 23; i < 50; ++
i) header_ok *=
read_float(&dummy, fin, swap);
1050 header_ok *=
read_float(&xorigin, fin, swap);
1051 header_ok *=
read_float(&yorigin, fin, swap);
1052 header_ok *=
read_float(&zorigin, fin, swap);
1054 if (header_ok == 0) {
1058 n_range_viols += (nc > 5000);
1059 n_range_viols += (nc < 0);
1060 n_range_viols += (nr > 5000);
1061 n_range_viols += (nr < 0);
1062 n_range_viols += (ns > 5000);
1063 n_range_viols += (ns < 0);
1064 n_range_viols += (ncstart > 5000);
1065 n_range_viols += (ncstart < -5000);
1066 n_range_viols += (nrstart > 5000);
1067 n_range_viols += (nrstart < -5000);
1068 n_range_viols += (nsstart > 5000);
1069 n_range_viols += (nsstart < -5000);
1070 n_range_viols += (mx > 5000);
1071 n_range_viols += (mx < 0);
1072 n_range_viols += (my > 5000);
1073 n_range_viols += (my < 0);
1074 n_range_viols += (mz > 5000);
1075 n_range_viols += (mz < 0);
1076 n_range_viols += (alpha > 360.0f);
1077 n_range_viols += (alpha < -360.0f);
1078 n_range_viols += (beta > 360.0f);
1079 n_range_viols += (beta < -360.0f);
1080 n_range_viols += (gamma > 360.0f);
1081 n_range_viols += (gamma < -360.0f);
1083 return n_range_viols;
1092 int i, header_ok = 1, n_range_viols = 0, headlen;
1094 float nslice, nrow, iform, imami,
fmax, fmin, av, sig, nsam, headrec;
1095 float iangle, phi,
theta,
gamma, xoff, yoff, zoff, scale, labbyt, lenbyt;
1096 float istack, inuse, maxim, kangle, phi1, theta1, psi1, phi2, theta2, psi2;
1098 fin = fopen(vol_file,
"rb");
1116 header_ok *=
read_float(&headrec, fin, swap);
1130 for (i = 0; i < 4; ++
i) header_ok *=
read_float(&dummy, fin, swap);
1139 if (header_ok == 0) {
1142 headlen = (int)(nsam * ceil(256 / (nsam * 1.0)));
1143 n_range_viols += (((int)
floor(100 * (4 * headlen - labbyt) + 0.5)) != 0);
1144 n_range_viols += (headrec > 1e10);
1145 n_range_viols += (headrec <= 1e-10);
1146 n_range_viols += (labbyt > 1e10);
1147 n_range_viols += (labbyt <= 1e-10);
1148 n_range_viols += (lenbyt > 1e10);
1149 n_range_viols += (lenbyt <= 1e-10);
1150 n_range_viols += (nslice > 1e10);
1151 n_range_viols += (nslice < 1e-10);
1152 n_range_viols += (nrow > 1e10);
1153 n_range_viols += (nrow < 1e-10);
1154 n_range_viols += (nsam > 1e10);
1155 n_range_viols += (nsam < 1e-10);
1156 return 2 * n_range_viols;
1166 unsigned nr,
unsigned ns)
1168 unsigned ic,
ir, is;
1169 unsigned long ncr, q;
1173 q = count - is * ncr;
1177 switch (ordermode) {
1179 return ic + ir * nc + is * nc * nr;
1181 return ic + is * nc + ir * nc * ns;
1183 return ir + ic * nr + is * nr * nc;
1185 return is + ic * ns + ir * ns * nc;
1187 return ir + is * nr + ic * nr * ns;
1189 return is + ir * ns + ic * ns * nr;
1198 void permute_map(
int ordermode,
unsigned nc,
unsigned nr,
unsigned ns,
1199 unsigned *nx,
unsigned *ny,
unsigned *nz,
int ncstart,
1200 int nrstart,
int nsstart,
int *nxstart,
int *nystart,
1201 int *nzstart,
double *phi,
double **pphi)
1204 unsigned long nvox, count;
1207 nvox = nc * nr * ns;
1209 for (count = 0; count < nvox; count++)
1210 *(*pphi +
permuted_index(ordermode, count, nc, nr, ns)) = *(phi + count);
1215 nsstart, nxstart, nystart, nzstart);
1221 unsigned *nx,
unsigned *ny,
unsigned *nz,
int ncstart,
1222 int nrstart,
int nsstart,
int *nxstart,
int *nystart,
1226 switch (ordermode) {
1282 void permute_print_info(
int ordermode,
unsigned nc,
unsigned nr,
unsigned ns,
unsigned nx,
unsigned ny,
unsigned nz,
int ncstart,
int nrstart,
int nsstart,
int nxstart,
int nystart,
int nzstart)
1285 if (ordermode > 1) {
1286 printf(
"lib_vio> Map parameters BEFORE selected axis permutation: \n");
1287 printf(
"lib_vio> NC = %8d (# columns)\n", nc);
1288 printf(
"lib_vio> NR = %8d (# rows)\n", nr);
1289 printf(
"lib_vio> NS = %8d (# sections)\n", ns);
1290 printf(
"lib_vio> NCSTART = %8d (index of first column, counting from 0)\n", ncstart);
1291 printf(
"lib_vio> NRSTART = %8d (index of first row, counting from 0)\n", nrstart);
1292 printf(
"lib_vio> NSSTART = %8d (index of first section, counting from 0)\n", nsstart);
1293 printf(
"lib_vio> Map parameters AFTER selected axis permutation: \n");
1294 printf(
"lib_vio> NX = %8d (# X fields)\n", nx);
1295 printf(
"lib_vio> NY = %8d (# Y fields)\n", ny);
1296 printf(
"lib_vio> NZ = %8d (# Z fields)\n", nz);
1297 printf(
"lib_vio> NXSTART = %8d (X index of first voxel, counting from 0)\n", nxstart);
1298 printf(
"lib_vio> NYSTART = %8d (Y index of first voxel, counting from 0)\n", nystart);
1299 printf(
"lib_vio> NZSTART = %8d (Z index of first voxel, counting from 0)\n", nzstart);
1301 printf(
"lib_vio> C,R,S = X,Y,Z (no axis permutation). \n");
1307 int set_origin_get_mode(
int nxstart,
int nystart,
int nzstart,
double widthx,
double widthy,
double widthz,
double xorigin,
double yorigin,
double zorigin,
double *origx,
double *origy,
double *origz)
1310 if (fabs(xorigin) < 0.0001 && fabs(yorigin) < 0.0001 && fabs(zorigin) < 0.0001) {
1311 printf(
"lib_vio> Using crystallographic (CCP4) style origin defined by unit cell start indices.\n");
1312 *origx = nxstart * widthx;
1313 *origy = nystart * widthy;
1314 *origz = nzstart * widthz;
1317 printf(
"lib_vio> Using MRC2000 style origin defined by [X,Y,Z]ORIGIN fields.\n");
1327 double gamma,
double widthx,
double widthy,
1328 double widthz,
unsigned nx,
unsigned ny,
unsigned nz,
1329 int nxstart,
int nystart,
int nzstart,
double xorigin,
1330 double yorigin,
double zorigin,
double *pwidth,
1331 double *porigx,
double *porigy,
double *porigz,
1332 unsigned *pextx,
unsigned *pexty,
unsigned *pextz,
1338 double origx, origy, origz;
1343 xorigin, yorigin, zorigin, porigx, porigy, porigz);
1344 printf(
"lib_vio> Cubic lattice present. \n");
1352 xorigin, yorigin, zorigin, &origx, &origy, &origz);
1353 printf(
"lib_vio> Orthogonal lattice with unequal spacings in x, y, z detected.\n");
1354 printf(
"lib_vio> Interpolating map to cubic lattice using smallest detected spacing rounded to nearest 0.1 Angstrom.\n");
1357 if (widthy < *pwidth) *pwidth = widthy;
1358 if (widthz < *pwidth) *pwidth = widthz;
1359 *pwidth =
floor(*pwidth * 10.0 + 0.5) / 10.0;
1362 *pwidth, *pwidth, *pwidth, *pphi, nx, ny, nz, origx,
1363 origy, origz, widthx, widthy, widthz);
1364 nvox = *pextx * *pexty * *pextz;
1367 printf(
"lib_vio> Skewed, non-orthogonal lattice detected.\n");
1368 printf(
"lib_vio> Interpolating skewed map to cubic lattice using smallest detected spacing rounded to nearest 0.1 Angstrom.\n");
1371 porigz, pwidth, *pphi, nx, ny, nz,
1372 nxstart, nystart, nzstart, xorigin,
1373 yorigin, zorigin, widthx, widthy, widthz,
1374 alpha, beta, gamma);
1375 nvox = *pextx * *pexty * *pextz;
1382 fprintf(stderr,
"lib_vio> Input grid not in register with origin of coordinate system.\n");
1389 unsigned *pexty,
unsigned *pextz,
1390 double *porigx,
double *porigy,
1391 double *porigz,
double *pwidth,
1392 double *phiin,
unsigned nx,
unsigned ny,
1393 unsigned nz,
int nxstart,
int nystart,
1394 int nzstart,
double xorigin,
1395 double yorigin,
double zorigin,
1396 double widthx,
double widthy,
1397 double widthz,
double alpha,
double beta,
1401 unsigned long pnvox;
1402 double ax, ay, az, bx, by, bz, cx, cy, cz, aix, aiy, aiz, bix, biy, biz, cix, ciy, ciz;
1403 double t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, cdet, scz;
1404 double ux[8], uy[8], uz[8], uxmin, uxmax, uymin, uymax, uzmin, uzmax;
1405 double xpos, ypos, zpos, gx, gy, gz,
a,
b,
c;
1406 int x0,
y0,
z0, x1, y1, z1;
1407 double endox, endoy, endoz;
1408 int i, indx, indy, indz;
1409 double origx, origy, origz;
1416 bx = cos(
PI * gamma / 180.0);
1417 by = sin(
PI * gamma / 180.0);
1419 t1x = cos(
PI * (gamma + alpha) / 180.0);
1420 t1y = sin(
PI * (gamma + alpha) / 180.0);
1421 t2x = cos(
PI * (gamma - alpha) / 180.0);
1422 t2y = sin(
PI * (gamma - alpha) / 180.0);
1423 t3x = cos(
PI * beta / 180.0);
1424 t3y = sin(
PI * beta / 180.0);
1425 t4x = cos(
PI * beta / 180.0);
1426 t4y = -1.0 * sin(
PI * beta / 180.0);
1427 cdet = (t4y - t3y) * (t2x - t1x) - (t2y - t1y) * (t4x - t3x);
1428 if (fabs(cdet) < 1E-15) {
1431 cx = ((t4x - t3x) * (t1y * t2x - t2y * t1x) -
1432 (t2x - t1x) * (t3y * t4x - t4y * t3x)) / cdet;
1433 cy = ((t4y - t3y) * (t1y * t2x - t2y * t1x) -
1434 (t2y - t1y) * (t3y * t4x - t4y * t3x)) / cdet;
1435 scz = 1.0 - (cx * cx + cy * cy);
1448 cix = (bx * cy - cx * by) / (by * cz);
1449 ciy = -cy / (by * cz);
1453 if (
set_origin_get_mode(nxstart, nystart, nzstart, widthx, widthy, widthz, xorigin, yorigin, zorigin, &origx, &origy, &origz)) {
1454 gx = aix * origx + bix * origy + cix * origz;
1455 gy = aiy * origx + biy * origy + ciy * origz;
1456 gz = aiz * origx + biz * origy + ciz * origz;
1463 endox = origx + (nx - 1) * widthx;
1464 endoy = origy + (ny - 1) * widthy;
1465 endoz = origz + (nz - 1) * widthz;
1466 ux[0] = ax * origx + bx * origy + cx * origz;
1467 uy[0] = ay * origx + by * origy + cy * origz;
1468 uz[0] = az * origx + bz * origy + cz * origz;
1469 ux[1] = ax * endox + bx * origy + cx * origz;
1470 uy[1] = ay * endox + by * origy + cy * origz;
1471 uz[1] = az * endox + bz * origy + cz * origz;
1472 ux[2] = ax * origx + bx * endoy + cx * origz;
1473 uy[2] = ay * origx + by * endoy + cy * origz;
1474 uz[2] = az * origx + bz * endoy + cz * origz;
1475 ux[3] = ax * origx + bx * origy + cx * endoz;
1476 uy[3] = ay * origx + by * origy + cy * endoz;
1477 uz[3] = az * origx + bz * origy + cz * endoz;
1478 ux[4] = ax * endox + bx * endoy + cx * origz;
1479 uy[4] = ay * endox + by * endoy + cy * origz;
1480 uz[4] = az * endox + bz * endoy + cz * origz;
1481 ux[5] = ax * origx + bx * endoy + cx * endoz;
1482 uy[5] = ay * origx + by * endoy + cy * endoz;
1483 uz[5] = az * origx + bz * endoy + cz * endoz;
1484 ux[6] = ax * endox + bx * origy + cx * endoz;
1485 uy[6] = ay * endox + by * origy + cy * endoz;
1486 uz[6] = az * endox + bz * origy + cz * endoz;
1487 ux[7] = ax * endox + bx * endoy + cx * endoz;
1488 uy[7] = ay * endox + by * endoy + cy * endoz;
1489 uz[7] = az * endox + bz * endoy + cz * endoz;
1491 for (i = 0; i < 8; i++)
if (ux[i] < uxmin) uxmin = ux[
i];
1493 for (i = 0; i < 8; i++)
if (uy[i] < uymin) uymin = uy[
i];
1495 for (i = 0; i < 8; i++)
if (uz[i] < uzmin) uzmin = uz[
i];
1497 for (i = 0; i < 8; i++) if (ux[i] > uxmax) uxmax = ux[i];
1499 for (i = 0; i < 8; i++) if (uy[i] > uymax) uymax = uy[i];
1501 for (i = 0; i < 8; i++) if (uz[i] > uzmax) uzmax = uz[i];
1505 if ((widthy * by) < *pwidth) *pwidth = widthy * by;
1506 if ((widthz * cz) < *pwidth) *pwidth = widthz * cz;
1507 *pwidth =
floor(*pwidth * 10.0 + 0.5) / 10.0;
1511 *porigx = *pwidth *
floor(uxmin / *pwidth);
1512 *porigy = *pwidth *
floor(uymin / *pwidth);
1513 *porigz = *pwidth *
floor(uzmin / *pwidth);
1517 *pextx = (int)(ceil(uxmax / *pwidth) -
floor(uxmin / *pwidth) + 1.5);
1518 *pexty = (int)(ceil(uymax / *pwidth) -
floor(uymin / *pwidth) + 1.5);
1519 *pextz = (int)(ceil(uzmax / *pwidth) -
floor(uzmin / *pwidth) + 1.5);
1520 pnvox = *pextx * *pexty * *pextz;
1524 for (indz = 0; indz < (int)*pextz; indz++)
1525 for (indy = 0; indy < (int)*pexty; indy++)
1526 for (indx = 0; indx < (int)*pextx; indx++) {
1528 xpos = *porigx + indx * *pwidth;
1529 ypos = *porigy + indy * *pwidth;
1530 zpos = *porigz + indz * *pwidth;
1533 gx = (aix * xpos + bix * ypos + cix * zpos - origx) / widthx;
1534 gy = (aiy * xpos + biy * ypos + ciy * zpos - origy) / widthy;
1535 gz = (aiz * xpos + biz * ypos + ciz * zpos - origz) / widthz;
1536 x0 = (int)
floor(gx);
1537 y0 = (int)
floor(gy);
1538 z0 = (int)
floor(gz);
1544 if (x0 >= 0 && x1 < (
int)nx && y0 >= 0 && y1 < (int)ny &&
1545 z0 >= 0 && z1 < (
int)nz) {
1549 *(*pphiout +
gidz_general(indz, indy, indx, *pexty, *pextx)) =
1550 a * b * c * *(phiin +
gidz_general(z1, y1, x1, ny, nx)) +
1551 (1 - a) * b * c * *(phiin +
gidz_general(z1, y1, x0, ny, nx)) +
1552 a * (1 - b) * c * *(phiin +
gidz_general(z1, y0, x1, ny, nx)) +
1553 a * b * (1 - c) * *(phiin +
gidz_general(z0, y1, x1, ny, nx)) +
1554 a * (1 - b) * (1 -
c) * *(phiin +
gidz_general(z0, y0, x1, ny, nx)) +
1555 (1 - a) * b * (1 -
c) * *(phiin +
gidz_general(z0, y1, x0, ny, nx)) +
1556 (1 - a) * (1 -
b) * c * *(phiin +
gidz_general(z1, y0, x0, ny, nx)) +
1557 (1 - a) * (1 -
b) * (1 - c) * *(phiin +
gidz_general(z0, y0, x0, ny, nx));
1560 printf(
"lib_vio> Conversion to cubic lattice completed. \n");
1565 void write_xplor(
char *vol_file,
double pwidth,
double porigx,
double porigy,
1566 double porigz,
unsigned pextx,
unsigned pexty,
unsigned pextz,
double *pphi)
1569 long mxstart, mystart, mzstart, mxend, myend, mzend;
1570 unsigned indx, indy, indz;
1571 unsigned long count;
1572 unsigned extx2, exty2, extz2;
1573 double origx2, origy2, origz2;
1576 fout = fopen(vol_file,
"w");
1584 printf(
"lib_vio> Input grid not in register with origin of coordinate system.\n");
1585 printf(
"lib_vio> Data will be interpolated to fit crystallographic X-PLOR format.\n");
1587 interpolate_map(&pphi2, &extx2, &exty2, &extz2, &origx2, &origy2, &origz2,
1588 pwidth, pwidth, pwidth, pphi, pextx, pexty, pextz, porigx,
1589 porigy, porigz, pwidth, pwidth, pwidth);
1592 mxstart = (long)
floor((origx2 / pwidth) + 0.5);
1593 mystart = (long)
floor((origy2 / pwidth) + 0.5);
1594 mzstart = (long)
floor((origz2 / pwidth) + 0.5);
1595 mxend = mxstart + extx2 - 1;
1596 myend = mystart + exty2 - 1;
1597 mzend = mzstart + extz2 - 1;
1600 printf(
"lib_vio>\n");
1601 printf(
"lib_vio> Writing X-PLOR formatted (ASCII) volumetric map \n");
1603 fprintf(fout,
" 2 !NTITLE \n");
1604 fprintf(fout,
"REMARKS FILENAME=\"%s\" \n", vol_file);
1605 fprintf(fout,
"REMARKS created by the Situs lib_vio library\n");
1606 fprintf(fout,
"%8d%8ld%8ld%8d%8ld%8ld%8d%8ld%8ld\n", extx2, mxstart, mxend, exty2, mystart, myend, extz2, mzstart, mzend);
1607 fprintf(fout,
"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n", extx2 * pwidth, exty2 * pwidth, extz2 * pwidth, 90.0, 90.0, 90.0);
1609 for (indz = 0; indz < extz2; indz++) {
1612 for (indy = 0; indy < exty2; indy++)
for (indx = 0; indx < extx2; indx++) {
1613 if ((count + 1) % 6 == 0)
fprintf(fout,
"%12.5E \n", *(pphi2 +
gidz_general(indz, indy, indx, exty2, extx2)));
1617 if ((count) % 6 != 0)
fprintf(fout,
" \n");
1619 fprintf(fout,
"%8d\n", -9999);
1623 printf(
"lib_vio> X-PLOR map written to file %s \n", vol_file);
1624 printf(
"lib_vio> X-PLOR map indexing (counting from 0): \n");
1625 printf(
"lib_vio> NA = %8d (# of X intervals in unit cell) \n", extx2);
1626 printf(
"lib_vio> AMIN = %8ld (start index X) \n", mxstart);
1627 printf(
"lib_vio> AMAX = %8ld (end index X) \n", mxend);
1628 printf(
"lib_vio> NB = %8d (# of Y intervals in unit cell) \n", exty2);
1629 printf(
"lib_vio> BMIN = %8ld (start index Y) \n", mystart);
1630 printf(
"lib_vio> BMAX = %8ld (end index Y) \n", myend);
1631 printf(
"lib_vio> NC = %8d (# of Z intervals in unit cell) \n", extz2);
1632 printf(
"lib_vio> CMIN = %8ld (start index Z) \n", mzstart);
1633 printf(
"lib_vio> CMAX = %8ld (end index Z) \n", mzend);
1634 printf(
"lib_vio> X-PLOR unit cell (based on the extent of the input density map): \n");
1635 printf(
"lib_vio> A = %8.3f (unit cell dimension) \n", extx2 * pwidth);
1636 printf(
"lib_vio> B = %8.3f (unit cell dimension) \n", exty2 * pwidth);
1637 printf(
"lib_vio> C = %8.3f (unit cell dimension) \n", extz2 * pwidth);
1638 printf(
"lib_vio> ALPHA = %8.3f (unit cell angle) \n", 90.0);
1639 printf(
"lib_vio> BETA = %8.3f (unit cell angle) \n", 90.0);
1640 printf(
"lib_vio> GAMMA = %8.3f (unit cell angle) \n", 90.0);
1646 void write_mrc(
int automode,
char *vol_file,
double pwidth,
double porigx,
1647 double porigy,
double porigz,
unsigned pextx,
unsigned pexty,
1648 unsigned pextz,
double *pphi)
1651 long nxstart, nystart, nzstart;
1652 float xorigin, yorigin, zorigin;
1653 unsigned long count, pnvox;
1654 long nx, ny, nz, mx, my, mz;
1657 float xlen, ylen, zlen, alpha,
beta,
gamma;
1658 long mapc, mapr, maps;
1659 long ispg = 1, nlabl = 0;
1662 float fmax, fmin, fav, fsig;
1663 double dshift, dscale;
1664 char dummychar =
'\0';
1669 float skwmat11 = 1.0f, skwmat12 = 0.0f, skwmat13 = 0.0f, skwmat21 = 0.0f;
1670 float skwmat22 = 1.0f, skwmat23 = 0.0f, skwmat31 = 0.0f, skwmat32 = 0.0f;
1671 float skwmat33 = 1.0f, skwtrn1 = 0.0f, skwtrn2 = 0.0f, skwtrn3 = 0.0f;
1672 char endchar[4] = {1, 1, 0, 0};
1674 fout = fopen(vol_file,
"wb");
1680 pnvox = pextx * pexty * pextz;
1688 if (*((
int *)endchar) > 65536) {
1724 nxstart = (long)
floor((porigx / pwidth) + 0.5);
1725 nystart = (long)
floor((porigy / pwidth) + 0.5);
1726 nzstart = (long)
floor((porigz / pwidth) + 0.5);
1730 if (automode == 0) {
1731 printf(
"lib_vio> Manual override of MRC / CCP4 header fields.\n");
1732 printf(
"lib_vio> \n");
1733 printf(
"lib_vio> The currently assigned data type is 32-bit floats (MODE 2)\n");
1734 printf(
"lib_vio> Select one of the following options: \n");
1735 printf(
"lib_vio> 1: keep MODE 2, 32-bit floats \n");
1736 printf(
"lib_vio> 2: MODE 0, signed 8-bit (range -127...128 as in CCP4), clip any out-of range\n");
1737 printf(
"lib_vio> 3: MODE 0, signed 8-bit (range -127...128 as in CCP4), scale to fit\n");
1738 printf(
"lib_vio> 4: MODE 0, signed 8-bit (range -127...128 as in CCP4), shift and scale to fit\n");
1739 printf(
"lib_vio> 5: MODE 0, unsigned 8-bit (range 0...255 as in old MRC), clip any out-of range\n");
1740 printf(
"lib_vio> 6: MODE 0, unsigned 8-bit (range 0...255 as in old MRC), scale pos., clip neg.\n");
1741 printf(
"lib_vio> 7: MODE 0, unsigned 8-bit (range 0...255 as in old MRC), shift and scale to fit\n");
1742 printf(
"lib_vio> \n");
1743 printf(
"lib_vio> Enter selection number: ");
1745 if (modesel < 1 && modesel > 7) {
1748 printf(
"lib_vio> Did not recognize value, assuming MODE 2, 32-bit floats. \n");
1749 }
else if (modesel > 1) mode = 0;
1751 printf(
"lib_vio> Note: The written data is not affected by changing the following header values.\n");
1752 printf(
"lib_vio> The currently assigned CCP4-style NC field (# columns) is %ld \n", nx);
1753 printf(
"lib_vio> Enter the same or a new value: ");
1755 printf(
"lib_vio> The currently assigned CCP4-style NR field (# rows) is %ld \n", ny);
1756 printf(
"lib_vio> Enter the same or a new value: ");
1758 printf(
"lib_vio> The currently assigned CCP4-style NS field (# sections) is %ld \n", nz);
1759 printf(
"lib_vio> Enter the same or a new value: ");
1761 if (((
unsigned long)nx * (
unsigned long)ny * (
unsigned long)nz) != pnvox) {
1762 printf(
"lib_vio> \n");
1763 fprintf(stderr,
"lib_vio> Warning: NC * NR * NS does not match the total number of voxels in the map.\n");
1764 printf(
"lib_vio> \n");
1766 printf(
"lib_vio> The currently assigned CCP4-style NCSTART field (column start index) is %ld \n", nxstart);
1768 printf(
"lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1770 printf(
"lib_vio> Enter the same or a new value: ");
1772 printf(
"lib_vio> The currently assigned CCP4-style NRSTART field (row start index) is %ld \n", nystart);
1774 printf(
"lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1776 printf(
"lib_vio> Enter the same or a new value: ");
1778 printf(
"lib_vio> The currently assigned CCP4-style NSSTART field (section start index) is %ld \n", nzstart);
1780 printf(
"lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1782 printf(
"lib_vio> Enter the same or a new value: ");
1784 printf(
"lib_vio> The currently assigned MX field (# of X intervals in the unit cell) is %ld \n", mx);
1785 printf(
"lib_vio> Enter the same or a new value: ");
1787 printf(
"lib_vio> The currently assigned MY field (# of Y intervals in the unit cell) is %ld \n", my);
1788 printf(
"lib_vio> Enter the same or a new value: ");
1790 printf(
"lib_vio> The currently assigned MZ field (# of Z intervals in the unit cell) is %ld \n", mz);
1791 printf(
"lib_vio> Enter the same or a new value: ");
1793 printf(
"lib_vio> The currently assigned X unit cell dimension is %f Angstrom\n", xlen);
1794 printf(
"lib_vio> Enter the same or a new value: ");
1796 printf(
"lib_vio> The currently assigned Y unit cell dimension is %f Angstrom\n", ylen);
1797 printf(
"lib_vio> Enter the same or a new value: ");
1799 printf(
"lib_vio> The currently assigned Z unit cell dimension is %f Angstrom\n", zlen);
1800 printf(
"lib_vio> Enter the same or a new value: ");
1802 printf(
"lib_vio> The currently assigned unit cell angle alpha is %f degrees \n", alpha);
1803 printf(
"lib_vio> Enter the same or a new value: ");
1805 printf(
"lib_vio> The currently assigned unit cell angle beta is %f degrees \n", beta);
1806 printf(
"lib_vio> Enter the same or a new value: ");
1808 printf(
"lib_vio> The currently assigned unit cell angle gamma is %f degrees \n", gamma);
1809 printf(
"lib_vio> Enter the same or a new value: ");
1811 printf(
"lib_vio> The currently assigned MAPC field (axis order) is %ld \n", mapc);
1812 printf(
"lib_vio> Enter the same or a new value: ");
1815 if (mapc == 1 || mapc == 2 || mapc == 3) wmaperr = 0;
1820 printf(
"lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1822 printf(
"lib_vio> The currently assigned MAPR field (axis order) is %ld \n", mapr);
1823 printf(
"lib_vio> Enter the same or a new value: ");
1826 if (mapc == 1 && mapr == 2) wmaperr = 0;
1827 if (mapc == 1 && mapr == 3) wmaperr = 0;
1828 if (mapc == 2 && mapr == 1) wmaperr = 0;
1829 if (mapc == 2 && mapr == 3) wmaperr = 0;
1830 if (mapc == 3 && mapr == 1) wmaperr = 0;
1831 if (mapc == 3 && mapr == 2) wmaperr = 0;
1836 printf(
"lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1838 printf(
"lib_vio> The currently assigned MAPS field (axis order) is %ld \n", maps);
1839 printf(
"lib_vio> Enter the same or a new value: ");
1842 if (mapc == 1 && mapr == 2 && maps == 3) wmaperr = 0;
1843 if (mapc == 1 && mapr == 3 && maps == 2) wmaperr = 0;
1844 if (mapc == 2 && mapr == 1 && maps == 3) wmaperr = 0;
1845 if (mapc == 2 && mapr == 3 && maps == 1) wmaperr = 0;
1846 if (mapc == 3 && mapr == 1 && maps == 2) wmaperr = 0;
1847 if (mapc == 3 && mapr == 2 && maps == 1) wmaperr = 0;
1852 printf(
"lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1856 printf(
"lib_vio> The currently assigned MRC2000-style XORIGIN field is %f Angstrom\n", xorigin);
1857 printf(
"lib_vio> Enter the same or a new value: ");
1859 printf(
"lib_vio> The currently assigned MRC2000-style YORIGIN field is %f Angstrom\n", yorigin);
1860 printf(
"lib_vio> Enter the same or a new value: ");
1862 printf(
"lib_vio> The currently assigned MRC2000-style ZORIGIN field is %f Angstrom\n", zorigin);
1863 printf(
"lib_vio> Enter the same or a new value: ");
1872 if (
clipped(pphi, pnvox, 127.0, -128.0))
1873 printf(
"lib_vio> Density values were clipped to fit signed 8-bit data format.\n");
1875 printf(
"lib_vio> Density values already fit signed 8-bit data format, no clipping necessary.\n");
1878 dscale = dmax / 127.0;
1879 if (dmin / -128.0 > dscale)
1880 dscale = dmin / -128.0;
1881 printf(
"lib_vio> Density values will be rescaled to fit signed 8-bit data format.\n");
1882 printf(
"lib_vio> Scaling factor (by which values will be divided): %f .\n", dscale);
1886 dscale = (dmax -
dmin) / 255.0;
1887 dshift = (dmax +
dmin) / 2.0 + 0.5 * dscale;
1888 printf(
"lib_vio> Density values will be centered and rescaled to fit signed 8-bit data format.\n");
1889 printf(
"lib_vio> Center offset value (will be added to map): %f .\n", -dshift);
1890 printf(
"lib_vio> Scaling factor (by which values will be divided): %f .\n", dscale);
1895 if (
clipped(pphi, pnvox, 255.0, 0.0)) printf(
"lib_vio> Density values were clipped to fit unsigned 8-bit data format.\n");
1896 else printf(
"lib_vio> Density values already fit unsigned 8-bit data format, no clipping necessary.\n");
1899 dscale = dmax / 255.0;
1900 printf(
"lib_vio> Density values will be rescaled to fit unsigned 8-bit data format.\n");
1901 printf(
"lib_vio> Scaling factor (by which values will be divided): %f .\n", dscale);
1903 if (
clipped(pphi, pnvox, 256.0, 0.0)) printf(
"lib_vio> Negative density values were clipped to fit unsigned 8-bit data format.\n");
1906 dscale = (dmax -
dmin) / 255.0;
1907 dshift = (dmax +
dmin) / 2.0 - 127.5 * dscale;
1908 printf(
"lib_vio> Density values will be centered and rescaled to fit unsigned 8-bit data format.\n");
1909 printf(
"lib_vio> Center offset value (will be added to map): %f .\n", -dshift);
1910 printf(
"lib_vio> Scaling factor (by which values will be divided): %f .\n", dscale);
1927 printf(
"lib_vio> Converting the unsigned 8-bit (0...255) values to signed (-128...127) for storage (ISO/IEC 10967).\n");
1928 printf(
"lib_vio> The conversion subtracts 256 from density values above 127.\n");
1929 printf(
"lib_vio> Check your software documentation if the `unsigned` convention is supported for MODE 0 maps.\n");
1930 for (count = 0; count < pnvox; count++) if (*(pphi + count) >= 127.5) *(pphi + count) -= 256.0;
1937 printf(
"lib_vio> Writing MRC / CCP4 (binary) volumetric map \n");
1939 if (fwrite(&nx, 4, 1, fout) != 1) wmaperr = 1;
1940 if (fwrite(&ny, 4, 1, fout) != 1) wmaperr = 1;
1941 if (fwrite(&nz, 4, 1, fout) != 1) wmaperr = 1;
1942 if (fwrite(&mode, 4, 1, fout) != 1) wmaperr = 1;
1943 if (fwrite(&nxstart, 4, 1, fout) != 1) wmaperr = 1;
1944 if (fwrite(&nystart, 4, 1, fout) != 1) wmaperr = 1;
1945 if (fwrite(&nzstart, 4, 1, fout) != 1) wmaperr = 1;
1946 if (fwrite(&mx, 4, 1, fout) != 1) wmaperr = 1;
1947 if (fwrite(&my, 4, 1, fout) != 1) wmaperr = 1;
1948 if (fwrite(&mz, 4, 1, fout) != 1) wmaperr = 1;
1949 if (fwrite(&xlen, 4, 1, fout) != 1) wmaperr = 1;
1950 if (fwrite(&ylen, 4, 1, fout) != 1) wmaperr = 1;
1951 if (fwrite(&zlen, 4, 1, fout) != 1) wmaperr = 1;
1952 if (fwrite(&alpha, 4, 1, fout) != 1) wmaperr = 1;
1953 if (fwrite(&beta, 4, 1, fout) != 1) wmaperr = 1;
1954 if (fwrite(&gamma, 4, 1, fout) != 1) wmaperr = 1;
1955 if (fwrite(&mapc, 4, 1, fout) != 1) wmaperr = 1;
1956 if (fwrite(&mapr, 4, 1, fout) != 1) wmaperr = 1;
1957 if (fwrite(&maps, 4, 1, fout) != 1) wmaperr = 1;
1958 if (fwrite(&fmin, 4, 1, fout) != 1) wmaperr = 1;
1959 if (fwrite(&fmax, 4, 1, fout) != 1) wmaperr = 1;
1960 if (fwrite(&fav, 4, 1, fout) != 1) wmaperr = 1;
1961 if (fwrite(&ispg, 4, 1, fout) != 1) wmaperr = 1;
1962 if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
1963 if (fwrite(&lskflg, 4, 1, fout) != 1) wmaperr = 1;
1964 if (fwrite(&skwmat11, 4, 1, fout) != 1) wmaperr = 1;
1965 if (fwrite(&skwmat12, 4, 1, fout) != 1) wmaperr = 1;
1966 if (fwrite(&skwmat13, 4, 1, fout) != 1) wmaperr = 1;
1967 if (fwrite(&skwmat21, 4, 1, fout) != 1) wmaperr = 1;
1968 if (fwrite(&skwmat22, 4, 1, fout) != 1) wmaperr = 1;
1969 if (fwrite(&skwmat23, 4, 1, fout) != 1) wmaperr = 1;
1970 if (fwrite(&skwmat31, 4, 1, fout) != 1) wmaperr = 1;
1971 if (fwrite(&skwmat32, 4, 1, fout) != 1) wmaperr = 1;
1972 if (fwrite(&skwmat33, 4, 1, fout) != 1) wmaperr = 1;
1973 if (fwrite(&skwtrn1, 4, 1, fout) != 1) wmaperr = 1;
1974 if (fwrite(&skwtrn2, 4, 1, fout) != 1) wmaperr = 1;
1975 if (fwrite(&skwtrn3, 4, 1, fout) != 1) wmaperr = 1;
1976 for (i = 38; i < 50; ++
i)
if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
1977 if (fwrite(&xorigin, 4, 1, fout) != 1) wmaperr = 1;
1978 if (fwrite(&yorigin, 4, 1, fout) != 1) wmaperr = 1;
1979 if (fwrite(&zorigin, 4, 1, fout) != 1) wmaperr = 1;
1980 if (fwrite(mapstring, 4, 1, fout) != 1) wmaperr = 1;
1981 if (fwrite(machstring, 4, 1, fout) != 1) wmaperr = 1;
1982 if (fwrite(&fsig, 4, 1, fout) != 1) wmaperr = 1;
1983 if (fwrite(&nlabl, 4, 1, fout) != 1) wmaperr = 1;
1984 for (i = 0; i < 800; ++
i)
if (fwrite(&dummychar,
sizeof(
char), 1, fout) != 1) wmaperr = 1;
1989 for (count = 0; count < pnvox; count++) {
1990 dummy =
floor(*(pphi + count) + 0.5);
1991 if (dummy < -128 || dummy > 127) wmaperr = 1;
1992 dummychar = (char)dummy;
1993 if (fwrite(&dummychar, 1, 1, fout) != 1) wmaperr = 1;
1997 for (count = 0; count < pnvox; count++) {
1998 dummy = *(pphi + count);
1999 if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2012 printf(
"lib_vio> Volumetric map written to file %s \n", vol_file);
2013 printf(
"lib_vio> Header information: \n");
2014 printf(
"lib_vio> NC = %8ld (# columns)\n", nx);
2015 printf(
"lib_vio> NR = %8ld (# rows)\n", ny);
2016 printf(
"lib_vio> NS = %8ld (# sections)\n", nz);
2017 printf(
"lib_vio> MODE = %8ld (data type: 0: 8-bit char, 2: 32-bit float)\n", mode);
2018 printf(
"lib_vio> NCSTART = %8ld (index of first column, counting from 0)\n", nxstart);
2019 printf(
"lib_vio> NRSTART = %8ld (index of first row, counting from 0)\n", nystart);
2020 printf(
"lib_vio> NSSTART = %8ld (index of first section, counting from 0)\n", nzstart);
2021 printf(
"lib_vio> MX = %8ld (# of X intervals in unit cell)\n", mx);
2022 printf(
"lib_vio> MY = %8ld (# of Y intervals in unit cell)\n", my);
2023 printf(
"lib_vio> MZ = %8ld (# of Z intervals in unit cell)\n", mz);
2024 printf(
"lib_vio> X length = %8.3f (unit cell dimension)\n", xlen);
2025 printf(
"lib_vio> Y length = %8.3f (unit cell dimension)\n", ylen);
2026 printf(
"lib_vio> Z length = %8.3f (unit cell dimension)\n", zlen);
2027 printf(
"lib_vio> Alpha = %8.3f (unit cell angle)\n", alpha);
2028 printf(
"lib_vio> Beta = %8.3f (unit cell angle)\n", beta);
2029 printf(
"lib_vio> Gamma = %8.3f (unit cell angle)\n", gamma);
2030 printf(
"lib_vio> MAPC = %8ld (columns axis: 1=X,2=Y,3=Z)\n", mapc);
2031 printf(
"lib_vio> MAPR = %8ld (rows axis: 1=X,2=Y,3=Z)\n", mapr);
2032 printf(
"lib_vio> MAPS = %8ld (sections axis: 1=X,2=Y,3=Z)\n", maps);
2033 printf(
"lib_vio> DMIN = %8.3f (minimum density value)\n", fmin);
2034 printf(
"lib_vio> DMAX = %8.3f (maximum density value)\n", fmax);
2035 printf(
"lib_vio> DMEAN = %8.3f (mean density value)\n", fav);
2036 printf(
"lib_vio> ISPG = %8ld (space group number)\n", ispg);
2037 printf(
"lib_vio> NSYMBT = %8d (# bytes used for storing symmetry operators)\n", 0);
2038 printf(
"lib_vio> LSKFLG = %8ld (skew matrix flag: 0:none, 1:follows)\n", lskflg);
2040 printf(
"lib_vio> S11 = %8.3f (skew matrix element 11)\n", skwmat11);
2041 printf(
"lib_vio> S12 = %8.3f (skew matrix element 12)\n", skwmat12);
2042 printf(
"lib_vio> S13 = %8.3f (skew matrix element 13)\n", skwmat13);
2043 printf(
"lib_vio> S21 = %8.3f (skew matrix element 21)\n", skwmat21);
2044 printf(
"lib_vio> S22 = %8.3f (skew matrix element 22)\n", skwmat22);
2045 printf(
"lib_vio> S23 = %8.3f (skew matrix element 23)\n", skwmat23);
2046 printf(
"lib_vio> S31 = %8.3f (skew matrix element 31)\n", skwmat31);
2047 printf(
"lib_vio> S32 = %8.3f (skew matrix element 32)\n", skwmat32);
2048 printf(
"lib_vio> S33 = %8.3f (skew matrix element 33)\n", skwmat33);
2049 printf(
"lib_vio> T1 = %8.3f (skew translation element 1)\n", skwtrn1);
2050 printf(
"lib_vio> T2 = %8.3f (skew translation element 2)\n", skwtrn2);
2051 printf(
"lib_vio> T3 = %8.3f (skew translation element 3)\n", skwtrn3);
2053 printf(
"lib_vio> XORIGIN = %8.3f (X origin - MRC2000 only)\n", xorigin);
2054 printf(
"lib_vio> YORIGIN = %8.3f (Y origin - MRC2000 only)\n", yorigin);
2055 printf(
"lib_vio> ZORIGIN = %8.3f (Z origin - MRC2000 only)\n", zorigin);
2056 printf(
"lib_vio> MAP = %c%c%c%c (map string)\n", mapstring[0], mapstring[1], mapstring[2], mapstring[3]);
2057 printf(
"lib_vio> MACHST = %d %d %d %d (machine stamp)\n", machstring[0], machstring[1], machstring[2], machstring[3]);
2058 printf(
"lib_vio> RMS = %8.3f (density rmsd)\n", fsig);
2062 printf(
"lib_vio> Input grid not in register with origin of coordinate system.\n");
2063 printf(
"lib_vio> The origin information was saved only in the MRC2000 format fields.\n");
2064 printf(
"lib_vio> The CCP4 start indexing was set to zero.\n");
2065 printf(
"lib_vio> To invoke a crystallographic CCP4 indexing, you can force an interpolation \n");
2066 printf(
"lib_vio> by converting to an intermediate X-PLOR map using the map2map tool. \n");
2072 void write_spider(
char *vol_file,
double pwidth,
double porigx,
double porigy,
2073 double porigz,
unsigned pextx,
unsigned pexty,
unsigned pextz,
2077 unsigned long count, pnvox;
2078 float nslice, nrow, iform, imami, nsam, headrec;
2080 float fmax, fmin, fav, fsig;
2081 float iangle, phi,
theta, xoff, yoff, zoff, scale, labbyt, lenbyt, irec;
2082 float istack, inuse, maxim, kangle, phi1, theta1, psi1, phi2, theta2, psi2;
2085 int i, headlen, wmaperr;
2087 fout = fopen(vol_file,
"wb");
2093 pnvox = pextx * pexty * pextz;
2127 headrec = ceil(256.0
f / (pextx * 1.0
f));
2128 lenbyt = nsam * 4.0f;
2129 labbyt = headrec * lenbyt;
2130 irec = nslice * nrow + headrec;
2131 headlen = (int)(headrec * nsam);
2132 printf(
"lib_vio>\n");
2133 printf(
"lib_vio> Writing SPIDER (binary) volumetric map... \n");
2135 if (fwrite(&nslice, 4, 1, fout) != 1) wmaperr = 1;
2136 if (fwrite(&nrow, 4, 1, fout) != 1) wmaperr = 1;
2137 if (fwrite(&irec, 4, 1, fout) != 1) wmaperr = 1;
2138 if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2139 if (fwrite(&iform, 4, 1, fout) != 1) wmaperr = 1;
2140 if (fwrite(&imami, 4, 1, fout) != 1) wmaperr = 1;
2141 if (fwrite(&fmax, 4, 1, fout) != 1) wmaperr = 1;
2142 if (fwrite(&fmin, 4, 1, fout) != 1) wmaperr = 1;
2143 if (fwrite(&fav, 4, 1, fout) != 1) wmaperr = 1;
2144 if (fwrite(&fsig, 4, 1, fout) != 1) wmaperr = 1;
2145 if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2146 if (fwrite(&nsam, 4, 1, fout) != 1) wmaperr = 1;
2147 if (fwrite(&headrec, 4, 1, fout) != 1) wmaperr = 1;
2148 if (fwrite(&iangle, 4, 1, fout) != 1) wmaperr = 1;
2149 if (fwrite(&phi, 4, 1, fout) != 1) wmaperr = 1;
2150 if (fwrite(&theta, 4, 1, fout) != 1) wmaperr = 1;
2151 if (fwrite(&gamma, 4, 1, fout) != 1) wmaperr = 1;
2152 if (fwrite(&xoff, 4, 1, fout) != 1) wmaperr = 1;
2153 if (fwrite(&yoff, 4, 1, fout) != 1) wmaperr = 1;
2154 if (fwrite(&zoff, 4, 1, fout) != 1) wmaperr = 1;
2155 if (fwrite(&scale, 4, 1, fout) != 1) wmaperr = 1;
2156 if (fwrite(&labbyt, 4, 1, fout) != 1) wmaperr = 1;
2157 if (fwrite(&lenbyt, 4, 1, fout) != 1) wmaperr = 1;
2158 if (fwrite(&istack, 4, 1, fout) != 1) wmaperr = 1;
2159 if (fwrite(&inuse, 4, 1, fout) != 1) wmaperr = 1;
2160 if (fwrite(&maxim, 4, 1, fout) != 1) wmaperr = 1;
2161 for (i = 0; i < 4; ++
i)
if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2162 if (fwrite(&kangle, 4, 1, fout) != 1) wmaperr = 1;
2163 if (fwrite(&phi1, 4, 1, fout) != 1) wmaperr = 1;
2164 if (fwrite(&theta1, 4, 1, fout) != 1) wmaperr = 1;
2165 if (fwrite(&psi1, 4, 1, fout) != 1) wmaperr = 1;
2166 if (fwrite(&phi2, 4, 1, fout) != 1) wmaperr = 1;
2167 if (fwrite(&theta2, 4, 1, fout) != 1) wmaperr = 1;
2168 if (fwrite(&psi2, 4, 1, fout) != 1) wmaperr = 1;
2169 for (i = 0; i < (headlen - 37); ++
i)
if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2170 for (count = 0; count < pnvox; count++) {
2171 dummy = *(pphi + count);
2172 if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2181 printf(
"lib_vio> SPIDER map written to file %s \n", vol_file);
2182 printf(
"lib_vio> SPIDER map indexing: \n");
2183 printf(
"lib_vio> NSLICE = %8.f (# sections)\n", nslice);
2184 printf(
"lib_vio> NROW = %8.f (# rows)\n", nrow);
2185 printf(
"lib_vio> IFORM = %8.f (file type specifier)\n", iform);
2186 printf(
"lib_vio> IMAMI = %8.f (flag: =1 the maximum and minimum values are computed)\n", imami);
2187 printf(
"lib_vio> FMAX = %8.3f (maximum density value)\n", fmax);
2188 printf(
"lib_vio> FMIN = %8.3f (minimum density value)\n", fmin);
2189 printf(
"lib_vio> AV = %8.3f (average density value)\n", fav);
2190 printf(
"lib_vio> SIG = %8.3f (standard deviation of density distribution)\n", fsig);
2191 printf(
"lib_vio> NSAM = %8.f (# columns)\n", nsam);
2192 printf(
"lib_vio> HEADREC = %8.f (number of records in file header)\n", headrec);
2193 printf(
"lib_vio> IANGLE = %8.f (flag: =1 tilt angles filled)\n", iangle);
2194 printf(
"lib_vio> PHI = %8.3f (tilt angle)\n", phi);
2195 printf(
"lib_vio> THETA = %8.3f (tilt angle)\n", theta);
2196 printf(
"lib_vio> GAMMA = %8.3f (tilt angle)\n", gamma);
2197 printf(
"lib_vio> XOFF = %8.3f (X offset)\n", xoff);
2198 printf(
"lib_vio> YOFF = %8.3f (Y offset)\n", yoff);
2199 printf(
"lib_vio> ZOFF = %8.3f (Z offset)\n", zoff);
2200 printf(
"lib_vio> SCALE = %8.3f (scale factor)\n", scale);
2201 printf(
"lib_vio> LABBYT = %8.f (total number of bytes in header)\n", labbyt);
2202 printf(
"lib_vio> LENBYT = %8.f (record length in bytes)\n", lenbyt);
2203 printf(
"lib_vio> ISTACK = %8.f (flag; file contains a stack of images)\n", istack);
2204 printf(
"lib_vio> INUSE = %8.f (flag; this image in stack is used)\n", inuse);
2205 printf(
"lib_vio> MAXIM = %8.f (maximum image used in stack)\n", maxim);
2206 printf(
"lib_vio> KANGLE = %8.f (flag; additional angles set)\n", kangle);
2207 printf(
"lib_vio> PHI1 = %8.3f (additional rotation)\n", phi1);
2208 printf(
"lib_vio> THETA1 = %8.3f (additional rotation)\n", theta1);
2209 printf(
"lib_vio> PSI1 = %8.3f (additional rotation)\n", psi1);
2210 printf(
"lib_vio> PHI2 = %8.3f (additional rotation)\n", phi2);
2211 printf(
"lib_vio> THETA2 = %8.3f (additional rotation)\n", theta2);
2212 printf(
"lib_vio> PSI2 = %8.3f (additional rotation)\n", psi2);
2213 printf(
"lib_vio> Warning: The voxel spacing %f has not been saved to the SPIDER map.\n", pwidth);
2214 printf(
"lib_vio> Warning: The map origin %f,%f,%f (first voxel position) has not been saved to the SPIDER map.\n", porigx, porigy, porigz);
int test_registration(float origx, float origy, float origz, float width)
void xplor_skip_to_number(FILE **fin, char **nextline)
void error_xplor_file_map_section(int error_number, const char *program)
void write_spider(char *vol_file, double pwidth, double porigx, double porigy, double porigz, unsigned pextx, unsigned pexty, unsigned pextz, double *pphi)
double rms(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=nullptr, MultidimArray< double > *Contributions=nullptr)
void cp_vect_destroy(double **pvect1, double **pvect2, unsigned long len)
__host__ __device__ float2 floor(const float2 v)
int read_float(float *currfloat, FILE *fin, int swap)
void free_vect_and_zero_ptr(void **pvect)
int read_int(int *currlong, FILE *fin, int swap)
void error_EOF(int error_number, const char *program)
void error_xplor_file_indexing(int error_number, const char *program)
void error_option(int error_number, const char *program)
void assert_cubic_map(int orom, int cubic, double alpha, double beta, double gamma, double widthx, double widthy, double widthz, unsigned nx, unsigned ny, unsigned nz, int nxstart, int nystart, int nzstart, double xorigin, double yorigin, double zorigin, double *pwidth, double *porigx, double *porigy, double *porigz, unsigned *pextx, unsigned *pexty, unsigned *pextz, double **pphi)
void error_write_filename(int error_number, const char *program)
double beta(const double a, const double b)
void error_xplor_file_map_section_number(int error_number, const char *program)
void sqrt(Image< double > &op)
int test_spider(char *vol_file, int swap)
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
void permute_map(int ordermode, unsigned nc, unsigned nr, unsigned ns, unsigned *nx, unsigned *ny, unsigned *nz, int ncstart, int nrstart, int nsstart, int *nxstart, int *nystart, int *nzstart, double *phi, double **pphi)
void read_vol(char *vol_file, double *width, double *origx, double *origy, double *origz, unsigned *extx, unsigned *exty, unsigned *extz, double **phi)
unsigned long permuted_index(int ordermode, unsigned long count, unsigned nc, unsigned nr, unsigned ns)
void calc_map_info(double *phi, unsigned long nvox, double *maxdens, double *mindens, double *ave, double *sig)
int read_short_float(float *currfloat, FILE *fin, int swap)
unsigned long count_floats(FILE **fin)
void interpolate_map(double **outmap, unsigned *out_extx, unsigned *out_exty, unsigned *out_extz, double *out_origx, double *out_origy, double *out_origz, double out_widthx, double out_widthy, double out_widthz, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double in_widthx, double in_widthy, double in_widthz)
void error_open_filename(int error_number, const char *program, char *argv)
int set_origin_get_mode(int nxstart, int nystart, int nzstart, double widthx, double widthy, double widthz, double xorigin, double yorigin, double zorigin, double *origx, double *origy, double *origz)
void interpolate_skewed_map_to_cubic(double **pphiout, unsigned *pextx, unsigned *pexty, unsigned *pextz, double *porigx, double *porigy, double *porigz, double *pwidth, double *phiin, unsigned nx, unsigned ny, unsigned nz, int nxstart, int nystart, int nzstart, double xorigin, double yorigin, double zorigin, double widthx, double widthy, double widthz, double alpha, double beta, double gamma)
int read_char(char *currchar, FILE *fin)
void error_fscanf(const char *program, const char *file)
quaternion_type< T > normalize(quaternion_type< T > q)
void error_sqrt_negative(int error_number, const char *program)
void read_xplor(char *vol_file, int *orom, int *cubic, double *widthx, double *widthy, double *widthz, double *alpha, double *beta, double *gamma, int *nxstart, int *nystart, int *nzstart, unsigned *extx, unsigned *exty, unsigned *extz, double **fphi)
void error_divide_zero(int error_number, const char *program)
void floatshift(double *phi, unsigned long nvox, double dens)
void error_axis_assignment(int error_number, const char *program)
void error_unreadable_file_long(int error_number, const char *program, const char *filename)
int test_mrc(char *vol_file, int swap)
void write_xplor(char *vol_file, double pwidth, double porigx, double porigy, double porigz, unsigned pextx, unsigned pexty, unsigned pextz, double *pphi)
void error_EOF_ZYX_mode(int error_number, const char *program)
void error_unreadable_file_short(int error_number, const char *program, const char *filename)
int clipped(double *phi, unsigned long nvox, double max, double min)
void read_spider(char *vol_file, unsigned *extx, unsigned *exty, unsigned *extz, double **fphi)
void read_situs(char *vol_file, double *width, double *origx, double *origy, double *origz, unsigned *extx, unsigned *exty, unsigned *extz, double **phi)
void error_xplor_maker(const char *program)
void read_mrc(char *vol_file, int *orom, int *cubic, int *ordermode, unsigned *nc, unsigned *nr, unsigned *ns, int *ncstart, int *nrstart, int *nsstart, double *widthx, double *widthy, double *widthz, double *xorigin, double *yorigin, double *zorigin, double *alpha, double *beta, double *gamma, double **fphi)
void error_file_convert(int error_number, const char *program, const char *filename)
void do_vect(double **pvect, unsigned long len)
void read_ascii(char *vol_file, unsigned long nvox, double **fphi)
void * alloc_vect(unsigned int n, size_t elem_size)
void write_situs(char *vol_file, double width, double origx, double origy, double origz, unsigned extx, unsigned exty, unsigned extz, double *phi)
void write_mrc(int automode, char *vol_file, double pwidth, double porigx, double porigy, double porigz, unsigned pextx, unsigned pexty, unsigned pextz, double *pphi)
int read_float_empty(FILE *fin)
void permute_dimensions(int ordermode, unsigned nc, unsigned nr, unsigned ns, unsigned *nx, unsigned *ny, unsigned *nz, int ncstart, int nrstart, int nsstart, int *nxstart, int *nystart, int *nzstart)
fprintf(glob_prnt.io, "\)
int have_situs_suffix(char *vol_file)
int test_situs_header_and_suffix(char *vol_file)
void error_xplor_file_unit_cell(int error_number, const char *program)
int read_char_float(float *currfloat, FILE *fin)
void error_spider_header(int error_number, const char *program)
void write_vol(char *vol_file, double width, double origx, double origy, double origz, unsigned extx, unsigned exty, unsigned extz, double *phi)
void permute_print_info(int ordermode, unsigned nc, unsigned nr, unsigned ns, unsigned nx, unsigned ny, unsigned nz, int ncstart, int nrstart, int nsstart, int nxstart, int nystart, int nzstart)
void error_file_header(int error_number, const char *program, const char *filename)
void dump_binary_and_exit(char *in_file, char *out_file, int swap)
int test_situs(char *vol_file)
void error_file_float_mode(int error_number, const char *program, const char *filename)