From 52c54ed5305ca69b154607ad3da42702f2897601 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Sat, 24 Feb 2024 21:00:43 +0100 Subject: [PATCH] op_dimensionality_estimate --- 3rdparty/eigen | 2 +- voxec.cpp | 1 + voxec.h | 183 +++++++++++++++++++++++++++++++++++++++++++++++++ wrap/wrapper.i | 8 ++- 4 files changed, 191 insertions(+), 3 deletions(-) diff --git a/3rdparty/eigen b/3rdparty/eigen index 0fd6b4f..3147391 160000 --- a/3rdparty/eigen +++ b/3rdparty/eigen @@ -1 +1 @@ -Subproject commit 0fd6b4f71dd85b2009ee4d1aeb296e2c11fc9d68 +Subproject commit 3147391d946bb4b6c68edd901f2add6ac1f31f8c diff --git a/voxec.cpp b/voxec.cpp index 9b655c0..15b46c9 100644 --- a/voxec.cpp +++ b/voxec.cpp @@ -65,6 +65,7 @@ voxel_operation_map::map_t& voxel_operation_map::map() { m.insert(std::make_pair(std::string("repeat_slice"), &instantiate)); m.insert(std::make_pair(std::string("component_foreach"), &instantiate)); m.insert(std::make_pair(std::string("normal_estimate"), &instantiate)); + m.insert(std::make_pair(std::string("dimensionality_estimate"), &instantiate)); m.insert(std::make_pair(std::string("segment"), &instantiate)); m.insert(std::make_pair(std::string("keep_neighbours"), &instantiate)); m.insert(std::make_pair(std::string("free"), &instantiate)); diff --git a/voxec.h b/voxec.h index e3ce9a8..0a678f2 100644 --- a/voxec.h +++ b/voxec.h @@ -2496,6 +2496,189 @@ class op_normal_estimate : public voxel_operation { } }; + +class op_dimensionality_estimate : public voxel_operation { +public: + const std::vector& arg_names() const { + static std::vector nm_ = { { true, "input", "voxels" }, { false, "max_depth", "integer"}, { false, "max_depth_2", "integer"}, {false, "distance_2", "real"}, {false, "neighbourhood_size", "integer"}, {false, "reposition", "integer"}}; + return nm_; + } + symbol_value invoke(const scope_map& scope) const { + static normal_and_curvature_t fmt; + normal_and_curvature v; + + auto voxels = (regular_voxel_storage*)scope.get_value("input"); + auto result = voxels->empty_copy_as(&fmt); + + boost::optional max_depth, max_depth_2, neighbourhood_size; + boost::optional distance_2; + + try { + max_depth = scope.get_value("max_depth"); + } catch (scope_map::not_in_scope&) {} + try { + max_depth_2 = scope.get_value("max_depth_2"); + } catch (scope_map::not_in_scope&) {} + try { + neighbourhood_size = scope.get_value("neighbourhood_size"); + } catch (scope_map::not_in_scope&) {} + try { + distance_2 = scope.get_value("distance_2"); + } catch (scope_map::not_in_scope&) {} + + if (!max_depth && !neighbourhood_size) { + throw std::runtime_error("Either max_depth or neighbourhood_size needs to be provided"); + } + + int reposition = scope.get_value_or("reposition", 0) ? 2 : 1; + + std::vector coords; + + if (neighbourhood_size) { + coords.reserve(*neighbourhood_size); + } else if (max_depth) { + auto box_dim = 1 + (*max_depth) * 2; + coords.reserve(box_dim * box_dim * box_dim); + } + + for (auto it = voxels->begin(); it != voxels->end(); ++it) { + visitor<26> vis; + + auto seed = *it; + + for (int iter = 0; iter < reposition; ++iter) { + + int previous_coord_count = -1; + int max_depth_guess = ((int)std::ceil(std::pow(26, 1. / 3) - 1)) / 2; + + while (true) { + if (neighbourhood_size) { + vis.max_depth = max_depth_guess; + } else { + vis.max_depth = iter == 1 && max_depth_2 ? max_depth_2 : max_depth; + } + + coords.clear(); + + vis([&coords, neighbourhood_size](const tagged_index& pos) { + if (pos.which == tagged_index::VOXEL) { + if (!neighbourhood_size || (coords.size() / 3 < *neighbourhood_size)) { + coords.push_back(pos.pos.get(0)); + coords.push_back(pos.pos.get(1)); + coords.push_back(pos.pos.get(2)); + } + } else { + throw std::runtime_error("Unexpected"); + } + }, voxels, seed); + + if (neighbourhood_size) { + if (coords.size() / 3 == *neighbourhood_size) { + break; + } + if (previous_coord_count == (int)coords.size()) { + break; + } + if (max_depth) { + if (*max_depth == max_depth_guess) { + break; + } + } + previous_coord_count = (int)coords.size(); + } else if (max_depth) { + break; + } + + ++max_depth_guess; + } + + auto points = Eigen::Map>(coords.data(), 3, coords.size() / 3).transpose(); + + if (reposition == 2 && iter == 0) { +#if 1 + auto mean = points.colwise().mean(); +#elif 0 + // median + Eigen::RowVector3f mean; + size_t i = 0; + for (const auto& c : points.colwise()) { + std::vector r; + r.reserve(c.size()); + for (auto& v : c) { + r.push_back(v); + } + std::sort(r.begin(), r.end()); + auto v = r.size() % 2 == 0 + ? (r[(r.size() - 2) / 2] + r[(r.size() - 2) / 2 + 1]) / 2. + : r[r.size() / 2]; + mean(i++) = v; + } +#else + // bounding box center + auto mi = points.colwise().minCoeff(); + auto ma = points.colwise().maxCoeff(); + auto mean = (mi + ma) / 2.; +#endif + + // Find point in neighborhood closest to median + auto minl = std::numeric_limits::infinity(); + Eigen::RowVectorXf p; + for (const auto& r : points.rowwise()) { + auto l = (r - mean).norm(); + if (l < minl) { + minl = l; + p = r; + } + } + + // Assign to seed to reposition neighborhood around found median + seed.get<0>() = (size_t)p.data()[0]; + seed.get<1>() = (size_t)p.data()[1]; + seed.get<2>() = (size_t)p.data()[2]; + + // Continue within outer loop to restart traversal from new seed + continue; + } + + Eigen::MatrixXf centered = points.rowwise() - points.colwise().mean(); + + if (distance_2) { + auto norms = centered.rowwise().norm().array(); + auto mask = norms < (*distance_2); + Eigen::MatrixXf filtered(mask.count(), centered.cols()); + size_t idx = 0; + for (size_t i = 0; i < mask.size(); ++i) { + if (mask(i)) { + filtered.row(idx++) = centered.row(i); + } + } + centered = filtered; + } + + Eigen::MatrixXf cov = centered.adjoint() * centered; + Eigen::SelfAdjointEigenSolver eig(cov); + + std::array norm_eigv = { + eig.eigenvalues().data()[0] / eig.eigenvalues().norm(), + eig.eigenvalues().data()[1] / eig.eigenvalues().norm(), + eig.eigenvalues().data()[2] / eig.eigenvalues().norm() + }; + + v.nxyz_curv = { + (int16_t)(norm_eigv[0] * (float)(std::numeric_limits::max() - 1)), + (int16_t)(norm_eigv[1] * (float)(std::numeric_limits::max() - 1)), + (int16_t)(norm_eigv[2] * (float)(std::numeric_limits::max() - 1)), + 0 + }; + + result->Set(*it, &v); + } + } + + return result; + } +}; + namespace { class check_curvature_and_normal_deviation { private: diff --git a/wrap/wrapper.i b/wrap/wrapper.i index a9e24ff..882f2c1 100644 --- a/wrap/wrapper.i +++ b/wrap/wrapper.i @@ -47,7 +47,7 @@ PyObject* get(long i, long j, long k) const { auto acvs_voxels = dynamic_cast($self); if (!acvs_voxels) { - throw std::runtime_error("Unsupported"); + throw std::runtime_error("Unsupported storage"); } auto ijk_long = (make_vec(i, j, k) - (acvs_voxels->grid_offset() * acvs_voxels->chunk_size()).as()); auto ijk = ijk_long.as(); @@ -60,6 +60,10 @@ } if ($self->value_bits() == 1) { return PyBool_FromLong(self->Get(ijk)); + } else if ($self->value_bits() == 8) { + uint8_t v; + $self->Get(ijk, &v); + return PyLong_FromLong(v); } else if ($self->value_bits() == 32) { uint32_t v; $self->Get(ijk, &v); @@ -78,7 +82,7 @@ return tup; } } else { - throw std::runtime_error("Unsupported"); + throw std::runtime_error("Unsupported data type, size: " + std::to_string($self->value_bits())); } } PyObject* get(double x, double y, double z) const {