Xmipp  v3.23.11-Nereus
Functions
lib_vio.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void read_vol (char *, double *, double *, double *, double *, unsigned *, unsigned *, unsigned *, double **)
 
void write_vol (char *, double, double, double, double, unsigned, unsigned, unsigned, double *)
 
void read_situs (char *, double *, double *, double *, double *, unsigned *, unsigned *, unsigned *, double **)
 
void write_situs (char *, double, double, double, double, unsigned, unsigned, unsigned, double *)
 
void read_ascii (char *, unsigned long, double **)
 
void read_xplor (char *, int *, int *, double *, double *, double *, double *, double *, double *, int *, int *, int *, unsigned *, unsigned *, unsigned *, double **)
 
void xplor_skip_to_number (FILE **, char **)
 
void read_mrc (char *, int *, int *, int *, unsigned *, unsigned *, unsigned *, int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double **)
 
void read_spider (char *, unsigned *, unsigned *, unsigned *, double **)
 
void dump_binary_and_exit (char *, char *, int)
 
int read_float (float *, FILE *, int)
 
int read_short_float (float *, FILE *, int)
 
int read_float_empty (FILE *)
 
int read_char_float (float *, FILE *)
 
int read_char (char *, FILE *)
 
int read_int (int *, FILE *, int)
 
unsigned long count_floats (FILE **)
 
int test_registration (float, float, float, float)
 
int test_situs (char *)
 
int have_situs_suffix (char *)
 
int test_situs_header_and_suffix (char *)
 
int test_mrc (char *, int)
 
int test_spider (char *, int)
 
unsigned long permuted_index (int, unsigned long, unsigned, unsigned, unsigned)
 
void permute_map (int, unsigned, unsigned, unsigned, unsigned *, unsigned *, unsigned *, int, int, int, int *, int *, int *, double *, double **)
 
void permute_dimensions (int, unsigned, unsigned, unsigned, unsigned *, unsigned *, unsigned *, int, int, int, int *, int *, int *)
 
void permute_print_info (int, unsigned, unsigned, unsigned, unsigned, unsigned, unsigned, int, int, int, int, int, int)
 
int set_origin_get_mode (int, int, int, double, double, double, double, double, double, double *, double *, double *)
 
void assert_cubic_map (int, int, double, double, double, double, double, double, unsigned, unsigned, unsigned, int, int, int, double, double, double, double *, double *, double *, double *, unsigned *, unsigned *, unsigned *, double **)
 
void interpolate_skewed_map_to_cubic (double **, unsigned *, unsigned *, unsigned *, double *, double *, double *, double *, double *, unsigned, unsigned, unsigned, int, int, int, double, double, double, double, double, double, double, double, double)
 
void write_xplor (char *, double, double, double, double, unsigned, unsigned, unsigned, double *)
 
void write_mrc (int, char *, double, double, double, double, unsigned, unsigned, unsigned, double *)
 
void write_spider (char *, double, double, double, double, unsigned, unsigned, unsigned, double *)
 

Function Documentation

◆ assert_cubic_map()

void assert_cubic_map ( int  ,
int  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
int  ,
int  ,
int  ,
double  ,
double  ,
double  ,
double *  ,
double *  ,
double *  ,
double *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double **   
)

Definition at line 1326 of file lib_vio.cpp.

1334 {
1335 
1336  unsigned long nvox;
1337  double *pp2;
1338  double origx, origy, origz;
1339 
1340  /* distinguish between cubic, orthonormal, and skewed */
1341  if (cubic) {
1342  set_origin_get_mode(nxstart, nystart, nzstart, widthx, widthy, widthz,
1343  xorigin, yorigin, zorigin, porigx, porigy, porigz);
1344  printf("lib_vio> Cubic lattice present. \n");
1345  *pwidth = widthx;
1346  *pextx = nx;
1347  *pexty = ny;
1348  *pextz = nz;
1349  } else {
1350  if (orom) {
1351  set_origin_get_mode(nxstart, nystart, nzstart, widthx, widthy, widthz,
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");
1355  /* for interpolation, find minimum voxel spacing, rounded to nearest 0.1 Angstrom for good measure */
1356  *pwidth = widthx;
1357  if (widthy < *pwidth) *pwidth = widthy;
1358  if (widthz < *pwidth) *pwidth = widthz;
1359  *pwidth = floor(*pwidth * 10.0 + 0.5) / 10.0;
1360  /* the new map will have the maximum box dimensions that do not exceed the old map */
1361  interpolate_map(&pp2, pextx, pexty, pextz, porigx, porigy, porigz,
1362  *pwidth, *pwidth, *pwidth, *pphi, nx, ny, nz, origx,
1363  origy, origz, widthx, widthy, widthz);
1364  nvox = *pextx * *pexty * *pextz;
1365  cp_vect_destroy(pphi, &pp2, nvox);
1366  } else {
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");
1369  /* here, the new map will have the minimum box dimensions that fully enclose the old map */
1370  interpolate_skewed_map_to_cubic(&pp2, pextx, pexty, pextz, porigx, porigy,
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;
1376  cp_vect_destroy(pphi, &pp2, nvox);
1377  }
1378  }
1379 
1380  /* check for registration of lattice with origin of coordinate system */
1381  if (test_registration(*porigx, *porigy, *porigz, *pwidth) == 0)
1382  fprintf(stderr, "lib_vio> Input grid not in register with origin of coordinate system.\n");
1383 }
test_registration
int test_registration(float origx, float origy, float origz, float width)
Definition: lib_vio.cpp:917
cp_vect_destroy
void cp_vect_destroy(double **pvect1, double **pvect2, unsigned long len)
Definition: lib_vec.cpp:51
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
interpolate_map
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)
Definition: lib_vwk.cpp:82
set_origin_get_mode
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)
Definition: lib_vio.cpp:1307
interpolate_skewed_map_to_cubic
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)
Definition: lib_vio.cpp:1388
fprintf
fprintf(glob_prnt.io, "\)

◆ count_floats()

unsigned long count_floats ( FILE **  )

Definition at line 902 of file lib_vio.cpp.

903 {
904  unsigned long finl = 0;
905 
906  rewind(*fin);
907  for (; ;) {
908  if (fgetc(*fin) == EOF) break;
909  ++finl;
910  }
911  rewind(*fin);
912  return finl / 4;
913 }

◆ dump_binary_and_exit()

void dump_binary_and_exit ( char *  ,
char *  ,
int   
)

Definition at line 774 of file lib_vio.cpp.

775 {
776  FILE *fin, *fout;
777  float *phi;
778  unsigned long nfloat;
779  unsigned long count;
780 
781  fin = fopen(in_file, "rb");
782  if (fin == NULL) {
783  error_open_filename(70910, "lib_vio", in_file);
784  }
785 
786  nfloat = count_floats(&fin);
787 
788  phi = (float *) alloc_vect(nfloat, sizeof(float));
789 
790  for (count = 0; count < nfloat; count++)
791  read_float(phi + count, fin, swap);
792 
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");
795  fclose(fin);
796 
797  fout = fopen(out_file, "w");
798  if (fout == NULL) {
799  error_open_filename(70940, "lib_vio", out_file);
800  }
801 
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));
806  }
807  fclose(fout);
808 
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");
812  exit(1);
813 }
read_float
int read_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:817
count_floats
unsigned long count_floats(FILE **fin)
Definition: lib_vio.cpp:902
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
alloc_vect
void * alloc_vect(unsigned int n, size_t elem_size)
Definition: lib_vec.cpp:71
fprintf
fprintf(glob_prnt.io, "\)

◆ have_situs_suffix()

int have_situs_suffix ( char *  )

Definition at line 981 of file lib_vio.cpp.

982 {
983 
984  int vl;
985  int situs_format = 0;
986 
987  vl = strlen(vol_file);
988  if (vl > 4) {
989  if (strstr(vol_file + (vl - 4), ".sit") != NULL) situs_format = 1;
990  if (strstr(vol_file + (vl - 4), ".SIT") != NULL) situs_format = 1;
991  }
992  if (vl > 6) {
993  if (strstr(vol_file + (vl - 6), ".situs") != NULL) situs_format = 1;
994  if (strstr(vol_file + (vl - 6), ".SITUS") != NULL) situs_format = 1;
995  }
996  return situs_format;
997 }

◆ interpolate_skewed_map_to_cubic()

void interpolate_skewed_map_to_cubic ( double **  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
unsigned  ,
unsigned  ,
unsigned  ,
int  ,
int  ,
int  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double   
)

Definition at line 1388 of file lib_vio.cpp.

1399 {
1400 
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;
1410 
1411  /* forward transform skewed -> rectangular (ax,ay,az,bx,by,bz,cx,cy,cz) */
1412  /* compute unit cell vectors a b c by intersection of projections in a,b plane; Bronstein & Semendjajew 2.6.6.1 */
1413  ax = 1.0;
1414  ay = 0.0;
1415  az = 0.0;
1416  bx = cos(PI * gamma / 180.0);
1417  by = sin(PI * gamma / 180.0);
1418  bz = 0.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) {
1429  error_divide_zero(71330, "lib_vio");
1430  }
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);
1436  if (scz < 0.0) {
1437  error_sqrt_negative(71340, "lib_vio");
1438  }
1439  cz = sqrt(scz);
1440 
1441  /* inverse transform rectangular -> skewed (aix,aiy,aiz,bix,biy,biz,cix,ciy,ciz) */
1442  aix = 1.0;
1443  aiy = 0.0;
1444  aiz = 0.0;
1445  bix = -bx / by;
1446  biy = 1.0 / by;
1447  biz = 0.0;
1448  cix = (bx * cy - cx * by) / (by * cz);
1449  ciy = -cy / (by * cz);
1450  ciz = 1.0 / cz;
1451 
1452  /* assign origin and map it to skewed coordinates if it is not yet skewed */
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;
1457  origx = gx;
1458  origy = gy;
1459  origz = gz;
1460  }
1461 
1462  /* compute actual x y z extent of the skewed map: */
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;
1490  uxmin = 1E20;
1491  for (i = 0; i < 8; i++) if (ux[i] < uxmin) uxmin = ux[i];
1492  uymin = 1E20;
1493  for (i = 0; i < 8; i++) if (uy[i] < uymin) uymin = uy[i];
1494  uzmin = 1E20;
1495  for (i = 0; i < 8; i++) if (uz[i] < uzmin) uzmin = uz[i];
1496  uxmax = -1E20;
1497  for (i = 0; i < 8; i++) if (ux[i] > uxmax) uxmax = ux[i];
1498  uymax = -1E20;
1499  for (i = 0; i < 8; i++) if (uy[i] > uymax) uymax = uy[i];
1500  uzmax = -1E20;
1501  for (i = 0; i < 8; i++) if (uz[i] > uzmax) uzmax = uz[i];
1502 
1503  /* for interpolation, find minimum voxel spacing, rounded to nearest 0.1 Angstrom for good measure */
1504  *pwidth = widthx; /* ax = 1 */
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;
1508 
1509  /* compute output origin */
1510  /* we start pphiout at or below lower bound of phiin, and assert the new origin is in register with origin of the orthogonal coordinate system */
1511  *porigx = *pwidth * floor(uxmin / *pwidth);
1512  *porigy = *pwidth * floor(uymin / *pwidth);
1513  *porigz = *pwidth * floor(uzmin / *pwidth);
1514 
1515  /* compute output map dimensions */
1516  /* we end pphiout at or above upper bound of phiin */
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;
1521 
1522  /* create output map, and loop through its cubic lattice to perform interpolation */
1523  do_vect(pphiout, pnvox);
1524  for (indz = 0; indz < (int)*pextz; indz++)
1525  for (indy = 0; indy < (int)*pexty; indy++)
1526  for (indx = 0; indx < (int)*pextx; indx++) {
1527  /* determine position in orthogonal coordinates */
1528  xpos = *porigx + indx * *pwidth;
1529  ypos = *porigy + indy * *pwidth;
1530  zpos = *porigz + indz * *pwidth;
1531 
1532  /* compute position of probe cube within skewed map in voxel units */
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);
1539  x1 = x0 + 1;
1540  y1 = y0 + 1;
1541  z1 = z0 + 1;
1542 
1543  /* if probe cube is fully within skewed map, do interpolate */
1544  if (x0 >= 0 && x1 < (int)nx && y0 >= 0 && y1 < (int)ny &&
1545  z0 >= 0 && z1 < (int)nz) {
1546  a = gx - x0;
1547  b = gy - y0;
1548  c = gz - z0;
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));
1558  }
1559  }
1560  printf("lib_vio> Conversion to cubic lattice completed. \n");
1561 }
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
c
doublereal * c
Definition: numerical_recipes.cpp:2230
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
gidz_general
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
Definition: lib_vwk.cpp:32
z0
#define z0
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
b
doublereal * b
Definition: numerical_recipes.cpp:2230
set_origin_get_mode
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)
Definition: lib_vio.cpp:1307
y0
#define y0
x0
#define x0
error_sqrt_negative
void error_sqrt_negative(int error_number, const char *program)
Definition: lib_err.cpp:222
error_divide_zero
void error_divide_zero(int error_number, const char *program)
Definition: lib_err.cpp:217
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32
PI
#define PI
Definition: tools.h:43
a
doublereal * a
Definition: numerical_recipes.cpp:2230

◆ permute_dimensions()

void permute_dimensions ( int  ,
unsigned  ,
unsigned  ,
unsigned  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
int  ,
int  ,
int  ,
int *  ,
int *  ,
int *   
)

Definition at line 1220 of file lib_vio.cpp.

1224 {
1225 
1226  switch (ordermode) {
1227  case 1:
1228  *nx = nc;
1229  *ny = nr;
1230  *nz = ns;
1231  *nxstart = ncstart;
1232  *nystart = nrstart;
1233  *nzstart = nsstart;
1234  break;
1235  case 2:
1236  *nx = nc;
1237  *ny = ns;
1238  *nz = nr;
1239  *nxstart = ncstart;
1240  *nystart = nsstart;
1241  *nzstart = nrstart;
1242  break;
1243  case 3:
1244  *nx = nr;
1245  *ny = nc;
1246  *nz = ns;
1247  *nxstart = nrstart;
1248  *nystart = ncstart;
1249  *nzstart = nsstart;
1250  break;
1251  case 4:
1252  *nx = ns;
1253  *ny = nc;
1254  *nz = nr;
1255  *nxstart = nsstart;
1256  *nystart = ncstart;
1257  *nzstart = nrstart;
1258  break;
1259  case 5:
1260  *nx = nr;
1261  *ny = ns;
1262  *nz = nc;
1263  *nxstart = nrstart;
1264  *nystart = nsstart;
1265  *nzstart = ncstart;
1266  break;
1267  case 6:
1268  *nx = ns;
1269  *ny = nr;
1270  *nz = nc;
1271  *nxstart = nsstart;
1272  *nystart = nrstart;
1273  *nzstart = ncstart;
1274  break;
1275  default:
1276  error_option(70212, "lib_vio");
1277  }
1278 }
error_option
void error_option(int error_number, const char *program)
Definition: lib_err.cpp:61

◆ permute_map()

void permute_map ( int  ,
unsigned  ,
unsigned  ,
unsigned  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
int  ,
int  ,
int  ,
int *  ,
int *  ,
int *  ,
double *  ,
double **   
)

Definition at line 1198 of file lib_vio.cpp.

1202 {
1203 
1204  unsigned long nvox, count;
1205 
1206  /* create permuted map pphi */
1207  nvox = nc * nr * ns;
1208  do_vect(pphi, nvox);
1209  for (count = 0; count < nvox; count++)
1210  *(*pphi + permuted_index(ordermode, count, nc, nr, ns)) = *(phi + count);
1211  free_vect_and_zero_ptr((void**)&phi);
1212 
1213  /* permute the map dimensions */
1214  permute_dimensions(ordermode, nc, nr, ns, nx, ny, nz, ncstart, nrstart,
1215  nsstart, nxstart, nystart, nzstart);
1216 }
free_vect_and_zero_ptr
void free_vect_and_zero_ptr(void **pvect)
Definition: lib_vec.cpp:83
permuted_index
unsigned long permuted_index(int ordermode, unsigned long count, unsigned nc, unsigned nr, unsigned ns)
Definition: lib_vio.cpp:1165
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32
permute_dimensions
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)
Definition: lib_vio.cpp:1220

◆ permute_print_info()

void permute_print_info ( int  ,
unsigned  ,
unsigned  ,
unsigned  ,
unsigned  ,
unsigned  ,
unsigned  ,
int  ,
int  ,
int  ,
int  ,
int  ,
int   
)

Definition at line 1282 of file lib_vio.cpp.

1283 {
1284 
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);
1300  } else {
1301  printf("lib_vio> C,R,S = X,Y,Z (no axis permutation). \n");
1302  }
1303 }

◆ permuted_index()

unsigned long permuted_index ( int  ,
unsigned  long,
unsigned  ,
unsigned  ,
unsigned   
)

Definition at line 1165 of file lib_vio.cpp.

1167 {
1168  unsigned ic, ir, is;
1169  unsigned long ncr, q;
1170 
1171  ncr = nc * nr;
1172  is = count / ncr;
1173  q = count - is * ncr;
1174  ir = q / nc;
1175  ic = q - ir * nc;
1176 
1177  switch (ordermode) {
1178  case 1:
1179  return ic + ir * nc + is * nc * nr;
1180  case 2:
1181  return ic + is * nc + ir * nc * ns;
1182  case 3:
1183  return ir + ic * nr + is * nr * nc;
1184  case 4:
1185  return is + ic * ns + ir * ns * nc;
1186  case 5:
1187  return ir + is * nr + ic * nr * ns;
1188  case 6:
1189  return is + ir * ns + ic * ns * nr;
1190  default:
1191  error_option(70210, "lib_vio");
1192  return 0;
1193  }
1194 }
error_option
void error_option(int error_number, const char *program)
Definition: lib_err.cpp:61
ir
int ir
Definition: image_find_center.cpp:57

◆ read_ascii()

void read_ascii ( char *  ,
unsigned  long,
double **   
)

Definition at line 163 of file lib_vio.cpp.

164 {
165  unsigned long count;
166  FILE *fin;
167  float currfloat;
168 
169  fin = fopen(vol_file, "r");
170  if (fin == NULL) {
171  error_open_filename(70310, "lib_vio", vol_file);
172  }
173  printf("lib_vio> Reading ASCII data... \n");
174  do_vect(fphi, nvox);
175  for (count = 0; count < nvox; count++) {
176  if (fscanf(fin, "%e", &currfloat) != 1) {
177  error_unreadable_file_short(70330 , "lib_vio", vol_file);
178  } else {
179  *(*fphi + count) = currfloat;
180  }
181  }
182  if (fscanf(fin, "%e", &currfloat) != EOF) {
183  error_unreadable_file_long(70340 , "lib_vio", vol_file);
184  }
185  printf("lib_vio> Volumetric data read from file %s\n", vol_file);
186  fclose(fin);
187 }
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
error_unreadable_file_long
void error_unreadable_file_long(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:148
error_unreadable_file_short
void error_unreadable_file_short(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:143
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32

◆ read_char()

int read_char ( char *  ,
FILE *   
)

Definition at line 874 of file lib_vio.cpp.

875 {
876 
877  if (fread(currchar, 1, 1, fin) != 1) return 0;
878  return 1;
879 }

◆ read_char_float()

int read_char_float ( float *  ,
FILE *   
)

Definition at line 863 of file lib_vio.cpp.

864 {
865  char currchar;
866 
867  if (fread(&currchar, 1, 1, fin) != 1) return 0;
868  *currfloat = (float)currchar;
869  return 1;
870 }

◆ read_float()

int read_float ( float *  ,
FILE *  ,
int   
)

Definition at line 817 of file lib_vio.cpp.

818 {
819  unsigned char *cptr, tmp;
820 
821  if (fread(currfloat, 4, 1, fin) != 1) return 0;
822  if (swap == 1) {
823  cptr = (unsigned char *)currfloat;
824  tmp = cptr[0];
825  cptr[0] = cptr[3];
826  cptr[3] = tmp;
827  tmp = cptr[1];
828  cptr[1] = cptr[2];
829  cptr[2] = tmp;
830  }
831  return 1;
832 }

◆ read_float_empty()

int read_float_empty ( FILE *  )

Definition at line 854 of file lib_vio.cpp.

855 {
856  float currfloat;
857 
858  if (fread(&currfloat, 4, 1, fin) != 1) return 0;
859  return 1;
860 }

◆ read_int()

int read_int ( int *  ,
FILE *  ,
int   
)

Definition at line 883 of file lib_vio.cpp.

884 {
885  unsigned char *cptr, tmp;
886 
887  if (fread(currlong, 4, 1, fin) != 1) return 0;
888  if (swap == 1) {
889  cptr = (unsigned char *)currlong;
890  tmp = cptr[0];
891  cptr[0] = cptr[3];
892  cptr[3] = tmp;
893  tmp = cptr[1];
894  cptr[1] = cptr[2];
895  cptr[2] = tmp;
896  }
897  return 1;
898 }

◆ read_mrc()

void read_mrc ( char *  ,
int *  ,
int *  ,
int *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
int *  ,
int *  ,
int *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double **   
)

Definition at line 369 of file lib_vio.cpp.

376 {
377 
378  unsigned long count, nvox;
379  FILE *fin;
380  int nc_tmp, nr_tmp, ns_tmp, mx, my, mz;
381  int mode;
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;
395 
396  n_range_viol0 = test_mrc(vol_file, 0);
397  n_range_viol1 = test_mrc(vol_file, 1);
398 
399  if (n_range_viol0 < n_range_viol1) { /* guess endianism */
400  swap = 0;
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);
403  }
404  } else {
405  swap = 1;
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);
408  }
409  }
410 
411  /* read header */
412  fin = fopen(vol_file, "rb");
413  if (fin == NULL) {
414  error_open_filename(70620, "lib_vio", vol_file);
415  }
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);
420  *nc = nc_tmp;
421  *nr = nr_tmp;
422  *ns = ns_tmp;
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);
430  header_ok *= read_float(&xlen, fin, swap);
431  header_ok *= read_float(&ylen, fin, swap);
432  header_ok *= read_float(&zlen, fin, swap);
433  header_ok *= read_float(&a_tmp, fin, swap);
434  header_ok *= read_float(&b_tmp, fin, swap);
435  header_ok *= read_float(&g_tmp, fin, swap);
436  *alpha = a_tmp;
437  *beta = b_tmp;
438  *gamma = g_tmp;
439  header_ok *= read_int(&mapc, fin, swap);
440  header_ok *= read_int(&mapr, fin, swap);
441  header_ok *= read_int(&maps, fin, swap);
442  header_ok *= read_float(&dmin, fin, swap);
443  header_ok *= read_float(&dmax, fin, swap);
444  header_ok *= read_float(&dmean, 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);
457  header_ok *= read_float(&skwtrn1, fin, swap);
458  header_ok *= read_float(&skwtrn2, fin, swap);
459  header_ok *= read_float(&skwtrn3, fin, swap);
460  for (i = 38; i < 50; ++i) header_ok *= read_float(&dummy, fin, swap);
461  header_ok *= read_float(&x_tmp, fin, swap);
462  header_ok *= read_float(&y_tmp, fin, swap);
463  header_ok *= read_float(&z_tmp, fin, swap);
464  *xorigin = x_tmp;
465  *yorigin = y_tmp;
466  *zorigin = z_tmp;
467  header_ok *= read_char(&mapchar1, fin);
468  header_ok *= read_char(&mapchar2, fin);
469  header_ok *= read_char(&mapchar3, fin);
470  header_ok *= read_char(&mapchar4, fin);
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);
475  header_ok *= read_float(&rms, fin, swap);
476 
477  if (header_ok == 0) {
478  error_file_header(70650, "lib_vio", vol_file);
479  }
480 
481  /* print some info */
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);
507  if (lskflg != 0) {
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");
523  }
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");
533  }
534 
535  /* extract data based on file type mode, currently supports 8-bit, 16-bit, and 32-bit */
536  /* note: fphi will contain data in the order of the input file, no axis permutation yet */
537  nvox = *nc * *nr * *ns;
538  switch (mode) {
539  case 0: /* char - converted to float, testing for signed-ness*/
540  rewind(fin);
541  do_vect(fphi, nvox);
542  for (count = 0; count < 256; ++count) if (read_float_empty(fin) == 0) error_file_convert(70642, "lib_vio", vol_file);
543  for (count = 0; count < (unsigned long)nsymbt; ++count) if (read_char_float(&currfloat, fin) == 0) error_file_convert(70643, "lib_vio", vol_file);
544  for (count = 0; count < nvox; ++count) {
545  if (read_char_float(&currfloat, fin) == 0) error_file_convert(70651, "lib_vio", vol_file);
546  else {
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;
554  }
555  }
556  fclose(fin);
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");
561  }
562  break;
563  case 1: /* 16-bit float */
564  rewind(fin);
565  do_vect(fphi, nvox);
566  for (count = 0; count < 256; ++count) if (read_float_empty(fin) == 0) {
567  error_file_convert(70644, "lib_vio", vol_file);
568  }
569  for (count = 0; count < (unsigned long)nsymbt; ++count) if (read_char_float(&currfloat, fin) == 0) {
570  error_file_convert(70645, "lib_vio", vol_file);
571  }
572 
573  for (count = 0; count < nvox; ++count) {
574  if (read_short_float(&currfloat, fin, swap) == 0) {
575  error_file_convert(70652, "lib_vio", vol_file);
576  } else {
577  *(*fphi + count) = currfloat;
578  }
579  }
580  fclose(fin);
581  break;
582  case 2: /* 32-bit float */
583  rewind(fin);
584  do_vect(fphi, nvox);
585  for (count = 0; count < 256; ++count) if (read_float_empty(fin) == 0) {
586  error_file_convert(70646, "lib_vio", vol_file);
587  }
588  for (count = 0; count < (unsigned long)nsymbt; ++count) if (read_char_float(&currfloat, fin) == 0) {
589  error_file_convert(70647, "lib_vio", vol_file);
590  }
591  for (count = 0; count < nvox; ++count) {
592  if (read_float(&currfloat, fin, swap) == 0) {
593  error_file_convert(70653, "lib_vio", vol_file);
594  } else {
595  *(*fphi + count) = currfloat;
596  }
597  }
598  fclose(fin);
599  break;
600  default:
601  error_file_float_mode(70654, "lib_vio", vol_file);
602  }
603 
604  /* assign voxel spacing parameters */
605  *widthx = xlen / (double) mx;
606  *widthy = ylen / (double) my;
607  *widthz = zlen / (double) mz;
608 
609  /* test for orthogonal and cubic lattice */
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;
615 
616  /* set axis permutation mode */
617  *ordermode = 7;
618  if (mapc == 1 && mapr == 2 && maps == 3) {
619  *ordermode = 1;
620  }
621  if (mapc == 1 && mapr == 3 && maps == 2) {
622  *ordermode = 2;
623  }
624  if (mapc == 2 && mapr == 1 && maps == 3) {
625  *ordermode = 3;
626  }
627  if (mapc == 2 && mapr == 3 && maps == 1) {
628  *ordermode = 4;
629  }
630  if (mapc == 3 && mapr == 1 && maps == 2) {
631  *ordermode = 5;
632  }
633  if (mapc == 3 && mapr == 2 && maps == 1) {
634  *ordermode = 6;
635  }
636  if (*ordermode == 7) {
637  error_axis_assignment(70680, "lib_vio");
638  }
639  printf("lib_vio> Volumetric data read from file %s\n", vol_file);
640 }
rms
double rms(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=nullptr, MultidimArray< double > *Contributions=nullptr)
Definition: filters.h:725
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
read_float
int read_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:817
read_int
int read_int(int *currlong, FILE *fin, int swap)
Definition: lib_vio.cpp:883
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
read_short_float
int read_short_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:836
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
dmin
#define dmin(a, b)
Definition: numerical_recipes.cpp:2112
read_char
int read_char(char *currchar, FILE *fin)
Definition: lib_vio.cpp:874
mode
void mode
Definition: numerical_recipes.cpp:4411
dummy
double dummy
Definition: numerical_recipes.cpp:4635
error_axis_assignment
void error_axis_assignment(int error_number, const char *program)
Definition: lib_err.cpp:198
test_mrc
int test_mrc(char *vol_file, int swap)
Definition: lib_vio.cpp:1010
dmax
#define dmax(a, b)
Definition: numerical_recipes.cpp:2113
error_file_convert
void error_file_convert(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:183
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32
read_float_empty
int read_float_empty(FILE *fin)
Definition: lib_vio.cpp:854
read_char_float
int read_char_float(float *currfloat, FILE *fin)
Definition: lib_vio.cpp:863
error_file_header
void error_file_header(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:188
error_file_float_mode
void error_file_float_mode(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:193

◆ read_short_float()

int read_short_float ( float *  ,
FILE *  ,
int   
)

Definition at line 836 of file lib_vio.cpp.

837 {
838  unsigned char *cptr, tmp;
839  short currshort;
840 
841  if (fread(&currshort, 2, 1, fin) != 1) return 0;
842  if (swap == 1) {
843  cptr = (unsigned char *)&currshort;
844  tmp = cptr[0];
845  cptr[0] = cptr[1];
846  cptr[1] = tmp;
847  }
848  *currfloat = (float)currshort;
849  return 1;
850 }

◆ read_situs()

void read_situs ( char *  ,
double *  ,
double *  ,
double *  ,
double *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double **   
)

Definition at line 80 of file lib_vio.cpp.

83 {
84  unsigned long nvox, count;
85  double dorigx, dorigy, dorigz, dwidth;
86  double phitest, dtemp;
87  const char *program = "lib_vio";
88  FILE *fin;
89 
90  fin = fopen(vol_file, "r");
91  if (fin == NULL) {
92  error_open_filename(13010, program, vol_file);
93  }
94 
95  /* read header and print information */
96  if (7 != fscanf(fin, "%le %le %le %le %d %d %d", &dwidth, &dorigx, &dorigy, &dorigz, extx, exty, extz)) error_fscanf(program, vol_file);
97  *width = dwidth;
98  *origx = dorigx;
99  *origy = dorigy;
100  *origz = dorigz;
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);
105 
106  nvox = *extx * *exty * *extz;
107 
108  /* allocate memory and read data */
109  printf("lib_vio> Reading density data... \n");
110  *phi = (double *) alloc_vect(nvox, sizeof(double));
111 
112  for (count = 0; count < nvox; count++) {
113  if (fscanf(fin, "%le", &dtemp) != 1) {
114  error_unreadable_file_long(13030, program, vol_file);
115  } else *(*phi + count) = dtemp;
116  }
117  if (fscanf(fin, "%le", &phitest) != EOF) {
118  error_unreadable_file_long(13040, program, vol_file);
119  }
120  fclose(fin);
121  printf("lib_vio> Volumetric data read from file %s\n", vol_file);
122  return;
123 }
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
error_fscanf
void error_fscanf(const char *program, const char *file)
Definition: lib_err.cpp:46
error_unreadable_file_long
void error_unreadable_file_long(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:148
alloc_vect
void * alloc_vect(unsigned int n, size_t elem_size)
Definition: lib_vec.cpp:71

◆ read_spider()

void read_spider ( char *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double **   
)

Definition at line 644 of file lib_vio.cpp.

645 {
646  unsigned long nvox;
647  FILE *fin;
648  int i, swap, header_ok = 1, headlen;
649  unsigned long count;
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;
655 
656  n_range_viol0 = test_spider(vol_file, 0);
657  n_range_viol1 = test_spider(vol_file, 1);
658 
659  if (n_range_viol0 < n_range_viol1) { /* guess endianism */
660  swap = 0;
661  if (n_range_viol0 > 0) {
662  printf("lib_vio> Warning: %i header field range violations detected \n", n_range_viol0);
663  }
664  } else {
665  swap = 1;
666  if (n_range_viol1 > 0) {
667  printf("lib_vio> Warning: %i header field range violations detected \n", n_range_viol1);
668  }
669  }
670 
671  /* read header */
672  fin = fopen(vol_file, "rb");
673  if (fin == NULL) {
674  error_open_filename(70820, "lib_vio", vol_file);
675  }
676  printf("lib_vio> Reading header information from SPIDER file %s \n", vol_file);
677  header_ok *= read_float(&nslice, fin, swap);
678  header_ok *= read_float(&nrow, fin, swap);
679  header_ok *= read_float(&dummy, fin, swap);
680  header_ok *= read_float(&dummy, fin, swap);
681  header_ok *= read_float(&iform, fin, swap);
682  header_ok *= read_float(&imami, fin, swap);
683  header_ok *= read_float(&fmax, fin, swap);
684  header_ok *= read_float(&fmin, fin, swap);
685  header_ok *= read_float(&av, fin, swap);
686  header_ok *= read_float(&sig, fin, swap);
687  header_ok *= read_float(&dummy, fin, swap);
688  header_ok *= read_float(&nsam, fin, swap);
689  header_ok *= read_float(&headrec, fin, swap);
690  header_ok *= read_float(&iangle, fin, swap);
691  header_ok *= read_float(&phi, fin, swap);
692  header_ok *= read_float(&theta, fin, swap);
693  header_ok *= read_float(&gamma, fin, swap);
694  header_ok *= read_float(&xoff, fin, swap);
695  header_ok *= read_float(&yoff, fin, swap);
696  header_ok *= read_float(&zoff, fin, swap);
697  header_ok *= read_float(&scale, fin, swap);
698  header_ok *= read_float(&labbyt, fin, swap);
699  header_ok *= read_float(&lenbyt, fin, swap);
700  header_ok *= read_float(&istack, fin, swap);
701  header_ok *= read_float(&inuse, fin, swap);
702  header_ok *= read_float(&maxim, fin, swap);
703  for (i = 0; i < 4; ++i) header_ok *= read_float(&dummy, fin, swap);
704  header_ok *= read_float(&kangle, fin, swap);
705  header_ok *= read_float(&phi1, fin, swap);
706  header_ok *= read_float(&theta1, fin, swap);
707  header_ok *= read_float(&psi1, fin, swap);
708  header_ok *= read_float(&phi2, fin, swap);
709  header_ok *= read_float(&theta2, fin, swap);
710  header_ok *= read_float(&psi2, fin, swap);
711  if (header_ok == 0) {
712  error_file_header(70850, "lib_vio", vol_file);
713  }
714 
715  /* print some info */
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);
746 
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));
752  do_vect(fphi, nvox);
753  rewind(fin);
754  for (count = 0; count < (unsigned long)headlen; ++count) if (read_float_empty(fin) == 0) {
755  error_file_convert(70841, "lib_vio", vol_file);
756  }
757  for (count = 0; count < nvox; ++count) if (read_float(&currfloat, fin, swap) == 0) {
758  error_file_convert(70842, "lib_vio", vol_file);
759  } else {
760  *(*fphi + count) = currfloat;
761  }
762  fclose(fin);
763 
764  /* notes and error checks */
765  if (((int) floor(100 * (sizeof(float)*headlen - labbyt) + 0.5)) != 0) {
766  error_spider_header(70860, "lib_vio");
767  }
768  printf("lib_vio> Volumetric data read from file %s\n", vol_file);
769 
770 }
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
read_float
int read_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:817
test_spider
int test_spider(char *vol_file, int swap)
Definition: lib_vio.cpp:1089
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
theta
double theta
Definition: numerical_recipes.cpp:6214
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
dummy
double dummy
Definition: numerical_recipes.cpp:4635
error_file_convert
void error_file_convert(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:183
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32
read_float_empty
int read_float_empty(FILE *fin)
Definition: lib_vio.cpp:854
error_spider_header
void error_spider_header(int error_number, const char *program)
Definition: lib_err.cpp:207
fmax
double fmax
Definition: numerical_recipes.cpp:5166
error_file_header
void error_file_header(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:188

◆ read_vol()

void read_vol ( char *  ,
double *  ,
double *  ,
double *  ,
double *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double **   
)

Definition at line 25 of file lib_vio.cpp.

28 {
29  int n_range_viol0, n_range_viol1, n_range_viol2;
30  int ordermode = 7, cubic = 1, orom = 1;
31  unsigned nc, nr, ns;
32  unsigned nx, ny, nz;
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;
37  double alpha, beta, gamma;
38  double *phi_raw;
39 
40  n_range_viol0 = test_mrc(vol_file, 0);
41  n_range_viol1 = test_mrc(vol_file, 1);
42  n_range_viol2 = test_situs(vol_file);
43 
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);
46  else {
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);
52  permute_print_info(ordermode, nc, nr, ns, nx, ny, nz, ncstart, nrstart,
53  nsstart, nxstart, nystart, nzstart);
54  assert_cubic_map(orom, cubic, alpha, beta, gamma, widthx, widthy, widthz,
55  nx, ny, nz, nxstart, nystart, nzstart, xorigin, yorigin,
56  zorigin, width, origx, origy, origz, extx, exty, extz, phi);
57  printf("lib_vio> \n");
58  }
59 
60  return;
61 }
assert_cubic_map
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)
Definition: lib_vio.cpp:1326
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
permute_map
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)
Definition: lib_vio.cpp:1198
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
test_mrc
int test_mrc(char *vol_file, int swap)
Definition: lib_vio.cpp:1010
read_situs
void read_situs(char *vol_file, double *width, double *origx, double *origy, double *origz, unsigned *extx, unsigned *exty, unsigned *extz, double **phi)
Definition: lib_vio.cpp:80
read_mrc
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)
Definition: lib_vio.cpp:369
permute_print_info
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)
Definition: lib_vio.cpp:1282
test_situs
int test_situs(char *vol_file)
Definition: lib_vio.cpp:937

◆ read_xplor()

void read_xplor ( char *  ,
int *  ,
int *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int *  ,
int *  ,
int *  ,
unsigned *  ,
unsigned *  ,
unsigned *  ,
double **   
)

Definition at line 191 of file lib_vio.cpp.

196 {
197 
198  unsigned long nvox;
199  FILE *fin;
200  int idummy, done;
201  char *nextline;
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;
208 
209  nextline = (char *) alloc_vect(FLENGTH, sizeof(char));
210 
211  fin = fopen(vol_file, "r");
212  if (fin == NULL) {
213  error_open_filename(70420, "lib_vio", vol_file);
214  }
215 
216  /* ignore header length line */
217  xplor_skip_to_number(&fin, &nextline);
218 
219  /* read index line */
220  xplor_skip_to_number(&fin, &nextline);
221  if (sscanf(nextline, "%8ld%8ld%8ld%8ld%8ld%8ld%8ld%8ld%8ld", &mx, &mxstart, &mxend, &my, &mystart, &myend, &mz, &mzstart, &mzend) != 9) {
222  error_xplor_file_indexing(70430, "lib_vio");
223  }
224  *extx = mxend - mxstart + 1;
225  *exty = myend - mystart + 1;
226  *extz = mzend - mzstart + 1;
227  nvox = *extx * *exty * *extz;
228 
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);
239 
240  /* read unit cell info and determine grid width and origin */
241  xplor_skip_to_number(&fin, &nextline);
242  if (sscanf(nextline, "%12f%12f%12f%12f%12f%12f", &xlen, &ylen, &zlen, &a_tmp, &b_tmp, &g_tmp) != 6) {
243  error_xplor_file_unit_cell(70440, "lib_vio");
244  } else {
245  *alpha = a_tmp;
246  *beta = b_tmp;
247  *gamma = g_tmp;
248  }
249 
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);
257 
258  /* assign voxel spacing parameters */
259  *widthx = xlen / (double) mx;
260  *widthy = ylen / (double) my;
261  *widthz = zlen / (double) mz;
262  *nxstart = mxstart;
263  *nystart = mystart;
264  *nzstart = mzstart;
265 
266  /* test for orthogonal and cubic lattice */
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)
274  *cubic = 0;
275 
276  /* read ZYX info */
277  for (done = 0; done == 0;) { /* read next line and check if it contains ZYX */
278  if (fgets(nextline, FLENGTH, fin) == NULL) {
279  error_EOF_ZYX_mode(70450, "lib_vio");
280  }
281  if (*nextline == 'Z' && *(nextline + 1) == 'Y' && *(nextline + 2) == 'X') {
282  done = 1;
283  break;
284  }
285  if (*nextline == 'z' && *(nextline + 1) == 'y' && *(nextline + 2) == 'x') {
286  done = 1;
287  break;
288  }
289  }
290 
291  /* read sections */
292  do_vect(fphi, nvox);
293  for (indz = 0; indz < (int)*extz; indz++) {
294  /* skip section header and read section number */
295  xplor_skip_to_number(&fin, &nextline);
296  if (sscanf(nextline, "%8d", &idummy) != 1) {
297  error_xplor_file_map_section(70470, "lib_vio");
298  }
299  if (idummy != indz) {
300  error_xplor_file_map_section_number(70480, "lib_vio");
301  }
302 
303  /* read section data */
304  for (indy = 0; indy < (int)*exty; indy++) for (indx = 0; indx < (int)*extx; indx++) {
305  if (fscanf(fin, "%12f", &currfloat) != 1) {
306  error_unreadable_file_short(70490, "lib_vio", vol_file);
307  } else {
308  *(*fphi + gidz_general(indz, indy, indx, *exty, *extx)) = currfloat;
309  }
310  }
311  }
312 
313  /* read end of data marker */
314  xplor_skip_to_number(&fin, &nextline);
315  if (sscanf(nextline, "%8d", &idummy) != 1 || idummy != -9999) {
316  error_xplor_maker("lib_vio");
317  }
318 
319  fclose(fin);
320  free_vect_and_zero_ptr((void**)&nextline);
321  printf("lib_vio> Volumetric data read from file %s\n", vol_file);
322 }
xplor_skip_to_number
void xplor_skip_to_number(FILE **fin, char **nextline)
Definition: lib_vio.cpp:327
error_xplor_file_map_section
void error_xplor_file_map_section(int error_number, const char *program)
Definition: lib_err.cpp:163
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
free_vect_and_zero_ptr
void free_vect_and_zero_ptr(void **pvect)
Definition: lib_vec.cpp:83
error_xplor_file_indexing
void error_xplor_file_indexing(int error_number, const char *program)
Definition: lib_err.cpp:153
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
error_xplor_file_map_section_number
void error_xplor_file_map_section_number(int error_number, const char *program)
Definition: lib_err.cpp:168
gidz_general
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
Definition: lib_vwk.cpp:32
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
done
done
Definition: numerical_recipes.cpp:7636
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
FLENGTH
#define FLENGTH
Definition: lib_vio.cpp:22
error_EOF_ZYX_mode
void error_EOF_ZYX_mode(int error_number, const char *program)
Definition: lib_err.cpp:173
error_unreadable_file_short
void error_unreadable_file_short(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:143
error_xplor_maker
void error_xplor_maker(const char *program)
Definition: lib_err.cpp:179
do_vect
void do_vect(double **pvect, unsigned long len)
Definition: lib_vec.cpp:32
alloc_vect
void * alloc_vect(unsigned int n, size_t elem_size)
Definition: lib_vec.cpp:71
error_xplor_file_unit_cell
void error_xplor_file_unit_cell(int error_number, const char *program)
Definition: lib_err.cpp:158

◆ set_origin_get_mode()

int set_origin_get_mode ( int  ,
int  ,
int  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double  ,
double *  ,
double *  ,
double *   
)

Definition at line 1307 of file lib_vio.cpp.

1308 {
1309 
1310  if (fabs(xorigin) < 0.0001 && fabs(yorigin) < 0.0001 && fabs(zorigin) < 0.0001) { /* seem to have crystallographic origin */
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;
1315  return 0;
1316  } else { /* seem to have MRC2000 */
1317  printf("lib_vio> Using MRC2000 style origin defined by [X,Y,Z]ORIGIN fields.\n");
1318  *origx = xorigin;
1319  *origy = yorigin;
1320  *origz = zorigin;
1321  return 1;
1322  }
1323 }

◆ test_mrc()

int test_mrc ( char *  ,
int   
)

Definition at line 1010 of file lib_vio.cpp.

1011 {
1012  FILE *fin;
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;
1018  float alpha, beta, gamma;
1019  float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin;
1020 
1021  fin = fopen(vol_file, "rb");
1022  if (fin == NULL) {
1023  error_open_filename(71010, "lib_vio", vol_file);
1024  }
1025 
1026  /* read header info */
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);
1037  header_ok *= read_float(&xlen, fin, swap);
1038  header_ok *= read_float(&ylen, fin, swap);
1039  header_ok *= read_float(&zlen, fin, swap);
1040  header_ok *= read_float(&alpha, fin, swap);
1041  header_ok *= read_float(&beta, fin, swap);
1042  header_ok *= read_float(&gamma, 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);
1046  header_ok *= read_float(&dmin, fin, swap);
1047  header_ok *= read_float(&dmax, fin, swap);
1048  header_ok *= read_float(&dmean, 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);
1053  fclose(fin);
1054  if (header_ok == 0) {
1055  error_file_header(71020, "lib_vio", vol_file);
1056  }
1057 
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);
1082 
1083  return n_range_viols;
1084 }
read_float
int read_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:817
read_int
int read_int(int *currlong, FILE *fin, int swap)
Definition: lib_vio.cpp:883
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
dmin
#define dmin(a, b)
Definition: numerical_recipes.cpp:2112
mode
void mode
Definition: numerical_recipes.cpp:4411
dummy
double dummy
Definition: numerical_recipes.cpp:4635
dmax
#define dmax(a, b)
Definition: numerical_recipes.cpp:2113
error_file_header
void error_file_header(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:188

◆ test_registration()

int test_registration ( float  ,
float  ,
float  ,
float   
)

Definition at line 917 of file lib_vio.cpp.

918 {
919  float xreg, xreg1, yreg, yreg1, zreg, zreg1;
920 
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;
931  else return 1;
932 }

◆ test_situs()

int test_situs ( char *  )

Definition at line 937 of file lib_vio.cpp.

938 {
939 
940  FILE *fin;
941  double width, origx, origy, origz;
942  unsigned extx, exty, extz;
943  int header_ok = 1, n_range_viols = 0;
944 
945  fin = fopen(vol_file, "r");
946  if (fin == NULL) {
947  error_open_filename(71010, "lib_vio", vol_file);
948  }
949 
950  /* read header */
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;
958  fclose(fin);
959 
960  if (header_ok == 0) return 999999; /* worst case, can't read data */
961  else {
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;
977  }
978 }
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67

◆ test_situs_header_and_suffix()

int test_situs_header_and_suffix ( char *  )

Definition at line 1001 of file lib_vio.cpp.

1002 {
1003 
1004  if (have_situs_suffix(vol_file) == 0) return 999999;
1005  else return test_situs(vol_file);
1006 }
have_situs_suffix
int have_situs_suffix(char *vol_file)
Definition: lib_vio.cpp:981
test_situs
int test_situs(char *vol_file)
Definition: lib_vio.cpp:937

◆ test_spider()

int test_spider ( char *  ,
int   
)

Definition at line 1089 of file lib_vio.cpp.

1090 {
1091  FILE *fin;
1092  int i, header_ok = 1, n_range_viols = 0, headlen;
1093  float dummy;
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;
1097 
1098  fin = fopen(vol_file, "rb");
1099  if (fin == NULL) {
1100  error_open_filename(71210, "lib_vio", vol_file);
1101  }
1102 
1103  /* read header info */
1104  header_ok *= read_float(&nslice, fin, swap);
1105  header_ok *= read_float(&nrow, fin, swap);
1106  header_ok *= read_float(&dummy, fin, swap);
1107  header_ok *= read_float(&dummy, fin, swap);
1108  header_ok *= read_float(&iform, fin, swap);
1109  header_ok *= read_float(&imami, fin, swap);
1110  header_ok *= read_float(&fmax, fin, swap);
1111  header_ok *= read_float(&fmin, fin, swap);
1112  header_ok *= read_float(&av, fin, swap);
1113  header_ok *= read_float(&sig, fin, swap);
1114  header_ok *= read_float(&dummy, fin, swap);
1115  header_ok *= read_float(&nsam, fin, swap);
1116  header_ok *= read_float(&headrec, fin, swap);
1117  header_ok *= read_float(&iangle, fin, swap);
1118  header_ok *= read_float(&phi, fin, swap);
1119  header_ok *= read_float(&theta, fin, swap);
1120  header_ok *= read_float(&gamma, fin, swap);
1121  header_ok *= read_float(&xoff, fin, swap);
1122  header_ok *= read_float(&yoff, fin, swap);
1123  header_ok *= read_float(&zoff, fin, swap);
1124  header_ok *= read_float(&scale, fin, swap);
1125  header_ok *= read_float(&labbyt, fin, swap);
1126  header_ok *= read_float(&lenbyt, fin, swap);
1127  header_ok *= read_float(&istack, fin, swap);
1128  header_ok *= read_float(&inuse, fin, swap);
1129  header_ok *= read_float(&maxim, fin, swap);
1130  for (i = 0; i < 4; ++i) header_ok *= read_float(&dummy, fin, swap);
1131  header_ok *= read_float(&kangle, fin, swap);
1132  header_ok *= read_float(&phi1, fin, swap);
1133  header_ok *= read_float(&theta1, fin, swap);
1134  header_ok *= read_float(&psi1, fin, swap);
1135  header_ok *= read_float(&phi2, fin, swap);
1136  header_ok *= read_float(&theta2, fin, swap);
1137  header_ok *= read_float(&psi2, fin, swap);
1138  fclose(fin);
1139  if (header_ok == 0) {
1140  error_file_header(71220, "lib_vio", vol_file);
1141  }
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;
1157 }
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
read_float
int read_float(float *currfloat, FILE *fin, int swap)
Definition: lib_vio.cpp:817
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
theta
double theta
Definition: numerical_recipes.cpp:6214
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
dummy
double dummy
Definition: numerical_recipes.cpp:4635
fmax
double fmax
Definition: numerical_recipes.cpp:5166
error_file_header
void error_file_header(int error_number, const char *program, const char *filename)
Definition: lib_err.cpp:188

◆ write_mrc()

void write_mrc ( int  ,
char *  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
double *   
)

Definition at line 1646 of file lib_vio.cpp.

1649 {
1650  FILE *fout;
1651  long nxstart, nystart, nzstart;
1652  float xorigin, yorigin, zorigin;
1653  unsigned long count, pnvox;
1654  long nx, ny, nz, mx, my, mz;
1655  long mode = 2;
1656  int modesel = 1;
1657  float xlen, ylen, zlen, alpha, beta, gamma;
1658  long mapc, mapr, maps;
1659  long ispg = 1, nlabl = 0;
1660  float dummy = 0.0f;
1661  double dmax, dmin, dav, dsig;
1662  float fmax, fmin, fav, fsig;
1663  double dshift, dscale;
1664  char dummychar = '\0';
1665  char mapstring[4];
1666  char machstring[4];
1667  int i, wmaperr;
1668  long lskflg = 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};
1673 
1674  fout = fopen(vol_file, "wb");
1675  if (fout == NULL) {
1676  error_open_filename(71420, "lib_vio", vol_file);
1677  }
1678 
1679  /* set initial values of map parameters */
1680  pnvox = pextx * pexty * pextz;
1681  nx = pextx;
1682  ny = pexty;
1683  nz = pextz;
1684  mapstring[0] = 'M';
1685  mapstring[1] = 'A';
1686  mapstring[2] = 'P';
1687  mapstring[3] = ' ';
1688  if (*((int *)endchar) > 65536) { /* big endian */
1689  machstring[0] = 17;
1690  machstring[1] = 17;
1691  } else { /* little endian */
1692  machstring[0] = 68;
1693  machstring[1] = 65;
1694  }
1695  machstring[2] = 0;
1696  machstring[3] = 0;
1697  calc_map_info(pphi, pnvox, &dmax, &dmin, &dav, &dsig);
1698  fmax = dmax;
1699  fmin = dmin;
1700  fav = dav;
1701  fsig = dsig;
1702  xorigin = porigx;
1703  yorigin = porigy;
1704  zorigin = porigz;
1705  mx = pextx;
1706  my = pexty;
1707  mz = pextz;
1708  xlen = pwidth * mx;
1709  ylen = pwidth * my;
1710  zlen = pwidth * mz;
1711  alpha = 90.0f;
1712  beta = 90.0f;
1713  gamma = 90.0f;
1714  mapc = 1;
1715  mapr = 2;
1716  maps = 3;
1717 
1718  /* if grid is not in register with origin, we set the CCP4 start indices to zero */
1719  if (test_registration(porigx, porigy, porigz, pwidth) == 0) {
1720  nxstart = 0;
1721  nystart = 0;
1722  nzstart = 0;
1723  } else {
1724  nxstart = (long)floor((porigx / pwidth) + 0.5);
1725  nystart = (long)floor((porigy / pwidth) + 0.5);
1726  nzstart = (long)floor((porigz / pwidth) + 0.5);
1727  }
1728 
1729  /* optional manual override of variable map parameters */
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: ");
1744  modesel = readln_int();
1745  if (modesel < 1 && modesel > 7) {
1746  modesel = 1;
1747  mode = 2;
1748  printf("lib_vio> Did not recognize value, assuming MODE 2, 32-bit floats. \n");
1749  } else if (modesel > 1) mode = 0;
1750 
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: ");
1754  nx = readln_int();
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: ");
1757  ny = readln_int();
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: ");
1760  nz = readln_int();
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");
1765  }
1766  printf("lib_vio> The currently assigned CCP4-style NCSTART field (column start index) is %ld \n", nxstart);
1767  if (nxstart == 0 && test_registration(porigx, porigy, porigz, pwidth) == 0) {
1768  printf("lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1769  }
1770  printf("lib_vio> Enter the same or a new value: ");
1771  nxstart = readln_int();
1772  printf("lib_vio> The currently assigned CCP4-style NRSTART field (row start index) is %ld \n", nystart);
1773  if (nxstart == 0 && test_registration(porigx, porigy, porigz, pwidth) == 0) {
1774  printf("lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1775  }
1776  printf("lib_vio> Enter the same or a new value: ");
1777  nystart = readln_int();
1778  printf("lib_vio> The currently assigned CCP4-style NSSTART field (section start index) is %ld \n", nzstart);
1779  if (nxstart == 0 && test_registration(porigx, porigy, porigz, pwidth) == 0) {
1780  printf("lib_vio> Set to zero by default because input grid not in register with origin of coordinate system.\n");
1781  }
1782  printf("lib_vio> Enter the same or a new value: ");
1783  nzstart = readln_int();
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: ");
1786  mx = readln_int();
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: ");
1789  my = readln_int();
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: ");
1792  mz = readln_int();
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: ");
1795  xlen = readln_double();
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: ");
1798  ylen = readln_double();
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: ");
1801  zlen = readln_double();
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: ");
1804  alpha = readln_double();
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: ");
1807  beta = readln_double();
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: ");
1810  gamma = readln_double();
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: ");
1813  mapc = readln_int();
1814  wmaperr = 1;
1815  if (mapc == 1 || mapc == 2 || mapc == 3) wmaperr = 0;
1816  if (wmaperr) {
1817  mapc = 1;
1818  mapr = 2;
1819  maps = 3;
1820  printf("lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1821  } else {
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: ");
1824  mapr = readln_int();
1825  wmaperr = 1;
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;
1832  if (wmaperr) {
1833  mapc = 1;
1834  mapr = 2;
1835  maps = 3;
1836  printf("lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1837  } else {
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: ");
1840  maps = readln_int();
1841  wmaperr = 1;
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;
1848  if (wmaperr) {
1849  mapc = 1;
1850  mapr = 2;
1851  maps = 3;
1852  printf("lib_vio> Inconsistent axis order values, assuming mapc = 1, mapr = 2, maps = 3. \n");
1853  }
1854  }
1855  }
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: ");
1858  xorigin = readln_double();
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: ");
1861  yorigin = readln_double();
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: ");
1864  zorigin = readln_double();
1865  } /* end manual override */
1866 
1867  /* adjust MODE 0 maps as necessary */
1868  switch (modesel) {
1869  case 1:
1870  break;
1871  case 2:
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");
1874  else
1875  printf("lib_vio> Density values already fit signed 8-bit data format, no clipping necessary.\n");
1876  break;
1877  case 3:
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);
1883  normalize(pphi, pnvox, dscale);
1884  break;
1885  case 4:
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);
1891  floatshift(pphi, pnvox, dshift);
1892  normalize(pphi, pnvox, dscale);
1893  break;
1894  case 5:
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");
1897  break;
1898  case 6:
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);
1902  normalize(pphi, pnvox, 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");
1904  break;
1905  case 7:
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);
1911  floatshift(pphi, pnvox, dshift);
1912  normalize(pphi, pnvox, dscale);
1913  break;
1914  default:
1915  error_option(70211, "lib_vio");
1916  }
1917 
1918  /* recompute map statistics and then adjust unsigned MODE 0 values for storage as signed */
1919  if (modesel > 1) {
1920  calc_map_info(pphi, pnvox, &dmax, &dmin, &dav, &dsig);
1921  fmax = dmax;
1922  fmin = dmin;
1923  fav = dav;
1924  fsig = dsig;
1925  }
1926  if (modesel > 4) {
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;
1931  }
1932 
1933  /* Warning */
1934  /* Do not recompute map statistics after unsigned to signed conversion */
1935 
1936  /* write header */
1937  printf("lib_vio> Writing MRC / CCP4 (binary) volumetric map \n");
1938  wmaperr = 0;
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;
1985 
1986  /* write actual data */
1987  switch (mode) {
1988  case 0:
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;
1994  }
1995  break;
1996  case 2:
1997  for (count = 0; count < pnvox; count++) {
1998  dummy = *(pphi + count);
1999  if (fwrite(&dummy, 4, 1, fout) != 1) wmaperr = 1;
2000  }
2001  break;
2002  default:
2003  wmaperr = 1;
2004  break;
2005  }
2006  if (wmaperr != 0) {
2007  error_write_filename(71430, "lib_vio");
2008  }
2009  fclose(fout);
2010 
2011  /* print some info */
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);
2039  if (lskflg != 0) {
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);
2052  }
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);
2059 
2060  /* if auto mode and grid is not in register with origin, issue a warning */
2061  if (automode == 1 && test_registration(porigx, porigy, porigz, pwidth) == 0) {
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");
2067  }
2068 }
test_registration
int test_registration(float origx, float origy, float origz, float width)
Definition: lib_vio.cpp:917
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
error_option
void error_option(int error_number, const char *program)
Definition: lib_err.cpp:61
error_write_filename
void error_write_filename(int error_number, const char *program)
Definition: lib_err.cpp:227
alglib::beta
double beta(const double a, const double b)
Definition: specialfunctions.cpp:962
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
calc_map_info
void calc_map_info(double *phi, unsigned long nvox, double *maxdens, double *mindens, double *ave, double *sig)
Definition: lib_vwk.cpp:491
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
readln_double
double readln_double()
Definition: lib_std.cpp:38
readln_int
int readln_int()
Definition: lib_std.cpp:19
dmin
#define dmin(a, b)
Definition: numerical_recipes.cpp:2112
cif::normalize
quaternion_type< T > normalize(quaternion_type< T > q)
Definition: point.cpp:278
mode
void mode
Definition: numerical_recipes.cpp:4411
dummy
double dummy
Definition: numerical_recipes.cpp:4635
floatshift
void floatshift(double *phi, unsigned long nvox, double dens)
Definition: lib_vwk.cpp:594
dmax
#define dmax(a, b)
Definition: numerical_recipes.cpp:2113
clipped
int clipped(double *phi, unsigned long nvox, double max, double min)
Definition: lib_vwk.cpp:603
fprintf
fprintf(glob_prnt.io, "\)
fmax
double fmax
Definition: numerical_recipes.cpp:5166

◆ write_situs()

void write_situs ( char *  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
double *   
)

Definition at line 127 of file lib_vio.cpp.

130 {
131  unsigned long nvox, count;
132  const char *program = "lib_vio";
133  FILE *fout;
134 
135  nvox = extx * exty * extz;
136  fout = fopen(vol_file, "w");
137  if (fout == NULL) {
138  error_open_filename(13210, program, vol_file);
139  }
140 
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);
143  fprintf(fout, "\n");
144 
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));
148  }
149  fclose(fout);
150  printf("lib_vio> Volumetric data written in Situs format to file %s \n", vol_file);
151 
152  /* header information */
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);
157 
158  return;
159 }
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
fprintf
fprintf(glob_prnt.io, "\)

◆ write_spider()

void write_spider ( char *  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
double *   
)

Definition at line 2072 of file lib_vio.cpp.

2075 {
2076  FILE *fout;
2077  unsigned long count, pnvox;
2078  float nslice, nrow, iform, imami, nsam, headrec;
2079  double dmax, dmin, dav, dsig;
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;
2083  float gamma;
2084  float dummy;
2085  int i, headlen, wmaperr;
2086 
2087  fout = fopen(vol_file, "wb");
2088  if (fout == NULL) {
2089  error_open_filename(71460, "lib_vio", vol_file);
2090  }
2091 
2092  /* set map parameters */
2093  pnvox = pextx * pexty * pextz;
2094  nsam = pextx;
2095  nrow = pexty;
2096  nslice = pextz;
2097  iform = 3.0f;
2098  imami = 1.0f;
2099  istack = 0.0f;
2100  inuse = -1.0f;
2101  maxim = 0.0f;
2102  kangle = 0.0f;
2103  phi1 = 0.0f;
2104  theta1 = 0.0f;
2105  psi1 = 0.0f;
2106  phi2 = 0.0f;
2107  theta2 = 0.0f;
2108  psi2 = 0.0f;
2109  scale = 0.0f;
2110  dummy = 0.0f;
2111  phi = 0;
2112  theta = 0;
2113  gamma = 0;
2114  xoff = 0;
2115  yoff = 0;
2116  zoff = 0;
2117  iangle = 0;
2118 
2119  /* compute density info */
2120  calc_map_info(pphi, pnvox, &dmax, &dmin, &dav, &dsig);
2121  fmax = dmax;
2122  fmin = dmin;
2123  fav = dav;
2124  fsig = dsig;
2125 
2126  /* write map */
2127  headrec = ceil(256.0f / (pextx * 1.0f));
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");
2134  wmaperr = 0;
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;
2173  }
2174 
2175  if (wmaperr != 0) {
2176  error_write_filename(71470, "lib_vio");
2177  }
2178  fclose(fout);
2179 
2180  /* print some info */
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);
2215 }
error_write_filename
void error_write_filename(int error_number, const char *program)
Definition: lib_err.cpp:227
gamma
double * gamma
Definition: numerical_recipes.cpp:7597
i
#define i
Definition: numerical_recipes.cpp:2493
calc_map_info
void calc_map_info(double *phi, unsigned long nvox, double *maxdens, double *mindens, double *ave, double *sig)
Definition: lib_vwk.cpp:491
theta
double theta
Definition: numerical_recipes.cpp:6214
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
dmin
#define dmin(a, b)
Definition: numerical_recipes.cpp:2112
f
double * f
Definition: numerical_recipes.cpp:4414
dummy
double dummy
Definition: numerical_recipes.cpp:4635
dmax
#define dmax(a, b)
Definition: numerical_recipes.cpp:2113
fmax
double fmax
Definition: numerical_recipes.cpp:5166

◆ write_vol()

void write_vol ( char *  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
double *   
)

Definition at line 65 of file lib_vio.cpp.

68 {
69 
70  if (have_situs_suffix(vol_file))
71  write_situs(vol_file, width, origx, origy, origz, extx, exty, extz, phi);
72  else
73  write_mrc(1, vol_file, width, origx, origy, origz, extx, exty, extz, phi);
74 
75  return;
76 }
write_situs
void write_situs(char *vol_file, double width, double origx, double origy, double origz, unsigned extx, unsigned exty, unsigned extz, double *phi)
Definition: lib_vio.cpp:127
write_mrc
void write_mrc(int automode, char *vol_file, double pwidth, double porigx, double porigy, double porigz, unsigned pextx, unsigned pexty, unsigned pextz, double *pphi)
Definition: lib_vio.cpp:1646
have_situs_suffix
int have_situs_suffix(char *vol_file)
Definition: lib_vio.cpp:981

◆ write_xplor()

void write_xplor ( char *  ,
double  ,
double  ,
double  ,
double  ,
unsigned  ,
unsigned  ,
unsigned  ,
double *   
)

Definition at line 1565 of file lib_vio.cpp.

1567 {
1568  FILE *fout;
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;
1574  double *pphi2;
1575 
1576  fout = fopen(vol_file, "w");
1577  if (fout == NULL) {
1578  error_open_filename(71410, "lib_vio", vol_file);
1579  }
1580 
1581  /* X-PLOR does not support free origin definitions, must work within indexing constraints */
1582  /* if necessary, bring into register with coordinate system origin (for integer start indices) */
1583  if (test_registration(porigx, porigy, porigz, pwidth) == 0) {
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");
1586  }
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);
1590 
1591  /* compute indices */
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;
1598 
1599  /* write map */
1600  printf("lib_vio>\n");
1601  printf("lib_vio> Writing X-PLOR formatted (ASCII) volumetric map \n");
1602  fprintf(fout, " \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);
1608  fprintf(fout, "ZYX\n");
1609  for (indz = 0; indz < extz2; indz++) {
1610  fprintf(fout, "%8d\n", indz);
1611  count = 0;
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)));
1614  else fprintf(fout, "%12.5E", *(pphi2 + gidz_general(indz, indy, indx, exty2, extx2)));
1615  ++count;
1616  }
1617  if ((count) % 6 != 0) fprintf(fout, " \n");
1618  }
1619  fprintf(fout, "%8d\n", -9999);
1620  fclose(fout);
1621 
1622  /* print some info */
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);
1641 }
test_registration
int test_registration(float origx, float origy, float origz, float width)
Definition: lib_vio.cpp:917
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
gidz_general
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
Definition: lib_vwk.cpp:32
interpolate_map
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)
Definition: lib_vwk.cpp:82
error_open_filename
void error_open_filename(int error_number, const char *program, char *argv)
Definition: lib_err.cpp:67
fprintf
fprintf(glob_prnt.io, "\)

◆ xplor_skip_to_number()

void xplor_skip_to_number ( FILE **  ,
char **   
)

Definition at line 327 of file lib_vio.cpp.

328 {
329  int i, done, foundletter;
330 
331  for (done = 0; done == 0;) { /* read next line */
332  if (fgets(*nextline, FLENGTH, *fin) == NULL) {
333  error_EOF(70510, "lib_vio");
334  }
335  for (i = 0; i < FLENGTH; ++i)
336  if ((*(*nextline + i) == '!') || (*(*nextline + i) == '\n')) {
337  *(*nextline + i) = '\0';
338  break;
339  }
340  foundletter = 0;
341  for (i = 0; * (*nextline + i) != '\0'; ++i) { /* check if it contains ABC FGHIJKLMNOPQRSTUVWXYZ (or lower case) */
342  if (*(*nextline + i) > 64 && *(*nextline + i) < 68) {
343  foundletter = 1;
344  break;
345  }
346  if (*(*nextline + i) > 69 && *(*nextline + i) < 91) {
347  foundletter = 1;
348  break;
349  }
350  if (*(*nextline + i) > 96 && *(*nextline + i) < 100) {
351  foundletter = 1;
352  break;
353  }
354  if (*(*nextline + i) > 101 && *(*nextline + i) < 123) {
355  foundletter = 1;
356  break;
357  }
358  }
359  if (foundletter == 0) for (i = 0; * (*nextline + i) != '\0'; ++i)
360  if (*(*nextline + i) >= '0' && *(*nextline + i) <= '9') {
361  done = 1;
362  break;
363  }
364  }
365 }
error_EOF
void error_EOF(int error_number, const char *program)
Definition: lib_err.cpp:101
done
done
Definition: numerical_recipes.cpp:7636
i
#define i
Definition: numerical_recipes.cpp:2493
FLENGTH
#define FLENGTH
Definition: lib_vio.cpp:22
Generated by   doxygen 1.8.13