Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/percolation/flood_fill.cpp
983bd1a55ab0features/netcdf-mpi

Last Change
Raw File

flood_fill.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* 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 <queue>
/* -------------------------------------------------------------------------- */
namespace tamaas {
// Both i and n should be signed, otherwise weird things happen
Int wrap_pbc(Int i, Int n) { return ((i % n) + n) % n; }
template <std::size_t dim>
std::array<UInt, dim> wrap_pbc(const std::array<Int, dim>& t,
const std::array<Int, dim>& n) {
std::array<UInt, dim> index;
for (UInt i = 0; i < dim; ++i)
index[i] = static_cast<UInt>(wrap_pbc(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] = static_cast<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<Int, 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<Int, 1>& /*p*/) {
return std::vector<Point>();
}
template <>
auto Cluster<2>::getNextNeighbors(const std::array<Int, 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<Int, 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<Int, 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<Int, 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>
typename Cluster<dim>::BBox Cluster<dim>::boundingBox() const {
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>
std::array<Int, dim> Cluster<dim>::extent() const {
const auto bb = boundingBox();
std::array<Int, dim> extent;
std::transform(bb.first.begin(), bb.first.end(), bb.second.begin(),
extent.begin(),
[](auto min, auto max) { return max - min + 1; });
return extent;
}
template <UInt dim>
Cluster<dim>::Cluster(Point start, const Grid<bool, dim>& map,
Grid<bool, dim>& visited, bool diagonal) {
// Visiting queue: a queue here prioritizes visiting local neighbors
// instead of growing a stack with each new visited point it gets rid of older
// neighbors first this prevents the points positions to be shifted
// indefinitely across periodic boundaries
std::queue<Point> visiting;
visiting.push(start);
// Contact sizes
const auto& n = cast_int(map.sizes());
while (not visiting.empty()) {
auto p = visiting.front();
auto p_wrapped = wrap_pbc(p, n);
if (visited(p_wrapped)) {
visiting.pop();
continue;
}
visited(p_wrapped) = true;
points.push_back(visiting.front());
visiting.pop();
auto process = [&](const std::vector<Point>& neighbors,
const bool do_perimeter) {
for (const auto& p : neighbors) {
const auto index = wrap_pbc(p, n);
if (not visited(index) and map(index)) {
visiting.push(p);
}
perimeter += do_perimeter and (not map(index));
}
};
process(getNextNeighbors(p), true);
if (diagonal) {
process(getDiagonalNeighbors(p), false);
}
}
}
typename FloodFill::List<1> FloodFill::getSegments(const Grid<bool, 1>& map) {
auto n = map.sizes();
Grid<bool, 1> visited(n, 1);
visited = false;
List<1> clusters;
for (UInt i = 0; i < n[0]; ++i) {
if (map(i) && !visited(i)) {
clusters.emplace_back(std::array<Int, 1>{(Int)i}, map, visited, false);
}
}
return clusters;
}
typename FloodFill::List<2> FloodFill::getClusters(const Grid<bool, 2>& map,
bool diagonal) {
auto global_contact = Partitioner<2>::gather(map);
auto n = global_contact.sizes();
Grid<bool, 2> visited(n, 1);
visited = false;
List<2> clusters;
if (mpi::rank() == 0) {
Logger().get(LogLevel::debug) << "Flood fill start\n";
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);
}
}
}
Logger().get(LogLevel::debug)
<< "Flood fill end, found " << clusters.size() << " clusters\n";
}
return clusters;
}
typename FloodFill::List<3> FloodFill::getVolumes(const Grid<bool, 3>& map,
bool diagonal) {
auto n = map.sizes();
Grid<bool, 3> visited(n, 1);
visited = false;
List<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 (map(i, j, k) && !visited(i, j, k)) {
clusters.emplace_back(std::array<Int, 3>{(Int)i, (Int)j, (Int)k}, map,
visited, diagonal);
}
}
}
}
return clusters;
}
template class Cluster<1>;
template class Cluster<2>;
template class Cluster<3>;
} // namespace tamaas
c4science · Help