tconnected_components.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tconnected_components.cc (4428B)
---
1 /* Copyright (C) 2013, 2014, 2016 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include <vector>
21 #include <cmath>
22
23 void run_union(std::vector<unsigned int> &parent, unsigned int run1, unsigned int run2) {
24 if (parent[run1] == run2 || parent[run2] == run1) {
25 return;
26 }
27
28 while (parent[run1] != 0) {
29 run1 = parent[run1];
30 }
31
32 while (parent[run2] != 0) {
33 run2 = parent[run2];
34 }
35
36 if (run1 > run2) {
37 parent[run1] = run2;
38 } else if (run1 < run2) {
39 parent[run2] = run1;
40 } else {
41 return;
42 }
43
44 }
45
46 //! In-place labeling of connected components using a 2-scan algorithm with run-length encoding.
47 void label_connected_components(double *image, unsigned int n_rows, unsigned int n_cols, bool identify_icebergs, double mask_grounded) {
48 unsigned int max_runs = 2*n_rows;
49 const double eps = 1e-6;
50
51 std::vector<unsigned int> parents(max_runs), lengths(max_runs),
52 rows(max_runs), columns(max_runs), mask(max_runs);
53
54 unsigned int run_number = 0, r, c, parent;
55
56 // First scan
57 for (r = 0; r < n_rows; ++r) {
58 for (c = 0; c < n_cols; ++c) {
59
60 if (image[r*n_cols + c] > 0.0) {
61 // looking at a foreground pixel
62 if (c > 0 && image[r*n_cols + (c-1)] > 0.0) {
63 // one to the left is also a foreground pixel; continue the run
64 lengths[run_number] += 1;
65 } else {
66 // one to the left is a background pixel (or this is column 0); start a new run
67
68 // set the run just above as a parent (if present)
69 if (r > 0 && image[(r-1)*n_cols + c] > 0.0) {
70 parent = (unsigned int)image[(r-1)*n_cols + c];
71 } else {
72 parent = 0;
73 }
74
75 run_number += 1;
76
77 // allocate more storage (if needed)
78 if (run_number == max_runs) {
79 max_runs += n_rows;
80 parents.resize(max_runs);
81 lengths.resize(max_runs);
82 rows.resize(max_runs);
83 columns.resize(max_runs);
84 mask.resize(max_runs);
85 }
86
87 // Record this run
88 rows[run_number] = r;
89 columns[run_number] = c;
90 lengths[run_number] = 1;
91 parents[run_number] = parent;
92 mask[run_number] = 0;
93 }
94
95 if (r > 0 && image[(r-1)*n_cols + c] > 0.0) {
96 run_union(parents, (unsigned int)image[(r-1)*n_cols + c], run_number);
97 }
98
99 if (mask[run_number] == 0 && fabs(image[r*n_cols + c] - mask_grounded) < eps) {
100 mask[run_number] = 1;
101 }
102
103 image[r*n_cols + c] = run_number;
104 }
105 }
106 }
107
108 // Assign labels to runs.
109 // This uses the fact that children always follow parents,
110 // so we can do just one sweep: by the time we get to a node (run),
111 // its parent already has a final label.
112 //
113 // We use "parents" to store labels here, because once a run's label is computed
114 // we don't need to know its parent run any more.
115
116 unsigned int label = 0;
117 std::vector<unsigned int> grounded(run_number + 1);
118 for (r = 0; r <= run_number; ++r) {
119 if (parents[r] == 0) {
120 parents[r] = label;
121 label += 1;
122 } else {
123 parents[r] = parents[parents[r]];
124 }
125
126 // remember current blob (parents[r]) as "grounded" if the current run is
127 // "grounded"
128 if (mask[r] == 1) {
129 grounded[parents[r]] = 1;
130 }
131 }
132
133 // Second scan (re-label)
134 if (identify_icebergs) {
135 for (r = 0; r <= run_number; ++r) {
136 for (c = 0; c < lengths[r]; ++c) {
137 image[rows[r]*n_cols + columns[r] + c] = 1 - grounded[parents[r]];
138 }
139 }
140 } else {
141 for (r = 0; r <= run_number; ++r) {
142 for (c = 0; c < lengths[r]; ++c) {
143 image[rows[r]*n_cols + columns[r] + c] = parents[r];
144 }
145 }
146 }
147
148 // Done!
149 }