in layer2/ObjectMap.cpp [2381:2809]
static int ObjectMapCCP4StrToMap(ObjectMap * I, char *CCP4Str, int bytes, int state,
int quiet)
{
char *p;
int *i;
size_t bytes_per_pt;
char *q;
float dens;
int a, b, c, d, e;
float v[3], vr[3], maxd, mind;
int ok = true;
int little_endian = 1, map_endian;
/* CCP4 named from their docs */
int nc, nr, ns;
int map_mode;
int ncstart, nrstart, nsstart;
int nx, ny, nz;
float xlen, ylen, zlen, alpha, beta, gamma;
int ispg; // space group number
int sym_skip;
int mapc, mapr, maps;
int cc[3];
int n_pts;
double sum, sumsq;
float mean, stdev;
int normalize;
ObjectMapState *ms;
int expectation;
if(bytes < 256 * sizeof(int)) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
" ObjectMapCCP4: Map appears to be truncated -- aborting." ENDFB(I->Obj.G);
return (0);
}
/* state check */
if(state < 0)
state = I->NState;
/* alloc/init a new MapState */
if(I->NState <= state) {
VLACheck(I->State, ObjectMapState, state);
I->NState = state + 1;
}
ms = &I->State[state];
ObjectMapStateInit(I->Obj.G, ms);
normalize = SettingGetGlobal_b(I->Obj.G, cSetting_normalize_ccp4_maps);
p = CCP4Str;
little_endian = *((char *) &little_endian);
map_endian = (*p || *(p + 1)); // NOTE: this assumes 0x0 < NC < 0x10000
if(little_endian != map_endian) {
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: Map appears to be reverse endian, swapping...\n" ENDFB(I->Obj.G);
}
swap_endian(p, 256, sizeof(int));
}
i = (int *) p;
nc = *(i++); /* columns */
nr = *(i++); /* rows */
ns = *(i++); /* sections */
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: NC %d NR %d NS %d\n", nc, nr, ns ENDFB(I->Obj.G);
}
map_mode = *(i++); /* mode */
switch(map_mode) {
case 0:
bytes_per_pt = 1; // int8_t
break;
case 1:
bytes_per_pt = 2; // int16_t
break;
case 2:
bytes_per_pt = 4; // float
break;
default:
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
"ObjectMapCCP4-ERR: Only map mode 0-2 currently supported (this map is mode %d)",
map_mode ENDFB(I->Obj.G);
return (0);
}
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: Map is mode %d.\n", map_mode ENDFB(I->Obj.G);
}
ncstart = *(i++);
nrstart = *(i++);
nsstart = *(i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: NCSTART %d NRSTART %d NSSTART %d\n",
ncstart, nrstart, nsstart ENDFB(I->Obj.G);
}
nx = *(i++);
ny = *(i++);
nz = *(i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: NX %d NY %d NZ %d \n", nx, ny, nz ENDFB(I->Obj.G);
}
xlen = *(float *) (i++);
ylen = *(float *) (i++);
zlen = *(float *) (i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: X %8.3f Y %8.3f Z %8.3f \n", xlen, ylen, zlen ENDFB(I->Obj.G);
}
alpha = *(float *) (i++);
beta = *(float *) (i++);
gamma = *(float *) (i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: alpha %8.3f beta %8.3f gamma %8.3f \n",
alpha, beta, gamma ENDFB(I->Obj.G);
}
mapc = *(i++);
mapr = *(i++);
maps = *(i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: MAPC %d MAPR %d MAPS %d \n", mapc, mapr, maps ENDFB(I->Obj.G);
}
// AMIN, AMAX, AMEAN
mind = *(float *) (i++);
maxd = *(float *) (i++);
mean = *(float *) (i++);
ispg = *(i++);
sym_skip = *(i++);
// CCP4 Format:
// 25-37 skew transformation
// 38-52 future use "Storage space sometimes used by specific programs (default = 0)"
// MRC 2000 Format:
// 25-49 "Extra, user defined storage space (default = 0)"
// 50-52 ORIGIN
{
int skew = *(i++);
if(skew) {
double matrix[16];
auto i_float = (const float *) i;
// SKWMAT Skew matrix S (in order S11, S12, S13, S21 etc)
copy33f44d(i_float, matrix);
xx_matrix_invert(matrix, matrix, 4);
// SKWTRN Skew translation t
matrix[3] = i_float[9];
matrix[7] = i_float[10];
matrix[11] = i_float[11];
// Xo(map) = S * (Xo(atoms) - t)
ObjectStateSetMatrix(&ms->State, matrix);
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapCCP4: Applied skew transformation\n"
ENDFB(I->Obj.G);
}
}
// XORIGIN, YORIGIN, ZORIGIN (50 - 52)
// TODO See "Origin Conventions" in http://situs.biomachina.org/fmap.pdf
float * mrc2000origin = (float *)(i + 49 - 25);
if (lengthsq3f(mrc2000origin) > R_SMALL4) {
double matrix[] = {
1., 0., 0., mrc2000origin[0],
0., 1., 0., mrc2000origin[1],
0., 0., 1., mrc2000origin[2],
0., 0., 0., 1.};
ObjectStateSetMatrix(&ms->State, matrix);
if (!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapCCP4: Applied MRC 2000 ORIGIN %.2f %.2f %.2f\n",
mrc2000origin[0], mrc2000origin[1], mrc2000origin[2]
ENDFB(I->Obj.G);
}
}
i += 54 - 25;
stdev = *(float *) (i++);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: AMIN %f AMAX %f AMEAN %f ARMS %f\n", mind, maxd, mean, stdev ENDFB(I->Obj.G);
}
n_pts = nc * ns * nr;
/* at least one EM map encountered lacks NZ, so we'll try to guess it */
if((nx == ny) && (!nz)) {
nz = nx;
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapCCP4: NZ value is zero, but NX = NY, so we'll guess NZ = NX = NY.\n"
ENDFB(I->Obj.G);
}
}
expectation = sym_skip + sizeof(int) * 256 + bytes_per_pt * n_pts;
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: sym_skip %d bytes %d expectation %d\n",
sym_skip, bytes, expectation ENDFB(I->Obj.G);
}
if(bytes < expectation) {
if(bytes == (expectation - sym_skip)) {
/* accomodate bogus CCP4 map files with bad symmetry length information */
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" ObjectMapCCP4: Map has invalid symmetry length -- working around.\n"
ENDFB(I->Obj.G);
}
expectation -= sym_skip;
sym_skip = 0;
} else {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
" ObjectMapCCP4: Map appears to be truncated -- aborting.\n" ENDFB(I->Obj.G);
return (0);
}
}
q = p + (sizeof(int) * 256) + sym_skip;
if(little_endian != map_endian && bytes_per_pt > 1) {
swap_endian(q, n_pts, bytes_per_pt);
}
// with normalize == 2, use mean and stdev from file header
if(normalize == 1 && n_pts > 1) {
c = n_pts;
sum = 0.0;
sumsq = 0.0;
while(c--) {
dens = ccp4_next_value(&q, map_mode);
sumsq += dens * dens;
sum += dens;
}
mean = (float) (sum / n_pts);
stdev = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
if(stdev < 0.000001)
stdev = 1.0;
}
q = p + (sizeof(int) * 256) + sym_skip;
mapc--; /* convert to C indexing... */
mapr--;
maps--;
ms->Div[0] = nx;
ms->Div[1] = ny;
ms->Div[2] = nz;
ms->FDim[mapc] = nc;
ms->Min[mapc] = ncstart;
ms->Max[mapc] = nc + ncstart - 1;
ms->FDim[mapr] = nr;
ms->Min[mapr] = nrstart;
ms->Max[mapr] = nr + nrstart - 1;
ms->FDim[maps] = ns;
ms->Min[maps] = nsstart;
ms->Max[maps] = ns + nsstart - 1;
if(!quiet) {
if(Feedback(I->Obj.G, FB_ObjectMap, FB_Blather)) {
dump3i(ms->Div, " ObjectMapCCP4: ms->Div");
dump3i(ms->Min, " ObjectMapCCP4: ms->Min");
dump3i(ms->Max, " ObjectMapCCP4: ms->Max");
dump3i(ms->FDim, " ObjectMapCCP4: ms->FDim");
}
}
ms->Symmetry->Crystal->Dim[0] = xlen;
ms->Symmetry->Crystal->Dim[1] = ylen;
ms->Symmetry->Crystal->Dim[2] = zlen;
ms->Symmetry->Crystal->Angle[0] = alpha;
ms->Symmetry->Crystal->Angle[1] = beta;
ms->Symmetry->Crystal->Angle[2] = gamma;
if(ispg < n_space_group_numbers) {
UtilNCopy(ms->Symmetry->SpaceGroup, space_group_numbers[ispg], WordLength);
}
/* -- JV; for vol
ms->Mean = mean;
ms->SD = stdev;
ms->MaxValue = maxd;
ms->MinValue = mind;
*/
ms->FDim[3] = 3;
if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
ok = false;
else {
SymmetryUpdate(ms->Symmetry);
/* CrystalDump(ms->Crystal); */
ms->Field = IsosurfFieldAlloc(I->Obj.G, ms->FDim);
ms->MapSource = cMapSourceCCP4;
ms->Field->save_points = false;
maxd = -FLT_MAX;
mind = FLT_MAX;
for(cc[maps] = 0; cc[maps] < ms->FDim[maps]; cc[maps]++) {
v[maps] = (cc[maps] + ms->Min[maps]) / ((float) ms->Div[maps]);
for(cc[mapr] = 0; cc[mapr] < ms->FDim[mapr]; cc[mapr]++) {
v[mapr] = (cc[mapr] + ms->Min[mapr]) / ((float) ms->Div[mapr]);
for(cc[mapc] = 0; cc[mapc] < ms->FDim[mapc]; cc[mapc]++) {
v[mapc] = (cc[mapc] + ms->Min[mapc]) / ((float) ms->Div[mapc]);
dens = ccp4_next_value(&q, map_mode);
if(normalize)
dens = (dens - mean) / stdev;
F3(ms->Field->data, cc[0], cc[1], cc[2]) = dens;
if(maxd < dens)
maxd = dens;
if(mind > dens)
mind = dens;
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, vr);
for(e = 0; e < 3; e++)
F4(ms->Field->points, cc[0], cc[1], cc[2], e) = vr[e];
}
}
}
}
if(ok) {
d = 0;
for(c = 0; c < 2; c++) {
v[2] = (c * (ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
for(b = 0; b < 2; b++) {
v[1] = (b * (ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
for(a = 0; a < 2; a++) {
v[0] = (a * (ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, vr);
copy3f(vr, ms->Corner + 3 * d);
d++;
}
}
}
}
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapCCP4: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
ENDFB(I->Obj.G);
}
if(!quiet) {
if(n_pts > 1) {
if(normalize) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapCCP4: Normalizing with mean = %8.6f and stdev = %8.6f.\n",
mean, stdev ENDFB(I->Obj.G);
} else {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapCCP4: Map will not be normalized.\n ObjectMapCCP4: Current mean = %8.6f and stdev = %8.6f.\n",
mean, stdev ENDFB(I->Obj.G);
}
}
}
if(ok) {
v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, ms->ExtentMin);
v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, ms->ExtentMax);
}
#ifdef _UNDEFINED
printf("%d %d %d %d %d %d %d %d %d\n",
ms->Div[0],
ms->Min[0],
ms->Max[0],
ms->Div[1], ms->Min[1], ms->Max[1], ms->Div[2], ms->Min[2], ms->Max[2]);
printf("Okay? %d\n", ok);
fflush(stdout);
#endif
if(!ok) {
ErrMessage(I->Obj.G, "ObjectMap", "Error reading map");
} else {
ms->Active = true;
ObjectMapUpdateExtents(I);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Results)
" ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->Obj.G);
}
}
return (ok);
}