Commit eaa74c87 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Implement Van-Herk 1D inplace runnin max.

parent ea5d6300
#ifndef MLN_MORPHO_CANVAS_RUNNING_MAX_1D_HPP
# define MLN_MORPHO_CANVAS_RUNNING_MAX_1D_HPP
# include <mln/core/assert.hpp>
# include <algorithm>
namespace mln
{
namespace morpho
{
namespace internal
{
/// \brief The Van-herk algorithm to compute the maximum over a sliding windows
/// The algorithm runs inplace over the \p f array
///
/// It computes: f[x + offset] = Max f[y] y ∈ [x:x+k)
/// Intermediate arrays are used to store:
/// g[x] = Max f[y] y ∈ [⌊x/k⌋*k:x]
/// h[x] = Max f[y] y ∈ [x:⌈x/k⌉*k]
///
/// \param[in,out] f Input array of size \p n
/// \param[in,out] g Temporary array of size \p n
/// \param[in,out] h Temporary array of size \p n
/// \param[in] n Size of the array
/// \param[in] k Size of the sliding windows
/// \param[in] offset Offset of the sliding windows
/// \param[in] cmp Comparison function
/// \param
template <class T, class BinaryFunction>
void running_max_1d(T* f, T* g, T* h, int n, int k, int offset, BinaryFunction sup, T zero)
{
mln_precondition(n > 0);
mln_precondition(k > 0);
// Forward pass
// Compute g = Max in range (floor(x / k) * k, x)
{
for (int x = 0; x < n; x += k)
{
g[x] = f[x];
int r = std::min(n - x, k);
for (int i = 1; i < r; i++)
g[x + i] = sup(g[x + i - 1], f[x + i]);
}
}
// Backward pass
// Compute h = Max in range (x, k * ceil(x/k))
{
h[0] = f[0];
for (int x = 0; x < n; x += k)
{
int r = std::min(k, n - x - 1);
h[x + r] = f[x + r];
for (int i = r - 1; i > 0; i--)
h[x + i] = sup(h[x + i + 1], f[x + i]);
}
}
// Compute local maximum out[x] = Max f[y:y+k) with y = x + offset
// If x % k != 0:
// out[x] = Max( Max input[x : k * ceil(x/k))), Max input[x * ceil(x/k) : x + k)
// = Max( Max input[x : k * ceil(x/k))), Max input[k * floor((x+k) / k) : x + k)
// = Max( Max input[x : k * ceil(x/k))), Max input[k * floor((x+k) / k) : x + k)
// = Max( h[x], g[x+k-1] )
// Else if x % k = 0
// out[x] = Max input[x : x + k)
// = Max input[x : k * (x / k + 1)) = h[x]
// Offset the index and pointers
{
int x = 0;
int y = x + offset;
// We store in f[x] = Max f[y:y+k)
int ymax = n + offset; // Note: x < n <=> y < ymax
int yend;
// Case 1: while y + k ≤ 0 → Max f[y:y+k) = 0
{
yend = std::min(ymax, -k + 1);
for (;y < yend; x++, y++)
f[x] = zero;
}
// Case 2: While y < 0 && y + k < n
{
yend = std::min(ymax, std::min(0, n - k + 1));
for (;y < yend; y++, x++)
f[x] = g[y+k-1];
if (x == n) return;
}
// Case 2: While y < 0 && n ≤ y + k
{
yend = std::min(ymax, 0);
for (; y < yend; ++y, ++x)
f[x] = g[n-1];
if (x == n) return;
}
mln_assertion(y >= 0);
// Case 4: While 0 ≤ y && y + k ≤ n
{
yend = std::min(ymax, n - k + 1);
for (;y < yend; y++, x++)
f[x] = sup(h[y], g[y+k-1]);
if (x == n) return;
}
// Case 4: While 0 ≤ y ≤ ⌈n-1⌉ₖ < n ≤ y + k
{
yend = std::min(ymax, (n - 1) / k * k + 1);
for (;y < yend; y++, x++)
f[x] = sup(h[y], g[n-1]);
if (x == n) return;
}
// Case 4: While ⌈n-1⌉ₖ < y < n ≤ y + k
{
yend = std::min(ymax, n);
for (;y < yend; y++, x++)
f[x] = h[y];
if (x == n) return;
}
// Case 5: While n ≤ y
for (;x < n; ++x)
f[x] = zero;
}
}
} // end of namespace mln::morpho::internal
} // end of namespace mln::morpho
} // end of namespace mln
#endif //!MLN_MORPHO_CANVAS_RUNNING_MAX_1D_HPP
......@@ -5,6 +5,7 @@ add_subdirectory(maxtree)
add_subdirectory(alphatree)
add_subdirectory(component_tree)
add_core_test(${test_prefix}running_max_1d running_max_1d.cpp)
add_core_test(${test_prefix}saturate saturate.cpp)
add_core_test(${test_prefix}dilate dilate.cpp)
add_core_test(${test_prefix}erode erode.cpp)
......
#include <mln/morpho/canvas/private/running_max_1d.hpp>
#include <functional>
#include <vector>
#include <numeric>
#include <gtest/gtest.h>
template <class T, class Compare>
void running_max_1d_naive(const T* input, int size, int k, int offset, Compare cmp, T zero, T* output)
{
for (int i = 0; i < size; ++i)
{
// Compute the Max input on the range [i - offset, i - offset + k)
const T* begin = input + std::max(i + offset, 0);
const T* end = input + std::min(i + offset + k, size);
T v = (begin < end) ? *(std::max_element(begin, end, cmp)) : zero;
output[i] = v;
}
}
template <class T, class Compare>
void running_max_1d_g(const T* input, int size, int k, Compare cmp, T* output)
{
for (int i = 0; i < size; ++i)
{
// Compute the Max input on the range (floor(i / k) * k, i)
const T* begin = input + k * (i / k);
const T* end = input + std::min(i + 1, size);
const T* v = std::max_element(begin, end, cmp);
output[i] = *v;
}
}
template <class T, class Compare>
void running_max_1d_h(const T* input, int size, int k, Compare cmp, T* output)
{
for (int i = 0; i < size; ++i)
{
// Compute the Max input on the range [i, k * ceil(i / k))
const T* begin = input + i;
const int ceil = (i % k == 0) ? i : (i / k + 1) * k;
const T* end = input + std::min(ceil + 1, size);
const T* v = std::max_element(begin, end, cmp);
output[i] = *v;
}
}
class RunningMax1d : public ::testing::TestWithParam< std::tuple<int, int, int> >
{
public:
RunningMax1d(bool compute_max = true)
{
if (compute_max)
{
m_sup = [](int x, int y) { return std::max(x,y); };
m_zero = 0;
m_cmp = std::less<int>();
}
else
{
m_sup = [](int x, int y) { return std::min(x,y); };
m_zero = 255;
m_cmp = std::greater<int>();
}
}
void Check(bool increasing) const
{
int size = std::get<0>(this->GetParam());
int k = std::get<1>(this->GetParam());
int offset = std::get<2>(this->GetParam());
std::vector<int> input(size);
if (increasing)
std::iota(input.begin(), input.end(), 42);
else
std::iota(input.rbegin(), input.rend(), 42);
std::vector<int> f = input;
std::vector<int> g(size);
std::vector<int> h(size);
// Run algorithm
mln::morpho::internal::running_max_1d(f.data(), g.data(), h.data(), size, k, offset, m_sup, m_zero);
// Generate references
std::vector<int> gref(size);
std::vector<int> href(size);
std::vector<int> fref(size);
running_max_1d_g(input.data(), size, k, m_cmp, gref.data());
running_max_1d_h(input.data(), size, k, m_cmp, href.data());
running_max_1d_naive(input.data(), size, k, offset, m_cmp, m_zero, fref.data());
EXPECT_EQ(gref, g);
EXPECT_EQ(href, h);
EXPECT_EQ(fref, f);
}
private:
std::function<int (int, int)> m_sup;
int m_zero;
std::function<bool (int, int)> m_cmp;
};
class RunningMin1d : public RunningMax1d
{
public:
RunningMin1d() : RunningMax1d(false) {}
};
TEST_P(RunningMax1d, check)
{
this->Check(true);
this->Check(false);
}
TEST_P(RunningMin1d, check)
{
this->Check(true);
this->Check(false);
}
INSTANTIATE_TEST_CASE_P(se_leq_size, RunningMax1d, ::testing::Values(std::make_tuple(1, 3, 0), // Identity
std::make_tuple(12, 1, 0), // Identity
std::make_tuple(12, 3, 0), // n % k == 0
std::make_tuple(13, 3, 0), // n % k == 1
std::make_tuple(14, 3, 0), // n % 2 == 1
std::make_tuple(12, 12, 0))); // n == k
INSTANTIATE_TEST_CASE_P(se_ge_size_, RunningMax1d, ::testing::Values(std::make_tuple(12, 13, 0))); // Identity
INSTANTIATE_TEST_CASE_P(negative_offset, RunningMax1d, ::testing::Values(std::make_tuple(1, 3, -1), // Identity
std::make_tuple(12, 1, -2), // Identity
std::make_tuple(5, 7, -3), // n multiple of k
std::make_tuple(12, 3, -2), // n % k == 0
std::make_tuple(13, 3, -2), // n % k == 1
std::make_tuple(14, 3, -2), // n % k == 2
std::make_tuple(12, 12, -2))); // n == k
INSTANTIATE_TEST_CASE_P(positive_offset, RunningMax1d, ::testing::Values(std::make_tuple(1, 3, 1), // Identity
std::make_tuple(12, 1, 10), // Identity
std::make_tuple(12, 3, 10), // n % k = 0
std::make_tuple(13, 3, 10), // n % k = 1
std::make_tuple(14, 3, 10), // n % k = 2
std::make_tuple(12, 12, 10))); // n == k
INSTANTIATE_TEST_CASE_P(offset_bigger_than_size, RunningMax1d,
::testing::Values(std::make_tuple(13, 3, 14), std::make_tuple(13, 3, -14)));
INSTANTIATE_TEST_CASE_P(se_leq_size, RunningMin1d, ::testing::Values(std::make_tuple(1, 3, 0), // Identity
std::make_tuple(12, 1, 0), // Identity
std::make_tuple(12, 3, 0), // n % k = 0
std::make_tuple(13, 3, 0), // n % k = 1
std::make_tuple(14, 3, 0), // n % k = 2
std::make_tuple(12, 12, 0))); // n == k
INSTANTIATE_TEST_CASE_P(se_ge_size_, RunningMin1d, ::testing::Values(std::make_tuple(12, 13, 0))); // Identity
INSTANTIATE_TEST_CASE_P(negative_offset, RunningMin1d, ::testing::Values(std::make_tuple(1, 3, -1), // Identity
std::make_tuple(12, 1, -2), // Identity
std::make_tuple(5, 7, -3), // n multiple of k
std::make_tuple(12, 3, -2), // n multiple of k
std::make_tuple(13, 3, -2), // n multiple of k
std::make_tuple(14, 3, -2), // n multiple of k
std::make_tuple(12, 12, -2))); // n == k
INSTANTIATE_TEST_CASE_P(positive_offset, RunningMin1d, ::testing::Values(std::make_tuple(1, 3, 1), // Identity
std::make_tuple(12, 1, 10), // Identity
std::make_tuple(12, 3, 10), // n multiple of k
std::make_tuple(13, 3, 10), // n multiple of k
std::make_tuple(14, 3, 10), // n multiple of k
std::make_tuple(12, 12, 10))); // n == k
INSTANTIATE_TEST_CASE_P(offset_bigger_than_size, RunningMin1d,
::testing::Values(std::make_tuple(13, 3, 14), std::make_tuple(13, 3, -14)));
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