+
+
+
+
+
+
+
+
+
+
+
+
+
22inline std::string get_namespace()
+
+
24 std::string ret =
"GooseEYE";
+
25#ifdef GOOSEEYE_USE_XTENSOR_PYTHON
+
+
+
+
+
+
+
+
37inline std::string shape_to_string(
const T& shape)
+
+
39 std::string ret =
"[";
+
40 for (
size_t i = 0; i < shape.size(); ++i) {
+
41 ret += std::to_string(shape[i]);
+
42 if (i < shape.size() - 1) {
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
97 const std::vector<size_t>& shape,
+
+
+
+
101 bool periodic =
true)
+
+
+
+
+
+
+
+
+
110 for (
size_t i = 0; i < row.size(); ++i) {
+
111 for (
int di = -r(i); di <= r(i); ++di) {
+
112 for (
int dj = -r(i); dj <= r(i); ++dj) {
+
113 int dr = (int)(ceil(pow((
double)(pow(di, 2) + pow(dj, 2)), 0.5)));
+
+
115 out.periodic(row(i) + di, col(i) + dj) = 1;
+
+
+
+
+
+
+
+
+
124 for (
size_t i = 0; i < row.size(); ++i) {
+
125 for (
int di = -r(i); di <= r(i); ++di) {
+
126 for (
int dj = -r(i); dj <= r(i); ++dj) {
+
127 if (out.in_bounds(row(i) + di, col(i) + dj)) {
+
128 int dr = (int)(ceil(pow((
double)(pow(di, 2) + pow(dj, 2)), 0.5)));
+
+
130 out(row(i) + di, col(i) + dj) = 1;
+
+
+
+
+
+
+
+
+
+
+
+
+
148dummy_circles(
const std::vector<size_t>& shape,
bool periodic =
true, uint64_t seed = 0)
+
+
+
151 prrng::pcg32 rng(seed);
+
+
+
154 size_t N = (size_t)(0.05 * (
double)shape[0]);
+
155 size_t M = (size_t)(0.05 * (
double)shape[1]);
+
156 size_t R = (size_t)(pow((0.3 * (
double)(shape[0] * shape[1])) / (M_PI * (double)(N * M)), 0.5));
+
+
+
+
+
+
+
163 for (
size_t i = 0; i < N; i++) {
+
164 for (
size_t j = 0; j < M; j++) {
+
165 row[i * M + j] = (int)((
double)i * (double)shape[0] / (
double)N);
+
166 col[i * M + j] = (int)((
double)j * (double)shape[1] / (
double)M);
+
167 r[i * M + j] = (int)R;
+
+
+
+
+
172 int dN = (int)(0.5 * (
double)shape[0] / (double)N);
+
173 int dM = (int)(0.5 * (
double)shape[1] / (double)M);
+
+
+
176 for (
size_t i = 0; i < N * M; i++) {
+
177 row(i) += rng.randint(2 * dN) - dN;
+
178 col(i) += rng.randint(2 * dM) - dM;
+
179 r(i) = (double)r(i) * 2.0 * rng.random();
+
+
+
+
+
+
+
+
+
+
+
+
199 std::is_integral<typename T::value_type>::value &&
+
200 std::is_integral<typename S::value_type>::value,
+
+
+
+
+
+
206 bool periodic =
true);
+
+
212template <class T, std::enable_if_t<std::is_integral<typename T::value_type>::value,
int> = 0>
+
+
+
+
+
+
+
223 std::is_integral<typename T::value_type>::value &&
+
224 std::is_integral<typename S::value_type>::value,
+
+
226inline T
dilate(
const T& f,
const S& kernel,
size_t iterations = 1,
bool periodic =
true);
+
+
232template <class T, std::enable_if_t<std::is_integral<typename T::value_type>::value,
int> = 0>
+
233inline T
dilate(
const T& f,
size_t iterations = 1,
bool periodic =
true);
+
+
241template <
class T,
class S>
+
+
+
+
+
+
+
+
+
249 for (
size_t i = 0; i < a.size(); ++i) {
+
250 ret(a.flat(i)) = b.flat(i);
+
+
+
+
+
+
+
+
+
+
+
265 using value_type =
typename T::value_type;
+
266 std::map<value_type, value_type> map;
+
+
268 for (
size_t i = 0; i < a.size(); ++i) {
+
269 map.try_emplace(a.flat(i), b.flat(i));
+
+
+
+
+
274 xt::empty<typename T::value_type>(std::array<size_t, 2>{map.size(), 2});
+
+
276 for (
auto const& [key, val] : map) {
+
+
+
+
+
+
+
+
+
+
291template <
class L,
class A>
+
+
+
+
+
+
296 using value_type =
typename A::value_type;
+
297 std::map<value_type, value_type> map;
+
298 for (
size_t i = 0; i < rename.shape(0); ++i) {
+
299 map.emplace(rename(i, 0), rename(i, 1));
+
+
+
302#ifdef GOOSEEYE_ENABLE_ASSERT
+
303 auto l = xt::unique(labels);
+
304 for (
size_t i = 0; i < l.size(); ++i) {
+
+
+
+
+
309 L ret = xt::empty_like(labels);
+
310 for (
size_t i = 0; i < labels.size(); ++i) {
+
311 ret.flat(i) = map[labels.flat(i)];
+
+
+
+
+
+
+
+
+
+
+
327 using value_type =
typename T::value_type;
+
328 auto unq = xt::unique(labels);
+
329 bool background = xt::any(xt::equal(unq, 0));
+
+
331 std::array<size_t, 2> shape = {unq.size(), 2};
+
+
+
+
+
+
+
338 for (
size_t i = 1; i < unq.size(); ++i) {
+
339 rename(i, 0) = unq(i);
+
+
+
+
+
+
345 for (
size_t i = 0; i < unq.size(); ++i) {
+
+
+
+
349 rename(row, 0) = unq(i);
+
350 rename(row, 1) = row;
+
+
+
+
+
+
356 for (
size_t i = 0; i < unq.size(); ++i) {
+
357 rename(i, 0) = unq(i);
+
358 rename(i, 1) = i + 1;
+
+
+
+
+
+
+
370template <
class L,
class A>
+
+
+
+
373#ifdef GOOSEEYE_ENABLE_ASSERT
+
374 auto a = xt::unique(labels);
+
375 auto b = xt::unique(order);
+
+
+
+
+
380 auto maxlab = *std::max_element(order.begin(), order.end());
+
381 std::vector<typename A::value_type> renum(maxlab + 1);
+
+
383 for (
size_t i = 0; i < order.size(); ++i) {
+
+
+
+
387 L ret = xt::empty_like(labels);
+
388 for (
size_t i = 0; i < labels.size(); ++i) {
+
389 ret.flat(i) = renum[labels.flat(i)];
+
+
+
+
+
+
+
+
+
+
+
403 using value_type =
typename T::value_type;
+
404 std::map<value_type, value_type> map;
+
+
406 for (
size_t i = 0; i < labels.size(); ++i) {
+
407 if (map.count(labels.flat(i)) == 0) {
+
408 map.emplace(labels.flat(i), 1);
+
+
+
411 map[labels.flat(i)]++;
+
+
+
+
+
+
417 xt::empty<value_type>(std::array<size_t, 2>{map.size(), 2});
+
+
419 for (
auto const& [key, val] : map) {
+
+
+
+
+
+
+
+
+
+
+
+
435template <
size_t Dim,
class T>
+
+
+
438#ifdef GOOSEEYE_ENABLE_ASSERT
+
439 for (
size_t i = 0; i < Dim; ++i) {
+
+
+
+
+
444 std::array<size_t, Dim> mid;
+
445 for (
size_t i = 0; i < Dim; ++i) {
+
446 mid[i] = (kernel.shape(i) - 1) / 2;
+
+
+
449 for (
size_t i = 0; i < Dim; ++i) {
+
450 idx += mid[i] * kernel.strides()[i];
+
+
+
453 kernel.flat(idx) = 0;
+
+
455 if constexpr (Dim == 1) {
+
456 auto i = xt::flatten_indices(xt::argwhere(kernel)) - mid[0];
+
+
458 std::copy(i.begin(), i.end(), ret.begin());
+
+
+
+
462 auto ret = xt::from_indices(xt::argwhere(kernel));
+
463 for (
size_t i = 0; i < Dim; ++i) {
+
464 xt::view(ret, xt::all(), i) -= mid[i];
+
+
+
+
+
+
+
476template <
size_t Dimension,
bool Periodicity = true>
+
+
+
+
479 static constexpr size_t Dim = Dimension;
+
+
+
+
483 std::array<ptrdiff_t, Dim> m_shape;
+
+
+
486 ptrdiff_t m_new_label = 1;
+
+
488 std::array<ptrdiff_t, Dim> m_strides;
+
+
496 std::vector<ptrdiff_t> m_renum;
+
+
508 std::vector<ptrdiff_t> m_next;
+
509 std::vector<ptrdiff_t> m_connected;
+
+
+
+
+
+
+
+
+
+
521 if constexpr (
Dim == 1) {
+
+
+
+
525 else if constexpr (
Dim == 2) {
+
+
527 m_dx = {{-1, 0}, {0, -1}, {0, 1}, {1, 0}};
+
+
529 else if constexpr (
Dim == 3) {
+
530 m_dx = {{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}};
+
+
+
533 throw std::runtime_error(
"Please specify the kernel in dimensions > 3.");
+
+
+
+
+
+
542 template <
class T,
class K>
+
+
+
+
545 m_dx = detail::kernel_to_dx<Dim>(kernel);
+
+
+
+
+
+
+
551 void init(
const T&
shape)
+
+
553 m_label = xt::empty<ptrdiff_t>(
shape);
+
554 m_renum.resize(m_label.size() + 1);
+
555 m_next.resize(m_label.size() + 1);
+
556 for (
size_t i = 0; i <
Dim; ++i) {
+
557 m_shape[i] =
static_cast<ptrdiff_t
>(
shape[i]);
+
558 if constexpr (
Dim >= 2) {
+
559 m_strides[i] =
static_cast<ptrdiff_t
>(m_label.strides()[i]);
+
+
+
+
563 m_connected.resize(m_dx.shape(0));
+
+
+
+
567 if constexpr (
Dim == 2) {
+
568 if (m_shape[0] == 1) {
+
569 get_compare = &ClusterLabeller<Dimension, Periodicity>::get_compare_2d_1n;
+
+
571 else if (m_shape[1] == 1) {
+
572 get_compare = &ClusterLabeller<Dimension, Periodicity>::get_compare_2d_n1;
+
+
+
+
+
+
+
+
+
583 std::fill(m_label.begin(), m_label.end(), 0);
+
584 std::iota(m_renum.begin(), m_renum.end(), 0);
+
+
+
+
+
+
+
+
+
595 ptrdiff_t n =
static_cast<ptrdiff_t
>(m_new_label);
+
+
+
598 for (ptrdiff_t i = 1; i < n; ++i) {
+
599 if (m_renum[i] == i) {
+
600 m_renum[i] = m_new_label;
+
+
+
+
604 this->private_renumber(m_renum);
+
605 std::iota(m_renum.begin(), m_renum.begin() + n, 0);
+
+
+
+
+
+
+
614 std::fill(m_next.begin(), m_next.end(), -1);
+
+
+
+
622 void private_renumber(
const T& renum)
+
+
624 for (
size_t i = 0; i < m_label.size(); ++i) {
+
625 m_label.flat(i) = renum[m_label.flat(i)];
+
+
+
+
637 void merge_detail(ptrdiff_t a, ptrdiff_t b)
+
+
+
640 ptrdiff_t i = m_renum[b];
+
641 ptrdiff_t target = m_renum[a];
+
+
+
+
+
+
+
+
+
+
651 while (m_next[a] != -1) {
+
+
+
+
+
+
663 size_t unique(ptrdiff_t*
labels,
size_t nlabels)
+
+
+
+
+
+
676 ptrdiff_t merge(ptrdiff_t*
labels,
size_t nlabels)
+
+
678 ptrdiff_t target =
labels[0];
+
679 for (
size_t i = 1; i < nlabels; ++i) {
+
680 this->merge_detail(target,
labels[i]);
+
+
+
+
+
+
+
+
+
+
+
691 this->private_renumber(m_renum);
+
+
+
+
+
700 ptrdiff_t get_compare_2d_1n(
size_t idx,
size_t j)
+
+
+
703 return (m_shape[1] + idx + m_dx(j, 1)) % m_shape[1];
+
+
+
706 ptrdiff_t compare = idx + m_dx(j, 1);
+
707 if (compare < 0 || compare >= m_shape[1]) {
+
+
+
+
+
+
+
718 ptrdiff_t get_compare_2d_n1(
size_t idx,
size_t j)
+
+
+
721 return (m_shape[0] + idx + m_dx(j, 0)) % m_shape[0];
+
+
+
724 ptrdiff_t compare = idx + m_dx(j, 0);
+
725 if (compare < 0 || compare >= m_shape[0]) {
+
+
+
+
+
+
+
740 ptrdiff_t get_compare_default(
size_t idx,
size_t j)
+
+
+
743 return (m_shape[0] + idx + m_dx.flat(j)) % m_shape[0];
+
+
+
746 ptrdiff_t compare = idx + m_dx.flat(j);
+
747 if (compare < 0 || compare >= m_shape[0]) {
+
+
+
750 return idx + m_dx.flat(j);
+
+
+
753 ptrdiff_t ii = (m_shape[0] + (idx / m_strides[0]) + m_dx(j, 0)) % m_shape[0];
+
754 ptrdiff_t jj = (m_shape[1] + (idx % m_strides[0]) + m_dx(j, 1)) % m_shape[1];
+
755 return ii * m_shape[1] + jj;
+
+
+
758 ptrdiff_t ii = (idx / m_strides[0]) + m_dx(j, 0);
+
759 ptrdiff_t jj = (idx % m_strides[0]) + m_dx(j, 1);
+
760 if (ii < 0 || ii >= m_shape[0] || jj < 0 || jj >= m_shape[1]) {
+
+
+
763 return ii * m_shape[1] + jj;
+
+
+
766 auto index = xt::unravel_from_strides(idx, m_strides, xt::layout_type::row_major);
+
767 for (
size_t d = 0; d <
Dim; ++d) {
+
768 index[d] += m_dx(j, d);
+
+
770 if (index[d] < 0 || index[d] >= m_shape[d]) {
+
+
+
+
+
+
776 index[d] = (n + (index[d] % n)) % n;
+
+
+
779 return xt::ravel_index(index, m_shape, xt::layout_type::row_major);
+
+
+
+
783 void label_impl(
size_t idx)
+
+
785 size_t nconnected = 0;
+
+
787 for (
size_t j = 0; j < m_dx.shape(0); ++j) {
+
788 ptrdiff_t compare = (this->*get_compare)(idx, j);
+
+
+
+
+
+
794 if (m_label.flat(compare) != 0) {
+
795 m_connected[nconnected] = m_renum[m_label.flat(compare)];
+
+
+
+
+
800 if (nconnected == 0) {
+
801 m_label.flat(idx) = m_new_label;
+
+
+
+
+
806 if (nconnected == 1) {
+
807 m_label.flat(idx) = m_connected[0];
+
+
+
+
811 nconnected = this->unique(&m_connected[0], nconnected);
+
812 if (nconnected == 1) {
+
813 m_label.flat(idx) = m_connected[0];
+
+
+
+
+
+
+
820 m_label.flat(idx) = this->merge(&m_connected[0], nconnected);
+
+
+
+
+
+
+
+
+
832 GOOSEEYE_ASSERT(xt::has_shape(img, m_label.shape()), std::out_of_range);
+
+
834 for (
size_t idx = 0; idx < img.size(); ++idx) {
+
835 if (img.flat(idx) == 0) {
+
+
+
838 if (m_label.flat(idx) != 0) {
+
+
+
841 this->label_impl(idx);
+
+
+
+
+
+
+
+
848 bool legal_points(
const T& begin,
const T& end)
+
+
850 size_t n = m_label.size();
+
851 if constexpr (std::is_signed_v<typename T::value_type>) {
+
852 return !std::any_of(begin, end, [n](
size_t i) {
return i < 0 || i >= n; });
+
+
+
855 return !std::any_of(begin, end, [n](
size_t i) {
return i >= n; });
+
+
+
+
+
+
+
+
+
+
869 for (
auto it = begin; it != end; ++it) {
+
870 if (m_label.flat(*it) != 0) {
+
+
+
873 this->label_impl(*it);
+
+
+
+
+
+
+
+
+
+
+
886 return this->
add_points(idx.begin(), idx.end());
+
+
+
+
+
+
+
+
+
+
902 GOOSEEYE_ASSERT(this->legal_points(idx.begin(), idx.end()), std::out_of_range);
+
903 std::vector<size_t> ret;
+
+
+
906 auto nl = m_new_label;
+
+
908 auto lab = m_label.flat(idx(i));
+
+
910 for (; i < idx.size(); ++i) {
+
911 auto l = m_label.flat(idx(i));
+
912 if (l != lab && l != 0) {
+
+
+
+
+
+
+
919 this->label_impl(idx(i));
+
920 if (m_new_label != nl || m_nmerge != nm) {
+
+
+
+
+
+
926 if (i == idx.size()) {
+
+
+
+
+
+
932 ret.push_back(idx.size());
+
+
+
+
+
+
+
+
942 return detail::get_namespace() +
"ClusterLabeller" + std::to_string(
Dim) +
" " +
+
943 detail::shape_to_string(m_shape);
+
+
+
+
+
+
+
952 return m_label.shape();
+
+
+
+
+
+
+
961 return m_label.size();
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
991 :
Clusters(f, kernel::nearest(f.dimension()), periodic)
+
+
+
+
+
1002 template <
class T,
class S>
+
+
1003 Clusters(
const T& f,
const S& kernel,
bool periodic =
true) : m_periodic(periodic)
+
+
+
1006 "(please open a PR for missing functions)");
+
+
1008 static_assert(std::is_integral<typename T::value_type>::value,
"Integral labels required.");
+
1009 static_assert(std::is_integral<typename S::value_type>::value,
"Integral kernel required.");
+
+
1011 GOOSEEYE_ASSERT(xt::all(xt::equal(f, 0) || xt::equal(f, 1)), std::out_of_range);
+
1012 GOOSEEYE_ASSERT(xt::all(xt::equal(kernel, 0) || xt::equal(kernel, 1)), std::out_of_range);
+
1013 GOOSEEYE_ASSERT(f.dimension() == kernel.dimension(), std::out_of_range);
+
+
1015 m_shape = detail::shape(f);
+
1016 m_kernel = xt::atleast_3d(kernel);
+
1017 m_pad = detail::pad_width(m_kernel);
+
1018 m_l = xt::atleast_3d(f);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1033 xt::view(renum, xt::keep(labels)) = xt::arange<int>(
static_cast<int>(labels.size()));
+
1034 for (
auto& i : m_l) {
+
+
+
+
+
+
+
+
+
1045 return xt::adapt(m_l.data(), m_shape);
+
+
+
+
+
+
+
+
+
+
1057 for (
int l = 1; l < static_cast<int>(x.shape(0)); ++l) {
+
1058 c(x(l, 0), x(l, 1), x(l, 2)) = l;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1074 x = this->average_position_periodic();
+
+
+
1077 x = this->average_position(m_l);
+
+
+
+
+
+
+
+
1085 return xt::view(x, xt::all(), xt::keep(axes));
+
+
+
+
+
+
+
+
+
+
1097 for (
size_t h = 0; h < m_l.shape(0); ++h) {
+
1098 for (
size_t i = 0; i < m_l.shape(1); ++i) {
+
1099 for (
size_t j = 0; j < m_l.shape(2); ++j) {
+
1100 ret(m_l(h, i, j))++;
+
+
+
+
+
+
+
+
+
+
+
+
1111 xt::pad_mode pad_mode = xt::pad_mode::constant;
+
+
+
+
1115 pad_mode = xt::pad_mode::periodic;
+
+
+
1118 m_l = xt::pad(m_l, m_pad, pad_mode, pad_value);
+
+
+
+
+
+
+
+
+
1127 array_type::tensor<int, 1> renum = xt::arange<int>(
static_cast<int>(m_l.size()));
+
+
1129 for (
size_t h = m_pad[0][0]; h < m_l.shape(0) - m_pad[0][1]; ++h) {
+
1130 for (
size_t i = m_pad[1][0]; i < m_l.shape(1) - m_pad[1][1]; ++i) {
+
1131 for (
size_t j = m_pad[2][0]; j < m_l.shape(2) - m_pad[2][1]; ++j) {
+
+
1133 if (!m_l(h, i, j)) {
+
+
+
+
+
+
1139 xt::range(h - m_pad[0][0], h + m_pad[0][1] + 1),
+
1140 xt::range(i - m_pad[1][0], i + m_pad[1][1] + 1),
+
1141 xt::range(j - m_pad[2][0], j + m_pad[2][1] + 1));
+
+
1143 auto Ni = Li * m_kernel;
+
+
1145 int l = xt::amax(Ni)();
+
+
+
+
+
+
+
1152 Li = xt::where(xt::equal(Ni, 1), l, Li);
+
+
1154 if (xt::all(xt::equal(Li, l) || xt::equal(Li, 0) || xt::equal(Li, 1))) {
+
+
+
+
+
1159 array_type::array<int> merge = xt::where(xt::less_equal(Li, 1), l, Li);
+
1160 merge = xt::unique(merge);
+
+
1162 int linkto = xt::amin(xt::view(renum, xt::keep(merge)))[0];
+
1163 for (
auto& a : merge) {
+
1164 renum = xt::where(xt::equal(renum, renum(a)), linkto, renum);
+
+
+
+
+
+
+
+
+
1173 xt::range(m_pad[0][0], m_l.shape(0) - m_pad[0][1]),
+
1174 xt::range(m_pad[1][0], m_l.shape(1) - m_pad[1][1]),
+
1175 xt::range(m_pad[2][0], m_l.shape(2) - m_pad[2][1]));
+
+
+
1178 for (
auto& i : m_l) {
+
+
+
+
+
+
1184 array_type::tensor<double, 2> average_position(
const T& lab)
const
+
+
+
1187 size_t N = xt::amax(lab)() + 1;
+
+
+
1190 array_type::tensor<double, 2> x = xt::zeros<double>({N, size_t(3)});
+
1191 array_type::tensor<double, 1> n = xt::zeros<double>({N});
+
+
1193 for (
size_t h = 0; h < lab.shape(0); ++h) {
+
1194 for (
size_t i = 0; i < lab.shape(1); ++i) {
+
1195 for (
size_t j = 0; j < lab.shape(2); ++j) {
+
+
1197 int l = lab(h, i, j);
+
+
+
1200 x(l, 0) += (double)h;
+
1201 x(l, 1) += (double)i;
+
1202 x(l, 2) += (double)j;
+
+
+
+
+
+
+
+
1210 n = xt::where(xt::equal(n, 0), 1, n);
+
+
+
1213 for (
size_t i = 0; i < x.shape(1); ++i) {
+
1214 auto xi = xt::view(x, xt::all(), i);
+
+
+
+
+
+
+
1221 array_type::tensor<double, 2> average_position_periodic()
const
+
+
+
+
+
+
1227 auto x_np = this->average_position(m_l_np);
+
+
+
1230 auto mid = detail::half_shape(m_shape);
+
+
+
1233 array_type::tensor<double, 2> shift = xt::zeros<double>({x_np.shape(0),
size_t(3)});
+
+
+
1236 for (
size_t i = 0; i < shift.shape(0); ++i) {
+
1237 for (
size_t j = 0; j < shift.shape(1); ++j) {
+
1238 if (x_np(i, j) > mid[j]) {
+
1239 shift(i, j) = -(double)m_shape[j];
+
+
+
+
+
+
1245 size_t N = xt::amax(m_l)() + 1;
+
+
+
1248 array_type::tensor<double, 2> x = xt::zeros<double>({N, size_t(3)});
+
1249 array_type::tensor<double, 1> n = xt::zeros<double>({N});
+
+
1251 for (
size_t h = 0; h < m_l.shape(0); ++h) {
+
1252 for (
size_t i = 0; i < m_l.shape(1); ++i) {
+
1253 for (
size_t j = 0; j < m_l.shape(2); ++j) {
+
+
1255 int l = m_l_np(h, i, j);
+
+
+
1258 x(relabel(l), 0) += (double)h + shift(l, 0);
+
1259 x(relabel(l), 1) += (double)i + shift(l, 1);
+
1260 x(relabel(l), 2) += (double)j + shift(l, 2);
+
1261 n(relabel(l)) += 1.0;
+
+
+
+
+
+
+
1268 n = xt::where(xt::equal(n, 0), 1, n);
+
+
+
1271 for (
size_t i = 0; i < x.shape(1); ++i) {
+
1272 auto xi = xt::view(x, xt::all(), i);
+
+
1274 xi = xt::where(xi < 0, xi + m_shape[i], xi);
+
+
+
+
+
+
1280 static const size_t MAX_DIM = 3;
+
1281 std::vector<size_t> m_shape;
+
1282 std::vector<std::vector<size_t>> m_pad;
+
1283 array_type::tensor<int, 3> m_kernel;
+
+
1285 array_type::tensor<int, 3> m_l;
+
1286 array_type::tensor<int, 3> m_l_np;
+
+
+
+
+
+
1291template <
size_t Dimension,
bool Periodicity>
+
1292class ClusterLabellerOverload :
public ClusterLabeller<Dimension, Periodicity> {
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1318 GOOSEEYE_ASSERT(f.layout() == xt::layout_type::row_major, std::runtime_error);
+
+
1320 auto n = f.dimension();
+
1321 if (n == 1 && periodic) {
+
1322 return detail::ClusterLabellerOverload<1, true>(f).get();
+
+
1324 if (n == 1 && !periodic) {
+
1325 return detail::ClusterLabellerOverload<1, false>(f).get();
+
+
1327 if (n == 2 && periodic) {
+
1328 return detail::ClusterLabellerOverload<2, true>(f).get();
+
+
1330 if (n == 2 && !periodic) {
+
1331 return detail::ClusterLabellerOverload<2, false>(f).get();
+
+
1333 if (n == 3 && periodic) {
+
1334 return detail::ClusterLabellerOverload<3, true>(f).get();
+
+
1336 if (n == 3 && !periodic) {
+
1337 return detail::ClusterLabellerOverload<3, false>(f).get();
+
+
1339 throw std::runtime_error(
"Please use ClusterLabeller directly for dimensions > 3.");
+
+
+
+
1349template <
typename T,
typename U,
typename V>
+
1350[[deprecated(
"Will not be supported in the future. See Python warning for new API.")]]
inline T
+
+
1351pos2img(
const T& img,
const U& positions,
const V& labels)
+
+
+
1354 "i = ravel_multi_index(positions.T, img.shape); img.flat[i] = labels")
+
+
+
1357 GOOSEEYE_ASSERT(img.dimension() == positions.shape(1), std::out_of_range);
+
1358 GOOSEEYE_ASSERT(positions.shape(0) == labels.size(), std::out_of_range);
+
+
+
1361 using value_type =
typename T::value_type;
+
+
+
1364 if (res.dimension() == 1) {
+
1365 xt::view(res, xt::keep(positions)) = labels;
+
+
1367 else if (res.dimension() == 2) {
+
1368 for (
size_t i = 0; i < positions.shape(0); ++i) {
+
1369 res(positions(i, 0), positions(i, 1)) =
static_cast<value_type
>(labels(i));
+
+
+
+
1373 for (
size_t i = 0; i < positions.shape(0); ++i) {
+
1374 res(positions(i, 0), positions(i, 1), positions(i, 2)) =
+
1375 static_cast<value_type
>(labels(i));
+
+
+
+
+
+
+
+
+
+
+
+
1397 bool periodic =
true)
+
+
1399 if (positions.size() == 0) {
+
1400 return xt::zeros<double>({shape.size()});
+
+
+
+
1404 return xt::mean(positions, 0);
+
+
+
1407 double pi = xt::numeric_constants<double>::PI;
+
1408 auto theta = 2.0 * pi * positions / shape;
+
1409 auto xi = xt::cos(theta);
+
1410 auto zeta = xt::sin(theta);
+
1411 auto xi_bar = xt::mean(xi, 0);
+
1412 auto zeta_bar = xt::mean(zeta, 0);
+
1413 auto theta_bar = xt::atan2(-zeta_bar, -xi_bar) + pi;
+
1414 return shape * theta_bar / (2.0 * pi);
+
+
+
+
+
+
+
+
+
+
1426 bool periodic =
true)
+
+
1428 if (positions.size() == 0) {
+
1429 return xt::zeros<double>({shape.size()});
+
+
+
+
1433 return xt::average(positions, weights, 0);
+
+
+
1436 double pi = xt::numeric_constants<double>::PI;
+
1437 auto theta = 2.0 * pi * positions / shape;
+
1438 auto xi = xt::cos(theta);
+
1439 auto zeta = xt::sin(theta);
+
1440 auto xi_bar = xt::average(xi, weights, 0);
+
1441 auto zeta_bar = xt::average(zeta, weights, 0);
+
1442 auto theta_bar = xt::atan2(-zeta_bar, -xi_bar) + pi;
+
1443 return shape * theta_bar / (2.0 * pi);
+
+
+
+
+
+
+
+
+
+
+
+
1454 PositionList() =
default;
+
+
+
1457 void set(
const T& condition)
+
+
1459 positions = xt::from_indices(xt::argwhere(condition));
+
+
+
1462 template <
class T,
class W>
+
1463 void set(
const T& condition,
const W& w)
+
+
1465 auto pos = xt::argwhere(condition);
+
1466 weights = xt::empty<double>({pos.size()});
+
1467 for (
size_t i = 0; i < pos.size(); ++i) {
+
1468 weights(i) = w[pos[i]];
+
+
1470 positions = xt::from_indices(pos);
+
+
+
+
+
+
1484template <
class T,
class N>
+
+
+
+
+
1488 static_assert(std::is_integral<typename T::value_type>::value,
"Integral labels required.");
+
+
+
+
1492 size_t rank = labels.dimension();
+
+
+
1495 detail::PositionList plist;
+
+
1497 for (
size_t l = 0; l < names.size(); ++l) {
+
1498 plist.set(xt::equal(labels, names(l)));
+
1499 if (plist.positions.size() == 0) {
+
+
+
1502 xt::view(ret, l, xt::all()) =
center(shape, plist.positions, periodic);
+
+
+
+
+
+
+
1512template <
class T,
class W,
class N>
+
+
+
+
+
1516 static_assert(std::is_integral<typename T::value_type>::value,
"Integral labels required.");
+
1517 GOOSEEYE_ASSERT(xt::has_shape(labels, weights.shape()), std::out_of_range);
+
+
+
+
1521 size_t rank = labels.dimension();
+
+
+
1524 detail::PositionList plist;
+
+
1526 for (
size_t l = 0; l < names.size(); ++l) {
+
1527 plist.set(xt::equal(labels, names(l)), weights);
+
1528 if (plist.positions.size() == 0) {
+
+
+
1531 xt::view(ret, l, xt::all()) =
+
+
+
+
+
+
+
+
+
+
+
+
+
1561 Ensemble(
const std::vector<size_t>& roi,
bool periodic =
true,
bool variance =
true);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1630 void mean(
const T& f);
+
+
1637 template <
class T,
class M>
+
1638 void mean(
const T& f,
const M& fmask);
+
+
+
1646 void S2(
const T& f,
const T& g);
+
+
1655 template <
class T,
class M>
+
1656 void S2(
const T& f,
const T& g,
const M& fmask,
const M& gmask);
+
+
+
1664 void C2(
const T& f,
const T& g);
+
+
1673 template <
class T,
class M>
+
1674 void C2(
const T& f,
const T& g,
const M& fmask,
const M& gmask);
+
+
+
1682 void W2(
const T& w,
const T& f);
+
+
1690 template <
class T,
class M>
+
1691 void W2(
const T& w,
const T& f,
const M& fmask);
+
+
1700 template <
class C,
class T>
+
+
+
+
1712 template <
class C,
class T,
class M>
+
+
+
+
+
+
+
+
+
+
+
1734 template <
class T,
class M>
+
+
+
+
+
+
+
+
+
+
+
1750 Type m_stat = Type::Unset;
+
+
+
1753 static const size_t MAX_DIM = 3;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1770 std::vector<size_t> m_shape_orig;
+
+
+
1773 std::vector<size_t> m_shape;
+
+
+
1776 std::vector<std::vector<size_t>> m_pad;
+
+
+
+
+
+
+
+
1788inline auto distance(
const std::vector<size_t>& roi);
+
+
1796inline auto distance(
const std::vector<size_t>& roi,
size_t axis);
+
+
1804inline auto distance(
const std::vector<size_t>& roi,
const std::vector<double>& h);
+
+
1813inline auto distance(
const std::vector<size_t>& roi,
const std::vector<double>& h,
size_t axis);
+
+
+
1823inline auto S2(
const std::vector<size_t>& roi,
const T& f,
const T& g,
bool periodic =
true);
+
+
1834template <
class T,
class M>
+
+
1836S2(
const std::vector<size_t>& roi,
+
+
+
+
+
1841 bool periodic =
true);
+
+
+
1851inline auto C2(
const std::vector<size_t>& roi,
const T& f,
const T& g,
bool periodic =
true);
+
+
1862template <
class T,
class M>
+
+
1864C2(
const std::vector<size_t>& roi,
+
+
+
+
+
1869 bool periodic =
true);
+
+
+
1879inline auto W2(
const std::vector<size_t>& roi,
const T& w,
const T& f,
bool periodic =
true);
+
+
1889template <
class T,
class M>
+
+
1891W2(
const std::vector<size_t>& roi,
const T& w,
const T& f,
const M& fmask,
bool periodic =
true);
+
+
1902template <
class C,
class T>
+
+
1904W2c(
const std::vector<size_t>& roi,
+
+
+
+
+
1909 bool periodic =
true);
+
+
1921template <
class C,
class T,
class M>
+
+
1923W2c(
const std::vector<size_t>& roi,
+
+
+
+
+
+
1929 bool periodic =
true);
+
+
+
1938inline auto heightheight(
const std::vector<size_t>& roi,
const T& f,
bool periodic =
true);
+
+
1947template <
class T,
class M>
+
+
1949heightheight(
const std::vector<size_t>& roi,
const T& f,
const M& fmask,
bool periodic =
true);
+
+
+
+
1960L(
const std::vector<size_t>& roi,
+
+
1962 bool periodic =
true,
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
(Incrementally) Label clusters (0 as background, 1..n as labels).
+
void prune()
Prune: renumber labels to lowest possible label, see also AvalancheSegmenter::nlabels.
+
void add_image(const T &img)
Add image.
+
auto size() const
Size of ClusterLabeller::s and ClusterLabeller::labels (== prod(shape)).
+
ClusterLabeller(const T &shape)
+
ClusterLabeller(const T &shape, const K &kernel)
+
std::string repr() const
Basic class info.
+
const auto & shape() const
Shape of ClusterLabeller::s and ClusterLabeller::labels.
+
static constexpr bool Periodic
Periodicity of the system.
+
void add_points(const T &begin, const T &end)
Add sequence of points.
+
const auto & labels() const
Per block, the label (0 for background).
+
std::vector< size_t > add_sequence(const T &idx)
Add a sequence of points.
+
void add_points(const T &idx)
Add sequence of points.
+
void reset()
Reset labels to zero.
+
static constexpr size_t Dim
Dimensionality of the system.
+
Compute clusters and obtain certain characteristic about them.
+
array_type::array< int > centers() const
Return label only in the center of gravity (1..n).
+
array_type::tensor< size_t, 1 > sizes() const
Return size per cluster.
+
Clusters(const T &f, const S &kernel, bool periodic=true)
Constructor.
+
Clusters(const T &f, bool periodic=true)
Constructor.
+
array_type::array< int > labels() const
Return labels (1..n).
+
array_type::tensor< double, 2 > center_positions(bool as3d=false) const
Return positions of the centers of gravity (in the original rank, or as 3-d).
+
Compute ensemble averaged statistics, by repetitively calling the member-function of a certain statis...
+
void mean(const T &f)
Add realization to arithmetic mean.
+
void L(const T &f, path_mode mode=path_mode::Bresenham)
Add realization to lineal-path function.
+
void S2(const T &f, const T &g)
Add realization to 2-point correlation: P(f(i) * g(i + di)).
+
array_type::array< double > data_second() const
Get raw-data: ensemble sum of the second moment: x_1^2 + x_2^2 + ...
+
Ensemble()=default
Constructor.
+
array_type::array< double > variance() const
Get ensemble variance.
+
void heightheight(const T &f)
Add realization to height-height correlation.
+
void W2c(const C &clusters, const C ¢ers, const T &f, path_mode mode=path_mode::Bresenham)
Add realization to collapsed weighted 2-point correlation.
+
array_type::array< double > data_first() const
Get raw-data: ensemble sum of the first moment: x_1 + x_2 + ...
+
void W2(const T &w, const T &f)
Add realization to weighted 2-point correlation.
+
array_type::array< double > distance() const
Get the relative distance of each pixel in the 'region-of-interest' to its center.
+
array_type::array< double > result() const
Get ensemble average.
+
array_type::array< double > norm() const
Get raw-data: normalisation (number of measurements per pixel).
+
void C2(const T &f, const T &g)
Add realization to 2-point cluster function: P(f(i) == g(i + di)).
+
+
#define GOOSEEYE_ASSERT(expr, assertion)
All assertions are implemented as:
+
#define GOOSEEYE_WARNING_PYTHON(message)
Warnings specific to the Python API are implemented as:
+
+
+
+
xt::xtensor< T, N > tensor
Fixed (static) rank array.
+
array_type::array< int > nearest(size_t ndim)
Return kernel with nearest neighbours.
+
Toolbox to compute statistics.
+
path_mode
Different methods to compute a pixel-path.
+
@ Bresenham
Bresenham algorithm.
+
+
@ full
Similar to actual selecting every voxel that is crossed.
+
array_type::tensor< double, 2 > labels_centers_of_mass(const T &labels, const W &weights, const N &names, bool periodic=true)
Get the position of the center of each label.
+
array_type::tensor< double, 1 > center_of_mass(const array_type::tensor< double, 1 > &shape, const array_type::tensor< double, 2 > &positions, const array_type::tensor< double, 1 > &weights, bool periodic=true)
Return the geometric center of a list of positions.
+
array_type::tensor< double, 1 > center(const array_type::tensor< double, 1 > &shape, const array_type::tensor< double, 2 > &positions, bool periodic=true)
Return the geometric center of a list of positions.
+
auto L(const std::vector< size_t > &roi, const T &f, bool periodic=true, path_mode mode=path_mode::Bresenham)
Lineal-path function.
+
T dilate(const T &f, const S &kernel, const array_type::tensor< size_t, 1 > &iterations, bool periodic)
Dilate image.
+
T labels_prune(const T &labels)
Prune labels: renumber labels to lowest possible label starting from 1.
+
auto W2c(const std::vector< size_t > &roi, const C &clusters, const C ¢ers, const T &f, path_mode mode=path_mode::Bresenham, bool periodic=true)
Collapsed weighted 2-point correlation.
+
array_type::tensor< typename T::value_type, 2 > labels_sizes(const T &labels)
Size per label.
+
array_type::tensor< typename T::value_type, 2 > labels_map(const T &a, const T &b)
Get a map to relabel from a to b.
+
auto C2(const std::vector< size_t > &roi, const T &f, const T &g, bool periodic=true)
2-point cluster function: P(f(i) == g(i + di)).
+
array_type::array< int > clusters(const T &f, bool periodic=true)
Compute clusters.
+
array_type::tensor< double, 2 > labels_centers(const T &labels, const N &names, bool periodic=true)
Get the position of the center of each label.
+
L labels_rename(const L &labels, const A &rename)
Rename labels.
+
T pos2img(const T &img, const U &positions, const V &labels)
Convert positions to an image.
+
auto W2(const std::vector< size_t > &roi, const T &w, const T &f, bool periodic=true)
Weighted 2-point correlation.
+
array_type::array< int > dummy_circles(const std::vector< size_t > &shape, const array_type::tensor< int, 1 > &row, const array_type::tensor< int, 1 > &col, const array_type::tensor< int, 1 > &r, bool periodic=true)
Dummy image with circles.
+
array_type::tensor< int, 2 > path(const array_type::tensor< int, 1 > &x0, const array_type::tensor< int, 1 > &x1, path_mode mode=path_mode::Bresenham)
Compute a path between two pixels.
+
auto distance(const std::vector< size_t > &roi)
Get the relative distance of each pixel in the 'region-of-interest' to its center.
+
array_type::tensor< size_t, 1 > relabel_map(const T &a, const S &b)
Find map to relabel from a to b.
+
auto S2(const std::vector< size_t > &roi, const T &f, const T &g, bool periodic=true)
2-point correlation: P(f(i) * g(i + di)).
+
auto heightheight(const std::vector< size_t > &roi, const T &f, bool periodic=true)
Height-height correlation.
+
L labels_reorder(const L &labels, const A &order)
Reorder labels.
+
+
1.9.8
+ 
');
+}
+
+function showSyncOn(n,relpath)
+{
+ n.html('
');
+}
+
+function toggleSyncButton(relpath)
+{
+ var navSync = $('#nav-sync');
+ if (navSync.hasClass('sync')) {
+ navSync.removeClass('sync');
+ showSyncOff(navSync,relpath);
+ storeLink(stripPath2(pathName())+hashUrl());
+ } else {
+ navSync.addClass('sync');
+ showSyncOn(navSync,relpath);
+ deleteLink();
+ }
+}
+
+var loadTriggered = false;
+var readyTriggered = false;
+var loadObject,loadToRoot,loadUrl,loadRelPath;
+
+$(window).on('load',function(){
+ if (readyTriggered) { // ready first
+ navTo(loadObject,loadToRoot,loadUrl,loadRelPath);
+ showRoot();
+ }
+ loadTriggered=true;
+});
+
+function initNavTree(toroot,relpath)
+{
+ var o = new Object();
+ o.toroot = toroot;
+ o.node = new Object();
+ o.node.li = document.getElementById("nav-tree-contents");
+ o.node.childrenData = NAVTREE;
+ o.node.children = new Array();
+ o.node.childrenUL = document.createElement("ul");
+ o.node.getChildrenUL = function() { return o.node.childrenUL; };
+ o.node.li.appendChild(o.node.childrenUL);
+ o.node.depth = 0;
+ o.node.relpath = relpath;
+ o.node.expanded = false;
+ o.node.isLast = true;
+ o.node.plus_img = document.createElement("span");
+ o.node.plus_img.className = 'arrow';
+ o.node.plus_img.innerHTML = arrowRight;
+
+ if (localStorageSupported()) {
+ var navSync = $('#nav-sync');
+ if (cachedLink()) {
+ showSyncOff(navSync,relpath);
+ navSync.removeClass('sync');
+ } else {
+ showSyncOn(navSync,relpath);
+ }
+ navSync.click(function(){ toggleSyncButton(relpath); });
+ }
+
+ if (loadTriggered) { // load before ready
+ navTo(o,toroot,hashUrl(),relpath);
+ showRoot();
+ } else { // ready before load
+ loadObject = o;
+ loadToRoot = toroot;
+ loadUrl = hashUrl();
+ loadRelPath = relpath;
+ readyTriggered=true;
+ }
+
+ $(window).bind('hashchange', function(){
+ if (window.location.hash && window.location.hash.length>1){
+ var a;
+ if ($(location).attr('hash')){
+ var clslink=stripPath(pathName())+':'+hashValue();
+ a=$('.item a[class$="'+clslink.replace(/1|%O$WD@{VPM$7~Ar*{o?;hlAFyLXmaDC0y
znK1_#cQqJWPES%4Uujug^TE?jMft$}Eq^WaR~)%f)vSNs&gek&x%A9X9sM