in layer2/ObjectMap.cpp [3767:4166]
static int ObjectMapBRIXStrToMap(ObjectMap * I, char *BRIXStr, int bytes, int state,
int quiet)
{
char *p, *pp;
float dens;
int a, b, c, d, e;
float v[3], vr[3], maxd, mind;
int ok = true;
char cc[MAXLINELEN];
ObjectMapState *ms;
int got_origin = false;
int got_extent = false;
int got_grid = false;
int got_cell = false;
int got_prod = false;
int got_plus = false;
int got_sigma = false;
int got_smiley = false;
float prod = 1.0F;
float plus = 0.0F;
float sigma = 0.0F;
float calc_sigma = 0.0F;
float calc_mean = 0.0F;
int normalize;
int swap_bytes;
normalize = SettingGetGlobal_b(I->Obj.G, cSetting_normalize_o_maps);
swap_bytes = SettingGetGlobal_b(I->Obj.G, cSetting_swap_dsn6_bytes);
if(state < 0)
state = I->NState;
if(I->NState <= state) {
VLACheck(I->State, ObjectMapState, state);
I->NState = state + 1;
}
ms = &I->State[state];
ObjectMapStateInit(I->Obj.G, ms);
maxd = -FLT_MAX;
mind = FLT_MAX;
p = BRIXStr;
ParseNCopy(cc, p, 3);
got_smiley = (strcmp(cc, ":-)") == 0);
if(got_smiley) {
/* ASCII BRIX format header */
while((*p) && (!(got_origin &&
got_extent &&
got_grid && got_cell && got_prod && got_plus && got_sigma))) {
if(!got_origin) {
pp = ParseWordCopy(cc, p, 6);
if(WordMatch(I->Obj.G, "origin", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%d", &ms->Min[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Min[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Min[2]) == 1) {
got_origin = true;
}
}
}
}
}
if(!got_extent) {
pp = ParseWordCopy(cc, p, 6);
if(WordMatch(I->Obj.G, "extent", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%d", &ms->Max[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Max[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Max[2]) == 1) {
got_extent = true;
ms->Max[0] += ms->Min[0] - 1;
ms->Max[1] += ms->Min[1] - 1;
ms->Max[2] += ms->Min[2] - 1;
}
}
}
}
}
if(!got_grid) {
pp = ParseWordCopy(cc, p, 4);
if(WordMatch(I->Obj.G, "grid", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%d", &ms->Div[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Div[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Div[2]) == 1) {
got_grid = true;
}
}
}
}
}
if(!got_cell) {
pp = ParseWordCopy(cc, p, 4);
if(WordMatch(I->Obj.G, "cell", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[2]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[2]) == 1) {
got_cell = true;
}
}
}
}
}
}
}
}
if(!got_plus) {
pp = ParseWordCopy(cc, p, 4);
if(WordMatch(I->Obj.G, "plus", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%f", &plus) == 1) {
got_plus = true;
}
}
}
if(!got_prod) {
pp = ParseWordCopy(cc, p, 4);
if(WordMatch(I->Obj.G, "prod", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%f", &prod) == 1) {
got_prod = true;
}
}
}
if(!got_sigma) {
pp = ParseWordCopy(cc, p, 5);
if(WordMatch(I->Obj.G, "sigma", cc, true) < 0) {
p = ParseWordCopy(cc, pp, 50);
if(sscanf(cc, "%f", &sigma) == 1) {
got_sigma = true;
}
}
}
p++;
}
} else {
/* Binary DSN6 format */
/* first, determine whether or not we need to byte-swap the header... */
int passed_endian_check = false;
short int *shint_ptr;
float scale_factor;
char tmp_char;
int a;
int start_swap_at = 0;
int stop_swap_at = bytes;
shint_ptr = (short int *) (p + 18 * 2);
if(*shint_ptr == 100) {
passed_endian_check = true;
start_swap_at = 512; /* don't byte-swap header */
}
if(!swap_bytes) {
stop_swap_at = 512;
}
/* for DSN6, always swap the bricks */
p = BRIXStr + start_swap_at;
for(a = start_swap_at; a < stop_swap_at; a += 2) {
tmp_char = *p;
*p = *(p + 1);
*(p + 1) = tmp_char;
p += 2;
}
p = BRIXStr;
if(*shint_ptr == 100) {
passed_endian_check = true;
}
if(!passed_endian_check) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
" Error: This looks like a DSN6 map file, but I can't match endianness.\n"
ENDFB(I->Obj.G);
} else {
shint_ptr = (short int *) p;
ms->Min[0] = *(shint_ptr++);
ms->Min[1] = *(shint_ptr++);
ms->Min[2] = *(shint_ptr++);
got_origin = true;
ms->Max[0] = *(shint_ptr++);
ms->Max[1] = *(shint_ptr++);
ms->Max[2] = *(shint_ptr++);
got_extent = true;
ms->Max[0] += ms->Min[0] - 1;
ms->Max[1] += ms->Min[1] - 1;
ms->Max[2] += ms->Min[2] - 1;
ms->Div[0] = *(shint_ptr++);
ms->Div[1] = *(shint_ptr++);
ms->Div[2] = *(shint_ptr++);
got_grid = true;
ms->Symmetry->Crystal->Dim[0] = (float) (*(shint_ptr++));
ms->Symmetry->Crystal->Dim[1] = (float) (*(shint_ptr++));
ms->Symmetry->Crystal->Dim[2] = (float) (*(shint_ptr++));
ms->Symmetry->Crystal->Angle[0] = (float) (*(shint_ptr++));
ms->Symmetry->Crystal->Angle[1] = (float) (*(shint_ptr++));
ms->Symmetry->Crystal->Angle[2] = (float) (*(shint_ptr++));
got_cell = true;
prod = (float) (*(shint_ptr++)) / 100.0F;
got_prod = true;
plus = (float) (*(shint_ptr++));
got_plus = true;
scale_factor = (float) (*(shint_ptr++));
for(a = 0; a < 3; a++) {
ms->Symmetry->Crystal->Dim[a] /= scale_factor;
ms->Symmetry->Crystal->Angle[a] /= scale_factor;
}
}
}
if(got_origin && got_extent && got_grid && got_cell && got_plus && got_prod) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Blather)
" BRIXStrToMap: Prod = %8.3f, Plus = %8.3f\n", prod, plus ENDFB(I->Obj.G);
ms->FDim[0] = ms->Max[0] - ms->Min[0] + 1;
ms->FDim[1] = ms->Max[1] - ms->Min[1] + 1;
ms->FDim[2] = ms->Max[2] - ms->Min[2] + 1;
ms->FDim[3] = 3;
if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
ok = false;
else {
SymmetryUpdate(ms->Symmetry);
ms->Field = IsosurfFieldAlloc(I->Obj.G, ms->FDim);
ms->MapSource = cMapSourceBRIX;
ms->Field->save_points = false;
{
int block_size = 8;
int xblocks = ((ms->FDim[0] - 1) / block_size) + 1;
int yblocks = ((ms->FDim[1] - 1) / block_size) + 1;
int zblocks = ((ms->FDim[2] - 1) / block_size) + 1;
int cc, bb, aa, ia, ib, ic, xa, xb, xc;
double sum = 0.0;
double sumsq = 0.0;
int n_pts = 0;
p = BRIXStr + 512;
for(cc = 0; cc < zblocks; cc++) {
for(bb = 0; bb < yblocks; bb++) {
for(aa = 0; aa < xblocks; aa++) {
for(c = 0; c < block_size; c++) {
xc = c + cc * block_size;
ic = xc + ms->Min[2];
v[2] = ic / ((float) ms->Div[2]);
for(b = 0; b < block_size; b++) {
xb = b + bb * block_size;
ib = xb + ms->Min[1];
v[1] = ib / ((float) ms->Div[1]);
for(a = 0; a < block_size; a++) {
xa = a + aa * block_size;
ia = xa + ms->Min[0];
v[0] = ia / ((float) ms->Div[0]);
dens = (((float) (*((unsigned char *) (p++)))) - plus) / prod;
if((ia <= ms->Max[0]) && (ib <= ms->Max[1]) && (ic <= ms->Max[2])) {
F3(ms->Field->data, xa, xb, xc) = dens;
sumsq += dens * dens;
sum += dens;
n_pts++;
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, xa, xb, xc, e) = vr[e];
}
}
}
}
}
}
}
}
if(n_pts > 1) {
calc_mean = (float) (sum / n_pts);
calc_sigma = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
}
}
if(ok) {
d = 0;
for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, vr);
copy3f(vr, ms->Corner + 3 * d);
d++;
}
}
}
}
if(ok && normalize) {
if(calc_sigma > R_SMALL8) {
for(c = 0; c < ms->FDim[2]; c++) {
for(b = 0; b < ms->FDim[1]; b++) {
for(a = 0; a < ms->FDim[0]; a++) {
F3(ms->Field->data, a, b, c) =
(F3(ms->Field->data, a, b, c) - calc_mean) / calc_sigma;
}
}
}
}
}
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);
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" BRIXStrToMap: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
ENDFB(I->Obj.G);
if(got_sigma) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" BRIXStrToMap: Reported Sigma = %8.3f\n", sigma ENDFB(I->Obj.G);
}
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" BRIXStrToMap: Range = %5.6f to %5.6f\n", mind, maxd ENDFB(I->Obj.G);
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" BRIXStrToMap: Calculated Mean = %8.3f, Sigma = %8.3f\n", calc_mean, calc_sigma
ENDFB(I->Obj.G);
if(normalize) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" BRIXStrToMap: Normalizing...\n" ENDFB(I->Obj.G);
}
ms->Active = true;
ObjectMapUpdateExtents(I);
}
} else {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
" Error: unable to read BRIX/DSN6 file.\n" ENDFB(I->Obj.G);
}
return (ok);
}