in layer0/Isosurf.cpp [1075:1528]
static int IsosurfGradients(PyMOLGlobals * G, CSetting * set1, CSetting * set2,
CIsosurf * II, Isofield * field,
int *range, float min_level, float max_level)
{
CIsosurf *I = II;
int ok = true;
/* use local copies for performance reasons */
int n_seg = I->NSeg;
int n_line = I->NLine;
float *i_line = I->Line;
CField *i_data = I->Data;
int *i_num = I->Num;
/* get cascaded state, object, or global settings */
int spacing = SettingGet_i(G, set1, set2, cSetting_gradient_spacing);
float step_size = SettingGet_f(G, set1, set2, cSetting_gradient_step_size);
float max_walk = SettingGet_f(G, set1, set2, cSetting_gradient_max_length);
float min_walk = SettingGet_f(G, set1, set2, cSetting_gradient_min_length);
float min_slope = SettingGet_f(G, set1, set2, cSetting_gradient_min_slope);
float min_dot = SettingGet_f(G, set1, set2, cSetting_gradient_normal_min_dot);
float symmetry = SettingGet_f(G, set1, set2, cSetting_gradient_symmetry);
int symmetry_flag = false; /* are we searching for symmetric segments? */
if(symmetry != 0.0F)
symmetry_flag = true; /* (very slow process) */
if(symmetry > 1.0F)
symmetry = 1.0F / symmetry;
/* clamp dangerous parameters */
if(step_size < 0.01F)
step_size = 0.01F;
if(min_slope < 0.00001F)
min_slope = 0.00001F;
/* make sure we have gradients available for map */
if(!field->gradients)
IsofieldComputeGradients(G, field);
/* and that map has a minimum size */
if(field->gradients) {
/* locals for performance */
CField *gradients = field->gradients;
CField *points = field->points;
/* flags marking excluded regions to avoid (currently wasteful) */
int *flag = NULL;
/* variable length array for recording segment paths */
int *active_cell = VLAlloc(int, 1000);
/* ordered list of coordinates for processing */
int *order = NULL;
int range_size; /* total points in region being drawn */
int range_dim[3]; /* dimension of drawn region */
int flag_stride[3]; /* stride values for flag array */
range_dim[0] = (range[3] - range[0]);
range_dim[1] = (range[4] - range[1]);
range_dim[2] = (range[5] - range[2]);
range_size = range_dim[0] * range_dim[1] * range_dim[2];
flag = Calloc(int, range_size);
flag_stride[0] = 1;
flag_stride[1] = range_dim[0];
flag_stride[2] = range_dim[0] * range_dim[1];
order = Calloc(int, 3 * range_size);
if(order && flag && (range_dim[0] > 1) && (range_dim[1] > 1) && (range_dim[2] > 1)) {
{
/* compute approximate cell spacing */
float average_cell_axis_dist;
float *pos[4];
pos[0] = Ffloat4p(I->Coord, 0, 0, 0, 0);
pos[1] = Ffloat4p(I->Coord, 1, 0, 0, 0);
pos[2] = Ffloat4p(I->Coord, 0, 1, 0, 0);
pos[3] = Ffloat4p(I->Coord, 0, 0, 1, 0);
average_cell_axis_dist = (float) ((diff3f(pos[0], pos[1]) +
diff3f(pos[0], pos[2]) +
diff3f(pos[0], pos[3])) / 3.0);
/* scale parameters into cell units */
max_walk /= average_cell_axis_dist;
min_walk /= average_cell_axis_dist;
step_size /= average_cell_axis_dist;
min_slope *= average_cell_axis_dist;
}
{
/* generate randomized list of cell coordinates */
/* always use same seed for same volume */
OVRandom *my_rand = OVRandom_NewBySeed(G->Context->heap, range_size);
{
/* fill */
int i, j, k, *p = order;
for(k = range[2]; k < range[5]; k++) {
for(j = range[1]; j < range[4]; j++) {
for(i = range[0]; i < range[3]; i++) {
p[0] = i;
p[1] = j;
p[2] = k;
p += 3;
}
}
}
}
{
/* shuffle */
int a;
for(a = 0; a < range_size; a++) {
int *p = order + 3 * (int) (range_size * OVRandom_Get_float64_exc1(my_rand));
int *q = order + 3 * (int) (range_size * OVRandom_Get_float64_exc1(my_rand));
int t0 = p[0], t1 = p[1], t2 = p[2];
p[0] = q[0];
p[1] = q[1];
p[2] = q[2];
q[0] = t0;
q[1] = t1;
q[2] = t2;
}
}
OVRandom_Del(my_rand);
}
{
/* now draw our lines */
int a;
int *start_locus = order;
float prev_grad_normal[3] = { 0.0F, 0.0F, 0.0F };
for(a = 0; a < range_size; a++) {
int n_active_cell = 0; /* how many cells have we traversed */
float walk = max_walk; /* distance remaining to travel */
int abort_n_line = n_line; /* for backtracking */
int abort_n_seg = n_seg;
int pass; /* the pass are we on */
float symmetry_max = FLT_MIN; /* if we're trying to get symmetric segments */
float symmetry_min = FLT_MAX;
for(pass = 0; pass < 2; pass++) { /* one pass down the gradient, one up */
int have_prev = false; /* flag & storage for previous gradient & locus */
int *prev_locus = NULL;
int locus[3]; /* what cell are we in? */
float fract[3] = { 0.0F, 0.0F, 0.0F }; /* where in the cell are we? */
int n_vert = 0;
locus[0] = start_locus[0];
locus[1] = start_locus[1];
locus[2] = start_locus[2];
for(;;) {
/* master segment extension loop, using "break" to exit */
{
/* normalize locus and fract before each new step */
int done = false;
int b;
for(b = 0; b < 3; b++) {
while(fract[b] < 0.0F) { /* force fract >= 0.0 */
fract[b] += 1.0F;
locus[b]--;
}
while(fract[b] >= 1.0F) { /* force fract into [0.0-1.0) */
fract[b] -= 1.0F;
locus[b]++;
}
while(locus[b] > (range[b + 3] - 2)) { /* above range? done */
if(fract[b] <= 0.0F) {
locus[b]--;
fract[b] += 1.0F;
if(locus[b] < range[b]) { /* below range? done */
done = true;
break;
}
} else {
done = true;
break;
}
}
while(locus[b] < range[b]) { /* below range? done */
if(fract[b] > 1.0F) {
locus[b]++;
fract[b] -= 1.0F;
if(locus[b] > (range[b + 3] - 2)) { /* above range? done */
done = true;
break;
}
} else {
done = true;
break;
}
}
}
if(done)
break;
}
/* have we moved cells? */
if((!have_prev) || (have_prev && ((locus[0] != prev_locus[0]) ||
(locus[1] != prev_locus[1]) ||
(locus[2] != prev_locus[2])))) {
/* above: prev_locus may be NULL, so relying upon shortcut logic eval */
/* stop if we hit a flagged cell (flag always in lower corner) */
if(*(flag + (((locus[0] - range[0]) * flag_stride[0]) +
((locus[1] - range[1]) * flag_stride[1]) +
((locus[2] - range[2]) * flag_stride[2])))) {
break;
}
}
{
/* stop if level exceeds desired ranges */
float level = FieldInterpolatef(i_data, locus[0], locus[1], locus[2],
fract[0], fract[1], fract[2]);
if((level < min_level) || (level > max_level))
break;
if(symmetry_flag) {
if(symmetry_min > level)
symmetry_min = level;
if(symmetry_max < level)
symmetry_max = level;
}
}
{
/* interpolate gradient relative to grid */
float interp_gradient[3];
FieldInterpolate3f(gradients, locus, fract, interp_gradient);
if(length3f(interp_gradient) < min_slope) {
/* if region is too flat, then bail */
break;
}
{
/* add a line vertex at this point */
{
float *f;
VLACheck(i_line, float, n_line * 3 + 2);
f = i_line + (n_line * 3);
FieldInterpolate3f(points, locus, fract, f);
n_line++;
n_vert++;
}
/* record locus for subsequent oblation */
if((!have_prev) || (have_prev && ((locus[0] != prev_locus[0]) ||
(locus[1] != prev_locus[1]) ||
(locus[2] != prev_locus[2])))) {
VLACheck(active_cell, int, n_active_cell * 3 + 2);
{
int *xrd = active_cell + (n_active_cell * 3);
xrd[0] = locus[0];
xrd[1] = locus[1];
xrd[2] = locus[2];
n_active_cell++;
prev_locus = xrd; /* warning: volatile pointer */
}
}
/* adjust length of gradient vector */
normalize3f(interp_gradient);
/* make sure gradient isn't too divergent to take another step */
if(have_prev) {
if(dot_product3f(interp_gradient, prev_grad_normal) < min_dot)
break;
}
/* take another step */
copy3f(interp_gradient, prev_grad_normal);
/* scale and flip sign */
if(pass) {
scale3f(interp_gradient, -step_size, interp_gradient);
} else {
scale3f(interp_gradient, step_size, interp_gradient);
}
/* record progress */
walk -= step_size;
/* leave if max_walk is reached */
if(walk < 0.0F) {
break;
} else {
/* otherwise move */
add3f(interp_gradient, fract, fract);
}
have_prev = true;
}
}
} /* for */
if(n_vert < 2) { /* quash isolated vertices */
if(n_vert) {
n_line = i_num[n_seg];
}
} else if(n_line != i_num[n_seg]) { /* or retain count of new line segment */
VLACheck(i_num, int, n_seg + 1);
i_num[n_seg] = n_line - i_num[n_seg];
n_seg++;
i_num[n_seg] = n_line;
}
}
{
int abort_segment = false;
if(symmetry_flag) {
if((symmetry_max * symmetry_min) >= 0.0F) /* abort if not both +/- pot. sampled */
abort_segment = true;
else {
float symmetry_ratio = (float) (fabs(symmetry_max) / fabs(symmetry_min));
if(symmetry_ratio > 1.0F)
symmetry_ratio = 1.0F / symmetry_ratio;
if(symmetry_ratio < symmetry) /* abort if +/- weren't close enough in magnitude */
abort_segment = true;
}
}
if((max_walk - walk) < min_walk) { /* ignore too-short segments */
abort_segment = true;
}
if(abort_segment) {
n_seg = abort_n_seg;
n_line = abort_n_line;
i_num[n_seg] = n_line;
} else {
/* otherwise, keep line and oblate neighborhood */
int *ac = active_cell;
int b;
int cutoff_sq = spacing * spacing;
for(b = 0; b < n_active_cell; b++) {
int ii = ac[0], jj = ac[1], kk = ac[2];
int i0 = ii - spacing;
int j0 = jj - spacing;
int k0 = kk - spacing;
int i1 = ii + spacing + 1;
int j1 = jj + spacing + 1;
int k1 = kk + spacing + 1;
if(i0 < range[0])
i0 = range[0];
if(i1 >= range[3])
i1 = range[3] - 1;
if(j0 < range[1])
j0 = range[1];
if(j1 >= range[4])
j1 = range[4] - 1;
if(k0 < range[2])
k0 = range[2];
if(k1 >= range[5])
k1 = range[5] - 1;
{
int i, j, k;
int *flag1 = flag + (((i0 - range[0]) * flag_stride[0]) +
((j0 - range[1]) * flag_stride[1]) +
((k0 - range[2]) * flag_stride[2]));
/* highly optimized spherical flag-fill routine */
for(k = k0; k < k1; k++) {
int *flag2 = flag1;
int kk_sq = (kk - k);
kk_sq = kk_sq * kk_sq;
for(j = j0; j < j1; j++) {
int *flag3 = flag2;
int jj_sq = (jj - j);
jj_sq = (jj_sq * jj_sq) + kk_sq;
if(!(jj_sq > cutoff_sq)) {
for(i = i0; i < i1; i++) {
if(!*flag3) {
int tot_sq = (ii - i);
tot_sq = (tot_sq * tot_sq) + jj_sq;
if(!(tot_sq > cutoff_sq)) {
*flag3 = true;
}
}
flag3++;
} /* for i */
}
flag2 += flag_stride[1];
} /* for j */
flag1 += flag_stride[2];
} /* for k */
}
ac += 3; /* advance to next active cell */
} /* for b in active_cell */
}
}
start_locus += 3;
} /* for a in range_size */
}
}
/* purge memory */
VLAFreeP(active_cell);
FreeP(order);
FreeP(flag);
}
/* restore modified local copies */
I->Line = i_line;
I->Num = i_num;
I->NLine = n_line;
I->NSeg = n_seg;
return (ok);
}