Commit e159f18d authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Updating marchingsquares to include march(const std::vector<data_type>& isovalues) method.

parent 363a6ae0
Pipeline #16254 passed with stages
in 10 minutes and 18 seconds
......@@ -44,6 +44,7 @@ class MarchingSquares
*/
bool starting_point(data_type isovalue, data_type wash_threshold,
std::pair<size_t, size_t>& pos) const;
/**
* @brief step Performs a step in the march
* @param r row aka y index
......@@ -88,7 +89,32 @@ class MarchingSquares
data_type isovalue, data_type wash_bit = {0},
data_type wash_threshold = std::numeric_limits<data_type>::max());
/**
* @brief march March on the data using multiple isovalue thresholds.
* @param isotvalues List of threshold values of the data.
* @return list of polygons for each isovalue [isovalue i][polygon i][point i]
*/
std::vector<std::vector<std::vector<std::pair<float, float>>>> march(
const std::vector<data_type>& isovalues);
void dump_component_map(std::ostream& os) const;
private:
/**
* @brief starting_point assumes labeling as occurred and finds a starting
* point which matches this label.
* @param pos std::pair<size_t, size_t> <row, column> aka <y index, x index>
* @param prev previous starting position. assumes clear_connected_components
* was called on prev
* @return bool true if a position was found
*/
bool starting_point(short label, std::pair<size_t, size_t>& pos,
std::pair<size_t, size_t> prev = {size_t(0),
size_t(0)}) const;
void step(size_t r, size_t c, short label);
bool accepts(size_t r, size_t c, short label) const;
void clear_isovalue_label(int row, int col, short label,
const std::vector<data_type>& ordered_isovalues);
}; // class
} // namespace radix
......
......@@ -3,6 +3,7 @@
#include <vector>
#include "radixalgorithm/marchingsquares.hh"
#include "radixalgorithm/ordering.hh"
#include "radixbug/bug.hh"
namespace radix
......@@ -153,6 +154,28 @@ bool MarchingSquares<data_type>::starting_point(
}
return false;
}
template <typename data_type>
bool MarchingSquares<data_type>::starting_point(
short label, std::pair<size_t, size_t>& pos,
std::pair<size_t, size_t> prev) const
{
size_t data_size = mBit.size();
size_t starti = mColumns * prev.first + prev.second;
for (size_t i = starti; i < data_size; ++i)
{
short value = mBit[i];
if (value == label)
{
// 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
{
......@@ -167,17 +190,24 @@ bool MarchingSquares<data_type>::accepts(size_t r, size_t c) const
template <typename data_type>
void MarchingSquares<data_type>::step(size_t r, size_t c)
{
// step using the default label of 1
step(r, c, 1);
}
template <typename data_type>
void MarchingSquares<data_type>::step(size_t r, size_t c, short label)
{
//
// 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);
bool u_left = accepts(r - 1, c - 1, label);
bool u_right = accepts(r - 1, c, label);
// lower relations
bool l_left = accepts(r, c - 1);
bool l_right = accepts(r, c);
bool l_left = accepts(r, c - 1, label);
bool l_right = accepts(r, c, label);
prev_step = next_step;
int state = 0;
......@@ -271,5 +301,206 @@ void MarchingSquares<data_type>::step(size_t r, size_t c)
: (next_step == StepDirection::Left) ? "Left"
: "Right")));
}
template <typename data_type>
bool MarchingSquares<data_type>::accepts(size_t r, size_t c, short label) 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
if (mBit[mColumns * r + c] >= label) return true;
return false;
}
template <typename data_type>
std::vector<std::vector<std::vector<std::pair<float, float>>>>
MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues)
{
std::pair<size_t, size_t> start, prev;
size_t isovalues_size = isovalues.size();
// [isovalue i][polygon i][point i]
std::vector<std::vector<std::vector<std::pair<float, float>>>> contours(
isovalues_size);
//
// Determine order to ensure largest to smallest ordering
// This ensures that we can deal with nested scenarios
auto comparator = [](data_type a, data_type b) { return (a > b); };
std::vector<size_t> perm = sort_permutation(isovalues, comparator);
// apply ordering so a > b > c > d ...
std::vector<data_type> ordered_isovalues = isovalues;
apply_permutation(ordered_isovalues, perm);
//
// Initialize bit field
std::vector<size_t> label_counts(isovalues_size + 1, 0);
size_t data_size = mData.size();
// data_type min_value = 9999999999;
// data_type max_value = -9999999999;
for (size_t p_i = 0; p_i < data_size; ++p_i)
{
// select data as 1 if greater than isovalue
data_type value = mData[p_i];
for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
{
if (value >= ordered_isovalues[isoi])
{
// min_value = std::min(min_value, value);
// max_value = std::max(max_value, value);
//
// save the category (aka there are only number of isovalue categories)
// specifically 1-N, 0 is background
mBit[p_i] = isovalues_size - isoi;
// label_counts[mBit[p_i]] = label_counts[mBit[p_i]] + 1;
// radix_tagged_line("assigning " << p_i << " " << value << " label("
// << mBit[p_i] << ")");
break;
}
}
}
// radix_tagged_line("Max value (" << max_value << ")");
// radix_tagged_line("Min value (" << min_value << ")");
// for (size_t i = 0; i < label_counts.size(); ++i)
//{
// radix_tagged_line("Label (" << i << ") count (" << label_counts[i] <<
// ")");
//}
for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
{
//
// since we are using permutation as an interface, to ensure a > b
// the label will be reversed to the size
short label = short(isovalues_size - isoi);
size_t contouri = size_t(label - 1);
//
// Get a starting point first, because if there isn't a place to start
// then we don't care about going any further for this isovalue
size_t polygoni = 0;
// reset previous for each contour search
prev = {0, 0};
while (starting_point(label, start, prev))
{
prev = start;
// we have a new polygon
contours[contouri].push_back(std::vector<std::pair<float, float>>());
radix_tagged_line("Starting point("
<< ordered_isovalues[isoi] << ") label (" << label
<< ")[" << start.first << "," << start.second << "]");
//
// Walk the perimeter
size_t row = start.first, column = start.second;
do
{
step(row, column, label);
// If our current point is within our image
// add it to the list of points
// 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)
{
contours[contouri][polygoni].push_back({row, column});
}
switch (next_step)
{
case StepDirection::Up:
--row;
break;
case StepDirection::Left:
--column;
break;
case StepDirection::Down:
++row;
break;
case StepDirection::Right:
++column;
break;
default:
break;
}
} while (row != start.first || column != start.second);
//
// Finish connecting the contour by making the last=first
contours[contouri][polygoni].push_back(start);
radix_tagged_line("clearing label ("
<< mBit[mColumns * start.first + start.second]
<< ") threshold (" << ordered_isovalues[isoi] << ")");
clear_isovalue_label(start.first, start.second,
mBit[mColumns * start.first + start.second],
ordered_isovalues);
// increment the polygon index
polygoni++;
}
}
return contours;
} // march
template <typename data_type>
void MarchingSquares<data_type>::clear_isovalue_label(
int row, int col, short label,
const std::vector<data_type>& ordered_isovalues)
{
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);
size_t isovalues_size = ordered_isovalues.size();
size_t isoi = isovalues_size - size_t(label) + 1;
while (!list.empty())
{
const auto& it = list.begin();
size_t c_i = *it;
// update the row
row = c_i / mColumns;
// upate the column
col = c_i % mColumns;
// re-initilize bit field
data_type value = mData[c_i];
if (label > 1) // not the last isovalue
{
radix_tagged_line("Comparing " << value
<< " >= " << ordered_isovalues[isoi]);
if (value >= ordered_isovalues[isoi])
{
//
// save the category (aka there are only number of isovalue categories)
// specifically 1-N, 0 is background
mBit[c_i] = label - 1;
radix_tagged_line("Found new isovalue (" << ordered_isovalues[isoi]
<< ") for [" << row << ","
<< col << "]");
}
}
// we didn't find a lower level contour, so zero it out
if (mBit[c_i] == label) mBit[c_i] = 0;
// search neighbors
for (int direction = 0; direction < 4; ++direction)
{
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
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);
}
}
}
list.erase(it);
}
}
} // namespace radix
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment