marchingsquares.i.hh 7.03 KB
Newer Older
1

2
#include <set>
3
4
#include <vector>

5
#include "radixalgorithm/marchingsquares.hh"
6
7
8
9
#include "radixbug/bug.hh"

namespace radix
{
10
template <typename data_type>
11
12
13
void MarchingSquares<data_type>::clear_connected_component(int row, int col,
                                                           size_t label,
                                                           data_type wash_bit)
14
{
15
16
17
18
19
  if (row < 0 || row == mColumns) return;  // out of bounds
  if (col < 0 || col == mRows) return;     // out of bounds
  std::set<size_t> list;
  list.insert(mColumns * row + col);
  while (!list.empty())
20
  {
21
22
23
24
25
26
27
28
    const auto& it = list.begin();
    size_t c_i     = *it;

    // update the row
    row = c_i / mColumns;
    // upate the column
    col = c_i % mColumns;
    // clear data
29
    mBit[c_i]  = 0;
30
31
    mData[c_i] = wash_bit;
    // search neighbors
32
    for (int direction = 0; direction < 4; ++direction)
33
    {
34
35
36
37
      int nc = col + dx[direction];
      int nr = row + dy[direction];
      if (nc < 0 || nc >= mColumns) continue;  // out of bounds
      if (nr < 0 || nr >= mRows) continue;     // out of bounds
38
39
40
41
42
43
44
45
46
47
      size_t nc_i = mColumns * nr + nc;
      if (mBit[nc_i] == label)
      {
        // if we already have this cell in the list to look at
        // don't add it again
        if (list.find(nc_i) == list.end())
        {
          list.insert(nc_i);
        }
      }
48
    }
49
    list.erase(it);
50
51
52
  }
}

53
54
template <typename data_type>
std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
55
    data_type isovalue, data_type wash_bit, data_type wash_threshold)
56
{
57
58
  std::pair<size_t, size_t> start;
  std::vector<std::pair<int, int>> out;
59
60

  //
61
62
  // Get a starting point first, because if there isn't a place to start
  // then we don't care about going any further
63
  if (starting_point(isovalue, wash_threshold, start) == false)
64
  {
65
66
67
68
69
70
71
72
73
74
75
76
    radix_tagged_line("Did not find a starting point.");
    // return empty listing
    return out;
  }

  radix_tagged_line("Starting point[" << start.first << "," << start.second
                                      << "]");
  //
  // Initialize bit field
  for (size_t p_i = 0; p_i < mData.size(); ++p_i)
  {
    // select data as 1 if greater than isovalue
77
    if (mData[p_i] >= isovalue && mData[p_i] < wash_threshold)
78
    {
79
      mBit[p_i] = 1;
80
    }
81
    else if (mData[p_i] >= wash_threshold)
82
    {
83
      mBit[p_i] = 2;
84
    }
85
  }
86
87
88
89

  //
  // Walk the perimeter
  size_t row = start.first, column = start.second;
90
  radix_tagged_line("Starting (" << row << ", " << column << ")");
91
92
93
94
95
  do
  {
    step(row, column);
    // If our current point is within our image
    // add it to the list of points
96
97
98
    // We have to allow for row and column being equal to the number
    // of rows and columns to allow for traveling the boundaries
    if (column >= 0 && column <= mColumns && row >= 0 && row <= mRows)
99
    {
100
      out.push_back({row, column});
101
    }
102
103

    switch (next_step)
104
    {
105
106
107
108
109
110
111
112
113
114
115
116
117
118
      case StepDirection::Up:
        --row;
        break;
      case StepDirection::Left:
        --column;
        break;
      case StepDirection::Down:
        ++row;
        break;
      case StepDirection::Right:
        ++column;
        break;
      default:
        break;
119
    }
120
121
  } while (row != start.first || column != start.second);

122
123
124
  clear_connected_component(start.first, start.second,
                            mBit[mColumns * start.first + start.second],
                            wash_bit);
125

126
127
128
129
130
  return out;
}  // march

template <typename data_type>
bool MarchingSquares<data_type>::starting_point(
131
    data_type isovalue, data_type wash_threshold,
132
133
134
135
    std::pair<size_t, size_t>& pos) const
{
  for (size_t i = 0; i < mData.size(); ++i)
  {
136
    if (mData[i] >= isovalue && mData[i] < wash_threshold)
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    {
      // save the row
      pos.first = i / mColumns;
      // save the column
      pos.second = i % mColumns;
      return true;
    }
  }
  return false;
}
template <typename data_type>
bool MarchingSquares<data_type>::accepts(size_t r, size_t c) const
{
  // Make sure we don't pick a point out of bounds
  if (c < 0 || r < 0 || c >= mColumns || r >= mRows) return false;

  // Check the data value
154
  if (mBit[mColumns * r + c] == 1) return true;
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177

  return false;
}

template <typename data_type>
void MarchingSquares<data_type>::step(size_t r, size_t c)
{
  //
  // The meat of the marching squares algorithm
  // See https://en.wikipedia.org/wiki/Marching_squares
  // for specifics on the algorithm
  // upper relations
  bool u_left  = accepts(r - 1, c - 1);
  bool u_right = accepts(r - 1, c);
  // lower relations
  bool l_left  = accepts(r, c - 1);
  bool l_right = accepts(r, c);

  prev_step = next_step;
  int state = 0;

  if (u_left) state |= 1;
  if (u_right) state |= 2;
178
179
  if (l_right) state |= 4;
  if (l_left) state |= 8;
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  // State now contains a number between 0 and 15
  // representing our state.
  // In binary, it looks like 0000-1111 (in binary)

  // An example. Let's say the top two pixels are filled,
  // and the bottom two are empty.
  // Stepping through the if statements above with a state
  // of 0b0000 initially produces:
  // upper Left == true ==>  0b0001
  // upper Right == true ==> 0b0011
  // The others are false, so 0b0011 is our state
  // (That's a result 3)

  // Looking at the chart above, we see that state
  // corresponds to a move right, so in our switch statement
  // below, we add a case for 3, and assign Right as the
  // direction of the next step. We repeat this process
  // for all 16 states.

  // So we can use a switch statement to determine our
  // next direction based on
  switch (state)
  {
    case 1:
204
      next_step = StepDirection::Left;
205
206
      break;
    case 2:
207
      next_step = StepDirection::Up;
208
209
      break;
    case 3:
210
      next_step = StepDirection::Left;
211
212
      break;
    case 4:
213
      next_step = StepDirection::Right;
214
215
      break;
    case 5:
216
      if (prev_step == StepDirection::Down)
217
        next_step = StepDirection::Right;
218
      else
219
220
221
222
        next_step = StepDirection::Left;
      break;
    case 6:
      next_step = StepDirection::Up;
223
224
      break;
    case 7:
225
      next_step = StepDirection::Left;
226
227
228
229
230
231
232
      break;
    case 8:
      next_step = StepDirection::Down;
      break;
    case 9:
      next_step = StepDirection::Down;
      break;
233
    case 10:
234
235
236
      if (prev_step == StepDirection::Left)
        next_step = StepDirection::Down;
      else
237
238
239
240
        next_step = StepDirection::Up;
      break;
    case 11:
      next_step = StepDirection::Down;
241
242
      break;
    case 12:
243
      next_step = StepDirection::Right;
244
245
      break;
    case 13:
246
      next_step = StepDirection::Right;
247
248
      break;
    case 14:
249
      next_step = StepDirection::Up;
250
251
252
253
      break;
    default:
      next_step = StepDirection::None;
      break;
254
  }
255
256
257
258
259
260
261
262
263

  radix_tagged_line(
      "(" << r << ", " << c << ") =" << state << " "
          << ((next_step == StepDirection::Up)
                  ? "Up"
                  : ((next_step == StepDirection::Down)
                         ? "Down"
                         : (next_step == StepDirection::Left) ? "Left"
                                                              : "Right")));
264
}
265
266

}  // namespace radix