Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/percolation/flood_fill.cpp
1fc7644bad52features/netcdf-mpi

Last Change
Raw File

flood_fill.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "flood_fill.hh"
#include "partitioner.hh"
#include <algorithm>
#include <limits>
#include <stack>
/* -------------------------------------------------------------------------- */
namespace tamaas {
Int unsigned_modulo(Int i, UInt n) { return (i % n + n) % n; }
template <std::size_t dim>
std::array<UInt, dim> unsigned_modulo(const std::array<Int, dim>& t,
const std::array<UInt, dim>& n) {
std::array<UInt, dim> index;
for (UInt i = 0; i < dim; ++i)
index[i] = unsigned_modulo(t[i], n[i]);
return index;
}
template <std::size_t dim>
std::array<Int, dim> cast_int(const std::array<UInt, dim>& a) {
std::array<Int, dim> r;
for (UInt i = 0; i < dim; ++i)
r[i] = (Int)a[i];
return r;
}
template <UInt dim>
Cluster<dim>::Cluster(const Cluster& other)
: points(other.points), perimeter(other.perimeter) {}
template <>
auto Cluster<1>::getNextNeighbors(const std::array<UInt, 1>& p) {
std::vector<Point> neighbors(2);
Int i = p[0];
neighbors[0] = {i - 1};
neighbors[1] = {i + 1};
return neighbors;
}
template <>
auto Cluster<1>::getDiagonalNeighbors(const std::array<UInt, 1>&) {
return std::vector<Point>();
}
template <>
auto Cluster<2>::getNextNeighbors(const std::array<UInt, 2>& p) {
std::vector<Point> neighbors(4);
Int i = p[0], j = p[1];
neighbors[0] = {i + 1, j};
neighbors[1] = {i - 1, j};
neighbors[2] = {i, j - 1};
neighbors[3] = {i, j + 1};
return neighbors;
}
template <>
auto Cluster<2>::getDiagonalNeighbors(const std::array<UInt, 2>& p) {
std::vector<Point> neighbors(4);
Int i = p[0], j = p[1];
neighbors[0] = {i + 1, j + 1};
neighbors[1] = {i - 1, j - 1};
neighbors[2] = {i - 1, j + 1};
neighbors[3] = {i + 1, j - 1};
return neighbors;
}
template <>
auto Cluster<3>::getNextNeighbors(const std::array<UInt, 3>& p) {
std::vector<Point> neighbors(6);
Int i = p[0], j = p[1], k = p[2];
neighbors[0] = {i + 1, j, k};
neighbors[1] = {i - 1, j, k};
neighbors[2] = {i, j - 1, k};
neighbors[3] = {i, j + 1, k};
neighbors[4] = {i, j, k - 1};
neighbors[5] = {i, j, k + 1};
return neighbors;
}
template <>
auto Cluster<3>::getDiagonalNeighbors(const std::array<UInt, 3>& p) {
std::vector<Point> neighbors(20);
Int i = p[0], j = p[1], k = p[2];
// 8 corners
neighbors[0] = {i + 1, j + 1, k + 1};
neighbors[1] = {i + 1, j + 1, k - 1};
neighbors[2] = {i + 1, j - 1, k + 1};
neighbors[3] = {i + 1, j - 1, k - 1};
neighbors[4] = {i - 1, j + 1, k + 1};
neighbors[5] = {i - 1, j + 1, k - 1};
neighbors[6] = {i - 1, j - 1, k + 1};
neighbors[7] = {i - 1, j - 1, k - 1};
// 4 diagonals in i = 0
neighbors[8] = {i, j + 1, k + 1};
neighbors[9] = {i, j + 1, k - 1};
neighbors[10] = {i, j - 1, k + 1};
neighbors[11] = {i, j - 1, k - 1};
// 4 diagonals in j = 0
neighbors[12] = {i + 1, j, k + 1};
neighbors[13] = {i + 1, j, k - 1};
neighbors[14] = {i - 1, j, k + 1};
neighbors[15] = {i - 1, j, k - 1};
// 4 diagonals in k = 0
neighbors[16] = {i + 1, j + 1, k};
neighbors[17] = {i + 1, j - 1, k};
neighbors[18] = {i - 1, j + 1, k};
neighbors[19] = {i - 1, j - 1, k};
return neighbors;
}
template <UInt dim>
std::pair<std::array<Int, dim>, std::array<Int, dim>>
Cluster<dim>::boundingBox() const {
// TODO make sure periodic boundaries are correctly unwrapped
std::array<Int, dim> mins, maxs;
mins.fill(std::numeric_limits<Int>::max());
maxs.fill(std::numeric_limits<Int>::min());
for (const auto& p : points)
for (UInt i = 0; i < dim; ++i) {
mins[i] = std::min(mins[i], p[i]);
maxs[i] = std::max(maxs[i], p[i]);
}
return std::make_pair(mins, maxs);
}
template <UInt dim>
Cluster<dim>::Cluster(Point start, const Grid<bool, dim>& contact,
Grid<bool, dim>& visited, bool diagonal) {
// Visiting stack
std::stack<Point, std::vector<Point>> visiting;
visiting.push(start);
// Contact sizes
const auto& n = contact.sizes();
while (not visiting.empty()) {
auto p = unsigned_modulo(visiting.top(), n);
if (visited(p)) {
visiting.pop();
continue;
}
visited(p) = true;
points.push_back(visiting.top());
visiting.pop();
auto process = [&](const std::vector<Point>& neighbors,
const bool do_perimeter) {
for (auto& p : neighbors) {
auto index = unsigned_modulo(p, n);
if (not visited(index) and contact(index)) {
visiting.push(cast_int(index));
}
if (do_perimeter and (not contact(index)))
++perimeter;
}
};
process(getNextNeighbors(p), true);
if (diagonal) {
process(getDiagonalNeighbors(p), false);
}
}
}
std::list<Cluster<1>> FloodFill::getSegments(const Grid<bool, 1>& contact) {
auto n = contact.sizes();
Grid<bool, 1> visited(n, 1);
visited = false;
std::list<Cluster<1>> clusters;
for (UInt i = 0; i < n[0]; ++i) {
if (contact(i) && !visited(i)) {
clusters.emplace_back(std::array<Int, 1>{(Int)i}, contact, visited,
false);
}
}
return clusters;
}
std::list<Cluster<2>> FloodFill::getClusters(const Grid<bool, 2>& contact,
bool diagonal) {
auto global_contact = Partitioner<2>::gather(contact);
auto n = global_contact.sizes();
Grid<bool, 2> visited(n, 1);
visited = false;
std::list<Cluster<2>> clusters;
if (mpi::rank() == 0)
for (UInt i = 0; i < n[0]; ++i) {
for (UInt j = 0; j < n[1]; ++j) {
if (global_contact(i, j) && !visited(i, j)) {
clusters.emplace_back(std::array<Int, 2>{(Int)i, (Int)j},
global_contact, visited, diagonal);
}
}
}
return clusters;
}
std::list<Cluster<3>> FloodFill::getVolumes(const Grid<bool, 3>& contact,
bool diagonal) {
auto n = contact.sizes();
Grid<bool, 3> visited(n, 1);
visited = false;
std::list<Cluster<3>> clusters;
for (UInt i = 0; i < n[0]; ++i) {
for (UInt j = 0; j < n[1]; ++j) {
for (UInt k = 0; k < n[2]; ++k) {
if (contact(i, j, k) && !visited(i, j, k)) {
clusters.emplace_back(std::array<Int, 3>{(Int)i, (Int)j, (Int)k},
contact, visited, diagonal);
}
}
}
}
return clusters;
}
template class Cluster<1>;
template class Cluster<2>;
template class Cluster<3>;
} // namespace tamaas
c4science · Help