20#include "pism/util/connected_components/label_components.hh"
21#include "pism/util/connected_components/label_components_impl.hh"
22#include "pism/util/Grid.hh"
23#include "pism/util/array/Scalar.hh"
37 int i_first = grid.xs();
38 return (i == i_first) and (i != 0);
43 int i_last = grid.xs() + grid.xm() - 1;
44 return (i == i_last) and (i != (
int)grid.Mx() - 1);
49 int j_last = grid.ys() + grid.ym() - 1;
50 return (j == j_last) and (j != (
int)grid.My() - 1);
55 int j_first = grid.ys();
56 return (j == j_first) and (j != 0);
64 auto grid = mask.
grid();
65 std::map<int, std::set<int> > graph;
71 auto maybe_add_edge = [&](
int a,
int b) {
72 if (a > 0 and b > 0) {
73 graph[std::min(a, b)].insert(std::max(a, b));
78 const int i = p.i(), j = p.j();
91 maybe_add_edge(M, mask.
as_int(i, j + 1));
95 maybe_add_edge(M, mask.
as_int(i + 1, j));
99 maybe_add_edge(M, mask.
as_int(i, j - 1));
103 maybe_add_edge(M, mask.
as_int(i - 1, j));
113static std::vector<int>
pack_data(
const std::map<
int, std::set<int> > &graph) {
116 std::vector<int> result{};
117 for (
const auto &e : graph) {
119 for (
auto b : e.second) {
113static std::vector<int>
pack_data(
const std::map<
int, std::set<int> > &graph) {
…}
130static std::map<int, std::set<int> >
unpack_data(
const std::vector<int> &buffer) {
131 std::map<int, std::set<int> > result;
132 for (
size_t n = 0;
n < buffer.size() / 2; ++
n) {
133 int a = buffer[2 *
n + 0];
134 int b = buffer[2 *
n + 1];
130static std::map<int, std::set<int> >
unpack_data(
const std::vector<int> &buffer) {
…}
150 const std::map<
int, std::set<int> > &subgraph) {
152 MPI_Comm_size(comm, &size);
154 MPI_Comm_rank(comm, &rank);
159 std::vector<int> message_length(size);
160 message_length[rank] = (
int)message.size();
161 MPI_Allgather(MPI_IN_PLACE, 1, MPI_INT, message_length.data(), 1, MPI_INT, comm);
164 std::vector<int> offsets(size);
166 for (
int k = 1;
k < size; ++
k) {
167 offsets[
k] = offsets[
k - 1] + message_length[
k - 1];
171 int buffer_size = offsets.back() + message_length.back();
172 std::vector<int> buffer(buffer_size);
175 MPI_Allgatherv(message.data(), (
int)message.size(), MPI_INT, buffer.data(), message_length.data(),
176 offsets.data(), MPI_INT, comm);
192static std::map<int, int>
resolve_labels(
const std::map<
int, std::set<int> > &graph,
193 bool mark_isolated_patches) {
196 std::map<int, bool> visited;
198 std::map<int, int> result;
201 const int attached = 1;
206 for (
const auto &edge : graph) {
216 std::set<int> equivalent{};
218 std::queue<int> queue;
220 while (not queue.empty()) {
221 int v = queue.front();
230 label = (v % 2 == 1) ? attached : std::min(label, v);
232 equivalent.insert(v);
235 auto it = graph.find(v);
236 if (it != graph.end()) {
237 for (
int w : it->second) {
245 for (
int v : equivalent) {
250 if (mark_isolated_patches) {
252 for (
auto &p : result) {
253 p.second = (p.second == attached) ? 0 : 1;
257 int current_label = 1;
258 std::map<int, int> new_labels;
259 for (
auto &p : result) {
260 if (new_labels.find(p.second) == new_labels.end()) {
261 new_labels[p.second] = current_label++;
263 p.second = new_labels[p.second];
192static std::map<int, int>
resolve_labels(
const std::map<
int, std::set<int> > &graph, {
…}
277 if (labels.empty()) {
284 const int i = p.i(), j = p.j();
286 int old_label = mask.
as_int(i, j);
288 auto it = labels.find(old_label);
289 if (it == labels.end()) {
293 mask(i, j) = it->second;
304 bool mark_isolated_patches) {
322 std::map<int, int> labels;
324 MPI_Comm comm = MPI_COMM_NULL;
325 MPI_Comm_split(input.
grid()->com, subdomain_is_not_empty ? 1 : 0, 0, &comm);
328 if (subdomain_is_not_empty) {
334 MPI_Comm_free(&comm);
345 int exponent =
static_cast<int>(std::ceil(std::log10(grid.max_patch_size())));
346 int labels_per_rank =
static_cast<int>(std::pow(10, exponent));
350 return grid.rank() * labels_per_rank + 1;
355namespace connected_components {
364 bool mark_isolated_patches =
false;
372 bool mark_isolated_patches =
true;
Describes the PISM grid and the distribution of data across processors.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
std::shared_ptr< const Grid > grid() const
void update_ghosts()
Updates ghost points.
int as_int(int i, int j) const
Adds "foregrdound" and "attached" concepts to an "array".
void label_isolated(array::Scalar1 &mask, int reachable)
void label(array::Scalar1 &mask)
connected_components::details::Mask< Array > Mask
static bool west_boundary(const Grid &grid, int i)
West boundary of a sub-domain.
int first_label(const Grid &grid)
static bool south_boundary(const Grid &grid, int j)
South boundary of a sub-domain.
void relabel(array::Scalar &mask, const std::map< int, int > &labels)
Assign new labels to elements of mask. Does not touch background grid cells.
static std::map< int, std::set< int > > unpack_data(const std::vector< int > &buffer)
static bool east_boundary(const Grid &grid, int i)
East boundary of a sub-domain.
std::map< int, int > final_labels(array::Scalar1 &input, bool subdomain_is_not_empty, bool mark_isolated_patches)
static std::map< int, int > resolve_labels(const std::map< int, std::set< int > > &graph, bool mark_isolated_patches)
static bool north_boundary(const Grid &grid, int j)
North boundary of a sub-domain.
static std::vector< int > pack_data(const std::map< int, std::set< int > > &graph)
static std::map< int, std::set< int > > assemble_graph(MPI_Comm comm, const std::map< int, std::set< int > > &subgraph)
static std::map< int, std::set< int > > detect_connections(array::Scalar1 &mask)
void label_components_impl(const T &input, bool mark_isolated_patches, array::Scalar1 &output)