Skip to content

Commit

Permalink
Passing flood fill tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbeilstenedmands
jbeilstenedmands committed Sep 6, 2024
1 parent 5ebbeab commit 24c78e4
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 142 deletions.
152 changes: 81 additions & 71 deletions src/dials/algorithms/indexing/fastindex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,8 +252,9 @@ scitbx::af::shared<scitbx::vec3<double>> do_floodfill(scitbx::af::shared<double>
double d_min = 1.8) {
auto start = std::chrono::system_clock::now();
int n_points = 256;
// double fft_cell_length = n_points * d_min / 2;
// first calc rmsd and use this to create a binary grid
// int n_points = 256;
// double fft_cell_length = n_points * d_min / 2;
// first calc rmsd and use this to create a binary grid
double sumg = 0.0;
for (int i = 0; i < grid.size(); ++i) {
sumg += grid[i];
Expand All @@ -264,7 +265,7 @@ scitbx::af::shared<scitbx::vec3<double>> do_floodfill(scitbx::af::shared<double>
sum_delta_sq += std::pow(grid[i] - meang, 2);
}
double rmsd = std::pow(sum_delta_sq / grid.size(), 0.5);
scitbx::af::shared<int> grid_binary(256 * 256 * 256, 0);
scitbx::af::shared<int> grid_binary(n_points * n_points * n_points, 0);
double cutoff = rmsd_cutoff * rmsd;

for (int i = 0; i < grid.size(); i++) {
Expand All @@ -276,8 +277,8 @@ scitbx::af::shared<scitbx::vec3<double>> do_floodfill(scitbx::af::shared<double>
// now do flood fill
int n_voids = 0;

std::stack<int> stack;
std::vector<std::vector<int>> accumulators;
std::stack<scitbx::vec3<int>> stack;
std::vector<std::vector<scitbx::vec3<int>>> accumulators;
int target = 1;
int replacement = 2;
std::vector<int> grid_points_per_void;
Expand All @@ -289,101 +290,110 @@ scitbx::af::shared<scitbx::vec3<double>> do_floodfill(scitbx::af::shared<double>

for (int i = 0; i < grid_binary.size(); i++) {
if (grid_binary[i] == target) {
stack.push(i);
// std::cout << i << std::endl;
int x = i % n_points;
int y = (i % n_sq) / n_points;
int z = i / n_sq;
scitbx::vec3<int> xyz = {x, y, z};
stack.push(xyz);
grid_binary[i] = replacement;
std::vector<int> this_accumulator;
std::vector<scitbx::vec3<int>> this_accumulator;
accumulators.push_back(this_accumulator);
n_voids++;
grid_points_per_void.push_back(0);
int x_plus = 0;
int x_minus = 0;
int y_plus = 0;
int y_minus = 0;
int z_plus = 0;
int z_minus = 0;

while (!stack.empty()) {
int index = stack.top();
scitbx::vec3<int> this_xyz = stack.top();
stack.pop();
accumulators[accumulator_index].push_back(index);
accumulators[accumulator_index].push_back(this_xyz);
grid_points_per_void[accumulator_index]++;
// when finding nearest neighbours, need to check we don't step over the edge in
// each dimension likely not very efficient right now!
if ((index + 1) % n_points != 0) {
x_plus = index + 1;
} else {
x_plus = index + 1 - n_points;
}
if (grid_binary[x_plus] == target) {
grid_binary[x_plus] = replacement;
stack.push(x_plus);
}
if (index % n_points != 0) {
x_minus = index - 1;
} else {
x_minus = index - 1 + n_points;
}
if (grid_binary[x_minus] == target) {
grid_binary[x_minus] = replacement;
stack.push(x_minus);
}
if ((index % n_sq) < n_sq_minus_n) {
y_plus = index + n_points;
} else {
y_plus = index + n_points - n_sq;
// std::cout << "index " << index << std::endl;
/*while (index < 0){
index += total;
}
if (grid_binary[y_plus] == target) {
grid_binary[y_plus] = replacement;
stack.push(y_plus);
}
if ((index % n_sq) >= n_points) {
y_minus = index - n_points;
} else {
y_minus = index - n_points + n_sq;
while (index >= total) {
index -= total;
}*/

// std::cout << "xyz " << x << " " << y << " " << z << std::endl;
// increment x and calculate the array index
// std::cout << "xyz " << this_xyz[0] << " " << this_xyz[1] << " " <<
// this_xyz[2] << std::endl;
int x_plus = this_xyz[0] + 1;
int modx = modulo(this_xyz[0], n_points);
int mody = modulo(this_xyz[1], n_points) * n_points;
int modz = modulo(this_xyz[2], n_points) * n_sq;
int array_index = modulo(x_plus, n_points) + mody + modz;
/// std::cout << "xplus idx " << array_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {x_plus, this_xyz[1], this_xyz[2]};
stack.push(new_xyz);
}
if (grid_binary[y_minus] == target) {
grid_binary[y_minus] = replacement;
stack.push(y_minus);
int x_minus = this_xyz[0] - 1;
array_index = modulo(x_minus, n_points) + mody + modz;
/// std::cout << "xminus idx " << marray_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {x_minus, this_xyz[1], this_xyz[2]};
stack.push(new_xyz);
}
if ((index % total) < nn_sq_minus_n) {
z_plus = index + n_sq;
} else {
z_plus = index + n_sq - total;

int y_plus = this_xyz[1] + 1;
array_index = modx + (modulo(y_plus, n_points) * n_points) + modz;
/// std::cout << "xplus idx " << array_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {this_xyz[0], y_plus, this_xyz[2]};
stack.push(new_xyz);
}
if (grid_binary[z_plus] == target) {
grid_binary[z_plus] = replacement;
stack.push(z_plus);
int y_minus = this_xyz[1] - 1;
array_index = modx + (modulo(y_minus, n_points) * n_points) + modz;
/// std::cout << "xminus idx " << marray_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {this_xyz[0], y_minus, this_xyz[2]};
stack.push(new_xyz);
}
if ((index % total) >= n_sq) {
z_minus = index - n_sq;
} else {
z_minus = index - n_sq + total;

int z_plus = this_xyz[2] + 1;
array_index = modx + mody + (modulo(z_plus, n_points) * n_sq);
/// std::cout << "xplus idx " << array_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {this_xyz[0], this_xyz[1], z_plus};
stack.push(new_xyz);
}
if (grid_binary[z_minus] == target) {
grid_binary[z_minus] = replacement;
stack.push(z_minus);
int z_minus = this_xyz[2] - 1;
array_index = modx + mody + (modulo(z_minus, n_points) * n_sq);
/// std::cout << "xminus idx " << marray_index << std::endl;
if (grid_binary[array_index] == target) {
grid_binary[array_index] = replacement;
scitbx::vec3<int> new_xyz = {this_xyz[0], this_xyz[1], z_minus};
stack.push(new_xyz);
}
}
replacement++;
accumulator_index++;
}
}
std::cout << "n voids cpp " << n_voids << std::endl;
// DIALS_ASSERT(n_voids > 3)

// want centres_of_mass_frac and grid_points_per_void
// now calc centres of mass.
scitbx::af::shared<scitbx::vec3<double>> centres_of_mass_frac(n_voids);
for (int i = 0; i < accumulators.size(); i++) {
std::vector<int> values = accumulators[i];
std::vector<scitbx::vec3<int>> values = accumulators[i];
int n = values.size();
int divisor = n * n_points;
double x = 0.0;
double y = 0.0;
double z = 0.0;
for (int j = 0; j < n; j++) {
x += (values[j] % n_points);
y += ((values[j] % n_sq) / n_points);
z += (values[j] / n_sq);
/*std::cout << values[j] << " " << (values[j] % n_points) << " "
<< ((values[j] % n_sq) / n_points) << " " << (values[j] / n_sq)
<< std::endl;*/
x += values[j][0]; // % n_points);
y += values[j][1]; // % n_sq) / n_points);
z += values[j][2]; // / n_sq);
}
x /= divisor;
y /= divisor;
Expand Down
Loading

0 comments on commit 24c78e4

Please sign in to comment.