op_dimensionality_estimate

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<op_repeat_slice>));
 		m.insert(std::make_pair(std::string("component_foreach"), &instantiate<op_component_foreach>));
 		m.insert(std::make_pair(std::string("normal_estimate"), &instantiate<op_normal_estimate>));
+		m.insert(std::make_pair(std::string("dimensionality_estimate"), &instantiate<op_dimensionality_estimate>));
 		m.insert(std::make_pair(std::string("segment"), &instantiate<op_segment>));
 		m.insert(std::make_pair(std::string("keep_neighbours"), &instantiate<op_keep_neighbours>));
 		m.insert(std::make_pair(std::string("free"), &instantiate<op_free>));

voxec.h

@@ -2496,6 +2496,189 @@ class op_normal_estimate : public voxel_operation {
 	}
 };
 
+
+class op_dimensionality_estimate : public voxel_operation {
+public:
+	const std::vector<argument_spec>& arg_names() const {
+		static std::vector<argument_spec> 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<int16_t> v;
+
+		auto voxels = (regular_voxel_storage*)scope.get_value<abstract_voxel_storage*>("input");
+		auto result = voxels->empty_copy_as(&fmt);
+
+		boost::optional<int> max_depth, max_depth_2, neighbourhood_size;
+		boost::optional<double> distance_2;
+
+		try {
+			max_depth = scope.get_value<int>("max_depth");
+		} catch (scope_map::not_in_scope&) {}
+		try {
+			max_depth_2 = scope.get_value<int>("max_depth_2");
+		} catch (scope_map::not_in_scope&) {}
+		try {
+			neighbourhood_size = scope.get_value<int>("neighbourhood_size");
+		} catch (scope_map::not_in_scope&) {}
+		try {
+			distance_2 = scope.get_value<double>("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<int>("reposition", 0) ? 2 : 1;
+
+		std::vector<float> 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<Eigen::Matrix<float, 3, Eigen::Dynamic>>(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<double> 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<double>::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<Eigen::MatrixXf> eig(cov);
+
+				std::array<double, 3> 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<int16_t>::max() - 1)),
+					(int16_t)(norm_eigv[1] * (float)(std::numeric_limits<int16_t>::max() - 1)),
+					(int16_t)(norm_eigv[2] * (float)(std::numeric_limits<int16_t>::max() - 1)),
+					0
+				};
+
+				result->Set(*it, &v);
+			}
+		}
+
+		return result;
+	}
+};
+
 namespace {
 	class check_curvature_and_normal_deviation {
 	private:

wrap/wrapper.i

@@ -47,7 +47,7 @@
 	PyObject* get(long i, long j, long k) const {
 		auto acvs_voxels = dynamic_cast<abstract_chunked_voxel_storage const*>($self);
 		if (!acvs_voxels) {
-			throw std::runtime_error("Unsupported");
+			throw std::runtime_error("Unsupported storage");
 		}
 		auto ijk_long = (make_vec<long>(i, j, k) - (acvs_voxels->grid_offset() * acvs_voxels->chunk_size()).as<long>());
 		auto ijk = ijk_long.as<size_t>();
@@ -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 {