Skip to content
GitLab
Commit cd4b5423 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Merge branch '14-marching-squares' into 'master' 

Resolve "marching squares"

Closes #14

See merge request !22
parents 93edcd28 b42c00c9
Pipeline #11655 passed with stages
in 5 minutes
Hide whitespace changes
Inline Side-by-side
radixalgorithm/CMakeLists.txt
View file @ cd4b5423
......@@ -5,6 +5,8 @@ ordering.cc
)
SET(HEADERS
ordering.hh
marchingsquares.hh
marchingsquares.i.hh
)
TRIBITS_ADD_LIBRARY(radixalgorithmlib
......
radixalgorithm/marchingsquares.hh 0 → 100644
View file @ cd4b5423
#ifndef RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
#define RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
#include <limits>
#include <memory>
#include <numeric>
#include <ostream>
#include <vector>
#include "radixbug/bug.hh"
namespace radix
{
template <typename data_type>
class MarchingSquares
{
private:
std::vector<int> dx;
std::vector<int> dy;
protected:
std::vector<data_type> mData;
std::vector<int> mBit;
std::vector<int> mComponent;
size_t mRows;
size_t mColumns;
enum class StepDirection
{
Up,
Down,
Right,
Left,
None
};
StepDirection next_step;
StepDirection prev_step;
/**
* @brief starting_point find a starting position which matches or exceeds
* isovalue
* @param pos std::pair<size_t, size_t> <row, column> aka <y index, xindex>
* @return bool true if a position was found
*/
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
* @param c column aka x index
*/
void step(size_t r, size_t c);
bool accepts(size_t r, size_t c) const;
void find_components();
void dfs(int x, int y, size_t current_label);
public:
/**
* @brief MarchingSquares
* @param data const std::vector<data_type> 1D representation of 2D data
* @param rows number of rows, aka y size
* @param columns number of columns, aka x size
*/
MarchingSquares(const std::vector<data_type>& data, size_t rows,
size_t columns)
: mData(data)
, mBit(data.size(), 0)
, mRows(rows)
, mColumns(columns)
{
// radix_line("MarchingSquares data size(" << data.size() << ")");
radix_check(data.size() == (rows * columns));
dx = {+1, 0, -1, 0};
dy = {0, +1, 0, -1};
}
/**
* @brief march March the data on the isovalue threshold. Can be called
* multiple times for multiple grouping at threshold.
* @param isovalue the threshold value of the data
* @param wash_bit the value that replace all data matchings isovalue
* @param wash_threshold the threshold to which washing shouldn't apply.
* Primary use in multi-contour requests.
* @return std::vector<std::pair<int,int>> listing of all x,y indices, aka
* a closed polygon
*/
std::vector<std::pair<int, int>> march(
data_type isovalue, data_type wash_bit = {0},
data_type wash_threshold = std::numeric_limits<data_type>::max());
void dump_component_map(std::ostream& os) const;
}; // class
} // namespace radix
/** Include implementation file */
#include "radixalgorithm/marchingsquares.i.hh"
#endif /** RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_ */
radixalgorithm/marchingsquares.i.hh 0 → 100644
View file @ cd4b5423
#include <array>
#include <unordered_set>
#include <vector>
#include "radixalgorithm/marchingsquares.hh"
#include "radixbug/bug.hh"
namespace radix
{
template <typename data_type>
void MarchingSquares<data_type>::dfs(int x, int y, size_t current_label)
{
if (x < 0 || x == mColumns) return; // out of bounds
if (y < 0 || y == mRows) return; // out of bounds
size_t c_i = mColumns * y + x;
if (mComponent[c_i] || !mBit[c_i])
return; // already labeled or not marked with 1 in m
// mark the current cell
mComponent[c_i] = current_label;
// recursively mark the neighbors
for (int direction = 0; direction < 4; ++direction)
dfs(x + dx[direction], y + dy[direction], current_label);
}
template <typename data_type>
void MarchingSquares<data_type>::find_components()
{
mComponent.resize(mBit.size(), 0);
size_t component = 0;
for (size_t c_i = 0; c_i < mComponent.size(); ++c_i)
{
// save the row
size_t y = c_i / mColumns;
// save the column
size_t x = c_i % mColumns;
// check if it is already labeled
if (!mComponent[c_i] && mBit[c_i]) dfs(x, y, ++component);
}
}
template <typename data_type>
std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
data_type isovalue, data_type wash_bit, data_type wash_threshold)
{
std::pair<size_t, size_t> start;
std::vector<std::pair<int, int>> out;
//
// Get a starting point first, because if there isn't a place to start
// then we don't care about going any further
if (starting_point(isovalue, wash_threshold, start) == false)
{
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
if (mData[p_i] >= isovalue && mData[p_i] < wash_threshold)
{
mBit[p_i] = 1;
}
else if (mData[p_i] >= wash_threshold)
{
mBit[p_i] = 2;
}
}
//
// connected component labeling is required for multiple groups
if (mComponent.empty()) find_components();
size_t group_i = mColumns * start.first + start.second;
radix_line("Starting point is a member of component(" << mComponent[group_i]
<< ")");
std::cerr << std::endl;
radix_block(
std::unordered_set<int> groups(mComponent.begin(), mComponent.end()));
radix_tagged_line("Number of groups(" << groups.size() << ")");
//
// Walk the perimeter
size_t row = start.first, column = start.second;
do
{
step(row, column);
// If our current point is within our image
// add it to the list of points
if (column >= 0 && column < mColumns && row >= 0 && row < mRows)
out.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);
//
// Now we need to wash out selected group within contour, leaving those at or
// above threhold
// Get the group the starting point is a member of
for (size_t i = 0; i < mData.size(); ++i)
{
if (mComponent[i] == mComponent[group_i])
{
mData[i] = wash_bit;
}
}
//
// TODO: Simply polygon being returned with Ramer-Douglas-Peuker algorithm
return out;
} // march
template <typename data_type>
bool MarchingSquares<data_type>::starting_point(
data_type isovalue, data_type wash_threshold,
std::pair<size_t, size_t>& pos) const
{
for (size_t i = 0; i < mData.size(); ++i)
{
if (mData[i] >= isovalue && mData[i] < wash_threshold)
{
// 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
if (mData[mColumns * r + c] > 0) return true;
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;
if (l_left) state |= 4;
if (l_right) state |= 8;
// 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:
next_step = StepDirection::Up;
break;
case 2:
next_step = StepDirection::Right;
break;
case 3:
next_step = StepDirection::Right;
break;
case 4:
next_step = StepDirection::Left;
break;
case 5:
next_step = StepDirection::Up;
break;
case 6:
if (prev_step == StepDirection::Up)
{
next_step = StepDirection::Left;
}
else
{
next_step = StepDirection::Right;
}
break;
case 7:
next_step = StepDirection::Right;
break;
case 8:
next_step = StepDirection::Down;
break;
case 9:
if (prev_step == StepDirection::Right)
{
next_step = StepDirection::Up;
}
else
{
next_step = StepDirection::Down;
}
break;
case 10:
next_step = StepDirection::Down;
break;
case 11:
next_step = StepDirection::Down;
break;
case 12:
next_step = StepDirection::Left;
break;
case 13:
next_step = StepDirection::Up;
break;
case 14:
next_step = StepDirection::Left;
break;
default:
next_step = StepDirection::None;
break;
}
}
template <typename data_type>
void MarchingSquares<data_type>::dump_component_map(std::ostream& os) const
{
for (size_t i = 0; i < mComponent.size(); ++i)
{
os << std::setw(4) << mComponent[i] << ", ";
if (((i + 1) % mColumns) == 0) os << std::endl;
}
os << std::endl;
}
} // namespace radix
radixalgorithm/tests/CMakeLists.txt
View file @ cd4b5423
INCLUDE(GoogleTest)
ADD_GOOGLE_TEST(tstOrdering.cc NP 1)
ADD_GOOGLE_TEST(tstMarchingSquares.cc NP 1)
radixalgorithm/tests/MarchingSquares.4EdgeGroups.txt 0 → 100644
View file @ cd4b5423
0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,5,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,5,0,0,0,
2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,
2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,4,
2,2,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,4,4,0,
2,2,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,4,4,0,
0,0,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,0,4,4,
0,0,0,3,3,3,3,3,0,0,0,0,0,0,0,0,0,0,0,0
radixalgorithm/tests/MarchingSquares.Simple.txt 0 → 100644
View file @ cd4b5423
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,0,2,2,2,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,0,0,
0,0,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
\ No newline at end of file
radixalgorithm/tests/MarchingSquares.TwoGroups.txt 0 → 100644
View file @ cd4b5423
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,2,0,0,2,2,2,0,0,0,0,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,
0,0,2,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,0,0,
0,0,2,2,2,2,2,2,2,2,2,2,2,0,0,0,2,2,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,2,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
\ No newline at end of file
radixalgorithm/tests/tstMarchingSquares.cc 0 → 100644
View file @ cd4b5423
#include <iostream>
#include "gtest/gtest.h"
#include "radixalgorithm/marchingsquares.hh"
#include "radixbug/bug.hh"
using namespace radix;
TEST(MarchingSquares, Simple)
{
//
// define 10x20 grid of values
size_t rows{10};
size_t columns{20};
std::vector<double> grid{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,
2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 2,
2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 2, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 2, 2, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 0, 0,
0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
MarchingSquares<double> ms(grid, rows, columns);
std::vector<std::pair<int, int>> contour = ms.march(2.);
// std::cout << "[ row , col ]" << std::endl;
// std::cout << "{" << std::endl;
// for (size_t i = 0; i < contour.size(); ++i)
// {
// const auto& pixel = contour[i];
// std::cout << "{" << pixel.first << "," << pixel.second << "}";
// if (i != contour.size() - 1)
// {
// std::cout << ",";
// }
// std::cout << std::endl;
// }
// std::cout << "}" << std::endl;
std::vector<std::pair<int, int>> blessed{
{1, 2}, {2, 2}, {3, 2}, {4, 2}, {5, 2}, {6, 2}, {7, 2}, {8, 2},
{8, 3}, {8, 4}, {8, 5}, {8, 6}, {8, 7}, {8, 8}, {8, 9}, {8, 10},
{8, 11}, {8, 12}, {8, 13}, {7, 13}, {6, 13}, {5, 13}, {5, 12}, {5, 11},
{4, 11}, {4, 10}, {4, 9}, {3, 9}, {3, 8}, {3, 7}, {2, 7}, {2, 6},
{2, 5}, {3, 5}, {3, 6}, {4, 6}, {4, 7}, {4, 8}, {5, 8}, {5, 9},
{5, 10}, {6, 10}, {6, 11}, {7, 11}, {7, 10}, {7, 9}, {7, 8}, {7, 7},
{7, 6}, {7, 5}, {7, 4}, {7, 3}, {6, 3}, {5, 3}, {4, 3}, {4, 4},
{3, 4}, {2, 4}, {1, 4}, {1, 3}};
ASSERT_EQ(blessed.size(), contour.size());
for (size_t i = 0; i < blessed.size(); ++i)
{
EXPECT_EQ(blessed[i], contour[i]);
}
}
TEST(MarchingSquares, TwoGroups)
{
//
// define 10x20 grid of values
size_t rows{10};
size_t columns{20};
std::vector<double> grid{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,
2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 2,
2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 2, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 2, 2, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 0, 0,
0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 2, 2, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
MarchingSquares<double> ms(grid, rows, columns);
{
std::vector<std::pair<int, int>> first_contour = ms.march(2.);
std::vector<std::pair<int, int>> first_blessed{
{1, 2}, {2, 2}, {3, 2}, {4, 2}, {5, 2}, {6, 2}, {7, 2}, {8, 2},
{8, 3}, {8, 4}, {8, 5}, {8, 6}, {8, 7}, {8, 8}, {8, 9}, {8, 10},
{8, 11}, {8, 12}, {8, 13}, {7, 13}, {6, 13}, {5, 13}, {5, 12}, {5, 11},
{4, 11}, {4, 10}, {4, 9}, {3, 9}, {3, 8}, {3, 7}, {2, 7}, {2, 6},
{2, 5}, {3, 5}, {3, 6}, {4, 6}, {4, 7}, {4, 8}, {5, 8}, {5, 9},
{5, 10}, {6, 10}, {6, 11}, {7, 11}, {7, 10}, {7, 9}, {7, 8}, {7, 7},
{7, 6}, {7, 5}, {7, 4}, {7, 3}, {6, 3}, {5, 3}, {4, 3}, {4, 4},
{3, 4}, {2, 4}, {1, 4}, {1, 3}};
ASSERT_EQ(first_blessed.size(), first_contour.size());
for (size_t i = 0; i < first_blessed.size(); ++i)
{
EXPECT_EQ(first_blessed[i], first_contour[i]);
}
}
{ // get the second grouping
std::vector<std::pair<int, int>> second_contour = ms.march(2.);
// std::cout << "[ row , col ]" << std::endl;
// std::cout << "{" << std::endl;
// for (size_t i = 0; i < second_contour.size(); ++i)
// {
// const auto& pixel = second_contour[i];
// std::cout << "{" << pixel.first << "," << pixel.second << "}";
// if (i != second_contour.size() - 1)
// {
// std::cout << ",";
// }
// std::cout << std::endl;
// }
// std::cout << "}" << std::endl;
std::vector<std::pair<int, int>> second_blessed{
{7, 16}, {8, 16}, {9, 16}, {9, 17}, {9, 18}, {8, 18}, {7, 18}, {7, 17}};
ASSERT_EQ(second_blessed.size(), second_contour.size());
for (size_t i = 0; i < second_blessed.size(); ++i)
{
EXPECT_EQ(second_blessed[i], second_contour[i]);
}
}
//
// Expect that there is no
EXPECT_EQ(0, ms.march(2.).size());
}
TEST(MarchingSquares, 4EdgeGroups)
{
//
// define 10x20 grid of values
size_t rows{10};
size_t columns{20};
std::vector<double> grid{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 0, 0, 0, 2, 2, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 4, 4, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 4, 4, 4, 2, 2, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,
4, 0, 2, 2, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0,
0, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 3,
3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
MarchingSquares<double> ms(grid, rows, columns);
{
std::vector<std::pair<int, int>> contour = ms.march(2., 0., 3.);
// ms.dump_component_map(std::cout);
// std::cout << "{" << std::endl;
// for (size_t i = 0; i < contour.size(); ++i)
// {
// const auto& pixel = contour[i];
// std::cout << "{" << pixel.first << "," << pixel.second << "}";
// if (i != contour.size() - 1)
// {
// std::cout << ",";
// }