diff --git a/docs/structural_lva.md b/docs/structural_lva.md new file mode 100644 index 000000000..99a4eae96 --- /dev/null +++ b/docs/structural_lva.md @@ -0,0 +1,103 @@ +# Structural LVA prototype + +The `feature/structural-lva` branch contains a native Polatory structural-LVA path. + +## Python API + +```python +import numpy as np +import polatory +from polatory import three as p3 + +rbf = p3.CovSpheroidal3([10.0, 100.0]) +model = p3.Model(rbf, 0) + +trend_input = polatory.StructuralTrendInput3( + vertices=mesh_vertices, + faces=mesh_faces.astype(np.int64), + strength=5.0, + range=50.0, +) + +structural = polatory.StructuralInterpolant3( + model, + outside_value=-1.0, + blend_power=7.0, +) + +domains = structural.fit_from_meshes( + points, + values, + [trend_input], + tolerance=1e-6, + trend_type=polatory.StructuralTrendType.STRONGEST_ALONG_INPUTS, +) + +predictions = structural.evaluate(query_points) +``` + +`fit_from_meshes` automatically: + +1. calculates equal-weight incident triangle normals at mesh vertices; +2. finds the nearest mesh vertex to each structural sample location; +3. applies the recovered single-input decay law; +4. constructs determinant-one local anisotropy matrices; +5. divides the interpolation data into overlapping spatial domains; +6. assigns every data point inside each expanded domain box as local support; +7. fits and blends the local Polatory interpolants. + +## Recovered single-input field + +For the nearest mesh-vertex normal `n`, distance `d`, input `strength`, and input `range`: + +```text +q = exp(-d / range) +r = 1 + (strength - 1) * q +M = r^(-1/3) * (I - n n^T) + r^(2/3) * (n n^T) +``` + +The equal-weight vertex normal and this matrix equation were verified against the supplied Leapfrog project. + +## Trend types + +The public API exposes Leapfrog-style names: + +- `STRONGEST_ALONG_INPUTS` +- `BLENDING` +- `NON_DECAYING` + +Multiple mesh inputs are accepted. `STRONGEST_ALONG_INPUTS` selects the input with the largest local decayed anisotropy contribution. + +`BLENDING` currently uses an axial weighted-normal blend. This is an experimental implementation. The Leapfrog manual confirms that multiple inputs are blended according to individual strength and that the result decays away from the meshes, but it does not publish the exact vector/tensor combination rule. + +## Parameters suitable for an RSGeo UI + +User-facing structural trend settings: + +- name; +- trend type; +- one or more mesh inputs; +- per-input strength; +- per-input range for decaying modes; +- optional global mean trend in a later parity update; +- compatibility mode in a later parity update. + +Advanced internal settings that should normally stay hidden: + +- domain size; +- domain overlap; +- minimum local support count; +- local blend power. + +## Known parity limits + +The following still require controlled Leapfrog A/B tests before claiming exact parity: + +1. multiple-input `BLENDING` orientation and strength equations; +2. transition behaviour where two inputs have equal influence in `STRONGEST_ALONG_INPUTS`; +3. global mean trend interaction; +4. Version 1 versus Version 2 compatibility; +5. Leapfrog's internal automatic domain decomposition; +6. whether one tuned local blend power generalises across unrelated datasets. + +The automatic domain builder is therefore ready for compilation and single-mesh validation, while multi-mesh blending remains explicitly experimental. diff --git a/include/polatory/isosurface/structural_rbf_field_function.hpp b/include/polatory/isosurface/structural_rbf_field_function.hpp new file mode 100644 index 000000000..d17343844 --- /dev/null +++ b/include/polatory/isosurface/structural_rbf_field_function.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace polatory::isosurface { + +class StructuralRbfFieldFunction : public FieldFunction { + static constexpr double kInfinity = std::numeric_limits::infinity(); + + public: + explicit StructuralRbfFieldFunction(structural::StructuralInterpolant3& interpolant, + double accuracy = kInfinity) + : interpolant_(interpolant), accuracy_(accuracy) {} + + VecX operator()(const geometry::Points3& points) const override { + return interpolant_.evaluate_impl(points); + } + + void set_evaluation_bbox(const geometry::Bbox3& bbox) override { + interpolant_.set_evaluation_bbox_impl(bbox, accuracy_); + } + + private: + structural::StructuralInterpolant3& interpolant_; + double accuracy_; +}; + +} // namespace polatory::isosurface diff --git a/include/polatory/point_cloud/sdf_data_generator.hpp b/include/polatory/point_cloud/sdf_data_generator.hpp index 61b1fd3b0..8cc985e58 100644 --- a/include/polatory/point_cloud/sdf_data_generator.hpp +++ b/include/polatory/point_cloud/sdf_data_generator.hpp @@ -15,14 +15,19 @@ class SdfDataGenerator { SdfDataGenerator(const geometry::Points3& points, const geometry::Vectors3& normals, double offset, const Mat3& aniso); + // Compatibility overload for the existing Python binding. The current C++ + // implementation uses a single offset; max_distance is used as that offset. + SdfDataGenerator(const geometry::Points3& points, const geometry::Vectors3& normals, + double min_distance, double max_distance, const Mat3& aniso); + const geometry::Points3& sdf_points() const; const VecX& sdf_values() const; private: static std::pair estimate_impl(const geometry::Points3& points, - const geometry::Vectors3& normals, - double offset); + const geometry::Vectors3& normals, + double offset); geometry::Points3 sdf_points_; VecX sdf_values_; diff --git a/include/polatory/polatory.hpp b/include/polatory/polatory.hpp index 4bf3717cd..010ca2223 100644 --- a/include/polatory/polatory.hpp +++ b/include/polatory/polatory.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -31,5 +32,8 @@ #include #include #include +#include +#include +#include #include #include diff --git a/include/polatory/structural/adaptive_domain_builder.hpp b/include/polatory/structural/adaptive_domain_builder.hpp new file mode 100644 index 000000000..d4c912dca --- /dev/null +++ b/include/polatory/structural/adaptive_domain_builder.hpp @@ -0,0 +1,650 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +// Experimental adaptive builder. It keeps the validated mesh sampling, +// anisotropy and overlap rules, but replaces the regular grid with recursive, +// data-balanced spatial partitions. Cells are subdivided when they are too large +// or when their sign-invariant structural normals are not sufficiently coherent. +class AdaptiveStructuralDomainBuilder3 { + using Point = geometry::Point3; + using Points = geometry::Points3; + + public: + explicit AdaptiveStructuralDomainBuilder3( + double overlap = 0.0, + double orientation_consistency = 0.97, + double minimum_core_size = 0.0, + double maximum_core_size = 0.0, + Index minimum_core_points = 24, + Index minimum_support_points = 4, + int maximum_depth = 20) + : overlap_(overlap), + orientation_consistency_(orientation_consistency), + minimum_core_size_(minimum_core_size), + maximum_core_size_(maximum_core_size), + minimum_core_points_(minimum_core_points), + minimum_support_points_(minimum_support_points), + maximum_depth_(maximum_depth) { + if (!(overlap_ >= 0.0)) { + throw std::invalid_argument("overlap must be non-negative"); + } + if (!(orientation_consistency_ > 0.0 && + orientation_consistency_ <= 1.0)) { + throw std::invalid_argument( + "orientation_consistency must be in (0, 1]"); + } + if (!(minimum_core_size_ >= 0.0)) { + throw std::invalid_argument("minimum_core_size must be non-negative"); + } + if (!(maximum_core_size_ >= 0.0)) { + throw std::invalid_argument("maximum_core_size must be non-negative"); + } + if (minimum_core_points_ <= 0) { + throw std::invalid_argument("minimum_core_points must be positive"); + } + if (minimum_support_points_ <= 0) { + throw std::invalid_argument("minimum_support_points must be positive"); + } + if (maximum_depth_ <= 0) { + throw std::invalid_argument("maximum_depth must be positive"); + } + } + + std::vector build( + const Points& points, + const std::vector& inputs, + StructuralTrendType trend_type = + StructuralTrendType::kStrongestAlongInputs, + const std::vector& model_parameters = {}) const { + if (points.rows() == 0) { + throw std::invalid_argument("points must not be empty"); + } + if (inputs.empty()) { + throw std::invalid_argument("structural trend inputs must not be empty"); + } + + auto prepared = prepare_inputs(inputs); + + auto maximum_range = 0.0; + for (const auto& input : inputs) { + maximum_range = std::max(maximum_range, input.range()); + } + + auto resolved_overlap = overlap_ > 0.0 ? overlap_ : maximum_range; + auto resolved_minimum_core_size = + minimum_core_size_ > 0.0 ? minimum_core_size_ : 0.75 * maximum_range; + auto resolved_maximum_core_size = + maximum_core_size_ > 0.0 ? maximum_core_size_ : 2.10 * maximum_range; + + if (!(resolved_maximum_core_size > 0.0)) { + resolved_maximum_core_size = 1.0; + } + if (!(resolved_minimum_core_size > 0.0)) { + resolved_minimum_core_size = 0.25 * resolved_maximum_core_size; + } + if (resolved_minimum_core_size > resolved_maximum_core_size) { + throw std::invalid_argument( + "minimum_core_size must not exceed maximum_core_size"); + } + + Points sampled_normals(points.rows(), 3); + for (Index point_i = 0; point_i < points.rows(); ++point_i) { + sampled_normals.row(point_i) = + evaluate_one(points.row(point_i), inputs, prepared, trend_type).normal; + } + + auto active_axes = detect_active_axes(points, sampled_normals); + + Node root; + root.minimum = points.colwise().minCoeff(); + root.maximum = points.colwise().maxCoeff(); + root.indices.resize(static_cast(points.rows())); + std::iota(root.indices.begin(), root.indices.end(), Index{0}); + + std::vector leaves; + split_recursive(root, points, sampled_normals, active_axes, + resolved_minimum_core_size, + resolved_maximum_core_size, 0, leaves); + + std::vector domains; + domains.reserve(leaves.size()); + + for (const auto& leaf : leaves) { + Point sample_point = Point::Zero(); + for (auto point_i : leaf.indices) { + sample_point += points.row(point_i); + } + sample_point /= static_cast(leaf.indices.size()); + + auto trend = evaluate_one(sample_point, inputs, prepared, trend_type); + + Point bbox_minimum = leaf.minimum.array() - resolved_overlap; + Point bbox_maximum = leaf.maximum.array() + resolved_overlap; + auto support_indices = + points_inside(points, bbox_minimum, bbox_maximum); + + if (static_cast(support_indices.size()) < + minimum_support_points_) { + expand_to_minimum_support(points, sample_point, + resolved_maximum_core_size, + support_indices, bbox_minimum, + bbox_maximum); + } + + domains.emplace_back(trend.anisotropy, bbox_minimum, bbox_maximum, + std::move(support_indices), model_parameters); + } + + return domains; + } + + double overlap() const { return overlap_; } + double orientation_consistency() const { + return orientation_consistency_; + } + double minimum_core_size() const { return minimum_core_size_; } + double maximum_core_size() const { return maximum_core_size_; } + Index minimum_core_points() const { return minimum_core_points_; } + Index minimum_support_points() const { return minimum_support_points_; } + int maximum_depth() const { return maximum_depth_; } + + private: + struct PreparedInput { + Points vertex_normals; + std::unique_ptr> tree; + }; + + struct TrendValue { + Point normal; + Mat3 anisotropy; + double ratio; + double distance; + Index dominant_input; + }; + + struct Node { + Point minimum; + Point maximum; + std::vector indices; + }; + + static Mat3 anisotropy_from_normal(const Point& normal, double ratio) { + Eigen::Vector3d n = normal.transpose(); + n.normalize(); + + auto tangent_scale = std::pow(ratio, -1.0 / 3.0); + auto normal_scale = std::pow(ratio, 2.0 / 3.0); + + Mat3 projector = n * n.transpose(); + return tangent_scale * (Mat3::Identity() - projector) + + normal_scale * projector; + } + + static Points compute_vertex_normals(const StructuralTrendInput3& input) { + Points normals = Points::Zero(input.vertices().rows(), 3); + + for (Index face_i = 0; face_i < input.faces().rows(); ++face_i) { + auto a = input.faces()(face_i, 0); + auto b = input.faces()(face_i, 1); + auto c = input.faces()(face_i, 2); + + Eigen::Vector3d va = input.vertices().row(a).transpose(); + Eigen::Vector3d vb = input.vertices().row(b).transpose(); + Eigen::Vector3d vc = input.vertices().row(c).transpose(); + + Eigen::Vector3d face_normal = (vb - va).cross(vc - va); + auto norm = face_normal.norm(); + if (!(norm > 0.0)) { + continue; + } + face_normal /= norm; + + normals.row(a) += face_normal.transpose(); + normals.row(b) += face_normal.transpose(); + normals.row(c) += face_normal.transpose(); + } + + for (Index vertex_i = 0; vertex_i < normals.rows(); ++vertex_i) { + auto norm = normals.row(vertex_i).norm(); + if (norm > 0.0) { + normals.row(vertex_i) /= norm; + } else { + normals.row(vertex_i) << 0.0, 0.0, 1.0; + } + } + + return normals; + } + + static std::vector prepare_inputs( + const std::vector& inputs) { + std::vector prepared; + prepared.reserve(inputs.size()); + + for (const auto& input : inputs) { + PreparedInput item; + item.vertex_normals = compute_vertex_normals(input); + item.tree = + std::make_unique>(input.vertices()); + prepared.push_back(std::move(item)); + } + + return prepared; + } + + TrendValue evaluate_one( + const Point& point, + const std::vector& inputs, + const std::vector& prepared, + StructuralTrendType trend_type) const { + std::vector normals(inputs.size()); + std::vector distances(inputs.size()); + std::vector q_values(inputs.size()); + std::vector contributions(inputs.size()); + + Index dominant_input = 0; + auto best_contribution = -1.0; + auto best_distance = std::numeric_limits::infinity(); + + for (Index input_i = 0; input_i < static_cast(inputs.size()); + ++input_i) { + std::vector indices; + std::vector nearest_distances; + prepared.at(static_cast(input_i)) + .tree->knn_search(point, 1, indices, nearest_distances); + + auto vertex_i = indices.at(0); + auto distance = nearest_distances.at(0); + auto normal = prepared.at(static_cast(input_i)) + .vertex_normals.row(vertex_i); + + auto q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-distance / + inputs.at(static_cast(input_i)) + .range()); + auto contribution = + (inputs.at(static_cast(input_i)).strength() - 1.0) * q; + + normals.at(static_cast(input_i)) = normal; + distances.at(static_cast(input_i)) = distance; + q_values.at(static_cast(input_i)) = q; + contributions.at(static_cast(input_i)) = contribution; + + if (contribution > best_contribution || + (contribution == best_contribution && distance < best_distance)) { + best_contribution = contribution; + best_distance = distance; + dominant_input = input_i; + } + } + + Point normal = normals.at(static_cast(dominant_input)); + auto ratio = + 1.0 + contributions.at(static_cast(dominant_input)); + + if (trend_type == StructuralTrendType::kBlending && inputs.size() > 1) { + Point blended = Point::Zero(); + auto weight_sum = 0.0; + auto q_sum = 0.0; + auto reference = normal; + + for (Index input_i = 0; + input_i < static_cast(inputs.size()); ++input_i) { + auto aligned = normals.at(static_cast(input_i)); + if (aligned.dot(reference) < 0.0) { + aligned *= -1.0; + } + + auto weight = contributions.at(static_cast(input_i)); + if (!(weight > 0.0)) { + weight = q_values.at(static_cast(input_i)); + } + + blended += weight * aligned; + weight_sum += weight; + q_sum += q_values.at(static_cast(input_i)); + } + + if (weight_sum > 0.0 && blended.norm() > 0.0) { + normal = blended / blended.norm(); + } + + auto contribution_sum = 0.0; + for (auto contribution : contributions) { + contribution_sum += contribution; + } + ratio = 1.0 + contribution_sum / std::max(1.0, q_sum); + } + + normal.normalize(); + auto minimum_distance = + *std::min_element(distances.begin(), distances.end()); + + return TrendValue{normal, + anisotropy_from_normal(normal, ratio), + ratio, + minimum_distance, + dominant_input}; + } + + static double axial_consistency(const std::vector& indices, + const Points& normals) { + if (indices.empty()) { + return 1.0; + } + + Mat3 tensor = Mat3::Zero(); + for (auto point_i : indices) { + Eigen::Vector3d normal = normals.row(point_i).transpose(); + tensor += normal * normal.transpose(); + } + tensor /= static_cast(indices.size()); + + Eigen::SelfAdjointEigenSolver solver(tensor); + if (solver.info() != Eigen::Success) { + return 1.0; + } + + auto trace = solver.eigenvalues().sum(); + if (!(trace > 0.0)) { + return 1.0; + } + + return solver.eigenvalues()(2) / trace; + } + + static std::array detect_active_axes( + const Points& points, const Points& normals) { + constexpr Index kNumBins = 8; + constexpr double kVariationThreshold = 0.02; + + std::array active_axes{false, false, false}; + + Mat3 global_tensor = Mat3::Zero(); + for (Index point_i = 0; point_i < normals.rows(); ++point_i) { + Eigen::Vector3d normal = normals.row(point_i).transpose(); + global_tensor += normal * normal.transpose(); + } + global_tensor /= static_cast(normals.rows()); + + for (Index axis = 0; axis < 3; ++axis) { + auto coordinate_minimum = points.col(axis).minCoeff(); + auto coordinate_maximum = points.col(axis).maxCoeff(); + auto span = coordinate_maximum - coordinate_minimum; + if (!(span > 0.0)) { + continue; + } + + std::array bin_tensors; + std::array bin_counts{}; + for (auto& tensor : bin_tensors) { + tensor.setZero(); + } + + for (Index point_i = 0; point_i < points.rows(); ++point_i) { + auto normalized = + (points(point_i, axis) - coordinate_minimum) / span; + auto bin_i = static_cast( + std::floor(normalized * static_cast(kNumBins))); + bin_i = std::clamp(bin_i, 0, kNumBins - 1); + + Eigen::Vector3d normal = normals.row(point_i).transpose(); + bin_tensors.at(static_cast(bin_i)) += + normal * normal.transpose(); + ++bin_counts.at(static_cast(bin_i)); + } + + auto maximum_variation = 0.0; + for (Index bin_i = 0; bin_i < kNumBins; ++bin_i) { + auto count = bin_counts.at(static_cast(bin_i)); + if (count == 0) { + continue; + } + auto local_tensor = + bin_tensors.at(static_cast(bin_i)) / + static_cast(count); + maximum_variation = + std::max(maximum_variation, + (local_tensor - global_tensor).norm()); + } + + active_axes.at(static_cast(axis)) = + maximum_variation > kVariationThreshold; + } + + return active_axes; + } + + static bool split_indices(const Node& node, Index axis, + const Points& points, Node& left, + Node& right) { + if (node.indices.size() < 2) { + return false; + } + + std::vector coordinates; + coordinates.reserve(node.indices.size()); + for (auto point_i : node.indices) { + coordinates.push_back(points(point_i, axis)); + } + + auto middle = coordinates.begin() + + static_cast(coordinates.size() / 2); + std::nth_element(coordinates.begin(), middle, coordinates.end()); + auto split_value = *middle; + + left.minimum = node.minimum; + left.maximum = node.maximum; + right.minimum = node.minimum; + right.maximum = node.maximum; + left.maximum(axis) = split_value; + right.minimum(axis) = split_value; + + left.indices.reserve(node.indices.size() / 2 + 1); + right.indices.reserve(node.indices.size() / 2 + 1); + + for (auto point_i : node.indices) { + if (points(point_i, axis) < split_value) { + left.indices.push_back(point_i); + } else { + right.indices.push_back(point_i); + } + } + + if (left.indices.empty() || right.indices.empty()) { + return false; + } + + return true; + } + + Index choose_split_axis( + const Node& node, const Points& points, const Points& normals, + const std::array& active_axes, bool oversized) const { + auto parent_consistency = axial_consistency(node.indices, normals); + auto best_score = -std::numeric_limits::infinity(); + Index best_axis = -1; + + for (Index axis = 0; axis < 3; ++axis) { + if (!active_axes.at(static_cast(axis))) { + continue; + } + + auto span = node.maximum(axis) - node.minimum(axis); + if (!(span > 0.0)) { + continue; + } + + Node left; + Node right; + if (!split_indices(node, axis, points, left, right)) { + continue; + } + if (static_cast(left.indices.size()) < minimum_core_points_ || + static_cast(right.indices.size()) < minimum_core_points_) { + continue; + } + + auto left_consistency = axial_consistency(left.indices, normals); + auto right_consistency = axial_consistency(right.indices, normals); + auto weighted_child_consistency = + (left_consistency * static_cast(left.indices.size()) + + right_consistency * static_cast(right.indices.size())) / + static_cast(node.indices.size()); + + auto improvement = weighted_child_consistency - parent_consistency; + auto score = oversized ? span + improvement * span : improvement * span; + + if (score > best_score) { + best_score = score; + best_axis = axis; + } + } + + return best_axis; + } + + void split_recursive( + const Node& node, const Points& points, const Points& normals, + const std::array& active_axes, + double minimum_core_size, double maximum_core_size, + int depth, std::vector& leaves) const { + auto maximum_active_span = 0.0; + for (Index axis = 0; axis < 3; ++axis) { + if (active_axes.at(static_cast(axis))) { + maximum_active_span = + std::max(maximum_active_span, + node.maximum(axis) - node.minimum(axis)); + } + } + + auto consistency = axial_consistency(node.indices, normals); + auto oversized = maximum_active_span > maximum_core_size; + auto inconsistent = + consistency < orientation_consistency_ && + maximum_active_span > minimum_core_size; + auto enough_points = + static_cast(node.indices.size()) >= + 2 * minimum_core_points_; + + if (depth >= maximum_depth_ || !enough_points || + (!oversized && !inconsistent)) { + leaves.push_back(node); + return; + } + + auto split_axis = + choose_split_axis(node, points, normals, active_axes, oversized); + if (split_axis < 0) { + leaves.push_back(node); + return; + } + + Node left; + Node right; + if (!split_indices(node, split_axis, points, left, right)) { + leaves.push_back(node); + return; + } + + split_recursive(left, points, normals, active_axes, + minimum_core_size, maximum_core_size, + depth + 1, leaves); + split_recursive(right, points, normals, active_axes, + minimum_core_size, maximum_core_size, + depth + 1, leaves); + } + + static std::vector points_inside(const Points& points, + const Point& bbox_minimum, + const Point& bbox_maximum) { + std::vector result; + result.reserve(static_cast(points.rows())); + + for (Index point_i = 0; point_i < points.rows(); ++point_i) { + auto inside = true; + for (Index axis = 0; axis < 3; ++axis) { + if (points(point_i, axis) < bbox_minimum(axis) || + points(point_i, axis) > bbox_maximum(axis)) { + inside = false; + break; + } + } + if (inside) { + result.push_back(point_i); + } + } + + return result; + } + + void expand_to_minimum_support( + const Points& points, const Point& sample_point, + double reference_size, std::vector& support_indices, + Point& bbox_minimum, Point& bbox_maximum) const { + auto expansion = 0.5 * reference_size; + + for (int attempt = 0; + attempt < 8 && + static_cast(support_indices.size()) < + minimum_support_points_; + ++attempt) { + bbox_minimum.array() -= expansion; + bbox_maximum.array() += expansion; + support_indices = + points_inside(points, bbox_minimum, bbox_maximum); + expansion *= 1.5; + } + + if (static_cast(support_indices.size()) >= + minimum_support_points_) { + return; + } + + std::vector> nearest; + nearest.reserve(static_cast(points.rows())); + for (Index point_i = 0; point_i < points.rows(); ++point_i) { + nearest.emplace_back( + (points.row(point_i) - sample_point).squaredNorm(), point_i); + } + std::sort(nearest.begin(), nearest.end()); + + support_indices.clear(); + auto count = std::min(minimum_support_points_, points.rows()); + for (Index i = 0; i < count; ++i) { + auto point_i = nearest.at(static_cast(i)).second; + support_indices.push_back(point_i); + bbox_minimum = bbox_minimum.cwiseMin(points.row(point_i)); + bbox_maximum = bbox_maximum.cwiseMax(points.row(point_i)); + } + } + + double overlap_; + double orientation_consistency_; + double minimum_core_size_; + double maximum_core_size_; + Index minimum_core_points_; + Index minimum_support_points_; + int maximum_depth_; +}; + +} // namespace polatory::structural diff --git a/include/polatory/structural/adaptive_domain_orientation_average.hpp b/include/polatory/structural/adaptive_domain_orientation_average.hpp new file mode 100644 index 000000000..efd7982e4 --- /dev/null +++ b/include/polatory/structural/adaptive_domain_orientation_average.hpp @@ -0,0 +1,133 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +// Experimental single-input refinement for the adaptive partitioner. +// +// The adaptive builder uses all sampled normals to decide where to split, but its +// standard build() assigns each final leaf the trend sampled at only the leaf +// centroid. This helper preserves the adaptive boxes and support memberships while +// replacing that centroid normal with a decay-weighted, sign-invariant axial +// average of the structural normals sampled at the interpolation points in the +// leaf core. +inline std::vector +build_adaptive_orientation_averaged_domains( + const geometry::Points3& points, + const std::vector& inputs, + StructuralTrendType trend_type = + StructuralTrendType::kStrongestAlongInputs, + const std::vector& model_parameters = {}, + double overlap = 0.0, + double orientation_consistency = 0.97, + double minimum_core_size = 0.0, + double maximum_core_size = 0.0, + Index minimum_core_points = 24, + Index minimum_support_points = 4, + int maximum_depth = 20) { + AdaptiveStructuralDomainBuilder3 builder( + overlap, orientation_consistency, minimum_core_size, + maximum_core_size, minimum_core_points, + minimum_support_points, maximum_depth); + + auto center_sampled_domains = + builder.build(points, inputs, trend_type, model_parameters); + + // Multiple-input orientation combination is still experimental. Keep the + // validated centre-sampled behaviour there and isolate this refinement to the + // current single-reference-mesh benchmarks. + if (inputs.size() != 1) { + return center_sampled_domains; + } + + const auto& input = inputs.front(); + auto vertex_normals = detail::averaged_vertex_normals(input); + + auto actual_overlap = overlap; + if (!(actual_overlap > 0.0)) { + actual_overlap = input.range(); + } + + std::vector averaged_domains; + averaged_domains.reserve(center_sampled_domains.size()); + + for (const auto& domain : center_sampled_domains) { + geometry::Point3 core_min = + domain.bbox().min().array() + actual_overlap; + geometry::Point3 core_max = + domain.bbox().max().array() - actual_overlap; + + std::vector core_indices; + core_indices.reserve(domain.support_indices().size()); + geometry::Point3 sample_point = geometry::Point3::Zero(); + + for (auto point_i : domain.support_indices()) { + if (detail::point_inside_box(points.row(point_i), core_min, core_max)) { + core_indices.push_back(point_i); + sample_point += points.row(point_i); + } + } + + if (!core_indices.empty()) { + sample_point /= static_cast(core_indices.size()); + } else { + sample_point = 0.5 * (core_min + core_max); + } + + auto centre_nearest = detail::nearest_trend_vertex(sample_point, input); + geometry::Point3 centre_normal = + vertex_normals.row(centre_nearest.index); + + auto centre_q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-centre_nearest.distance / input.range()); + auto ratio = 1.0 + (input.strength() - 1.0) * centre_q; + + Mat3 axial_tensor = Mat3::Zero(); + auto weight_sum = 0.0; + + for (auto point_i : core_indices) { + auto nearest = detail::nearest_trend_vertex(points.row(point_i), input); + Eigen::Vector3d normal = + vertex_normals.row(nearest.index).transpose(); + auto q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-nearest.distance / input.range()); + auto weight = std::max(q, 1e-12); + axial_tensor += weight * normal * normal.transpose(); + weight_sum += weight; + } + + geometry::Point3 representative_normal = centre_normal; + if (weight_sum > 0.0) { + axial_tensor /= weight_sum; + Eigen::SelfAdjointEigenSolver solver(axial_tensor); + if (solver.info() == Eigen::Success) { + Eigen::Vector3d axis = solver.eigenvectors().col(2); + Eigen::Vector3d reference = centre_normal.transpose(); + if (axis.dot(reference) < 0.0) { + axis *= -1.0; + } + representative_normal = axis.transpose(); + } + } + + auto anisotropy = + detail::anisotropy_from_axis(representative_normal, ratio); + + averaged_domains.emplace_back( + anisotropy, domain.bbox().min(), domain.bbox().max(), + domain.support_indices(), domain.model_parameters()); + } + + return averaged_domains; +} + +} // namespace polatory::structural diff --git a/include/polatory/structural/domain_builder.hpp b/include/polatory/structural/domain_builder.hpp new file mode 100644 index 000000000..52a1a4168 --- /dev/null +++ b/include/polatory/structural/domain_builder.hpp @@ -0,0 +1,579 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +using TriangleFaces3 = + Eigen::Matrix; + +enum class StructuralTrendType { + kStrongestAlongInputs, + kBlending, + kNonDecaying, +}; + +class StructuralTrendInput3 { + public: + StructuralTrendInput3(geometry::Points3 vertices, TriangleFaces3 faces, + double strength, double range) + : vertices_(std::move(vertices)), + faces_(std::move(faces)), + strength_(strength), + range_(range) { + if (vertices_.rows() == 0) { + throw std::invalid_argument("trend input vertices must not be empty"); + } + if (faces_.rows() == 0) { + throw std::invalid_argument("trend input faces must not be empty"); + } + if (!(strength_ >= 1.0)) { + throw std::invalid_argument("trend input strength must be at least 1"); + } + if (!(range_ > 0.0)) { + throw std::invalid_argument("trend input range must be positive"); + } + + for (Index i = 0; i < faces_.rows(); ++i) { + for (Index j = 0; j < 3; ++j) { + auto vertex_i = faces_(i, j); + if (vertex_i < 0 || vertex_i >= vertices_.rows()) { + throw std::out_of_range("trend input face index is outside vertices"); + } + } + } + } + + const TriangleFaces3& faces() const { return faces_; } + + double range() const { return range_; } + + double strength() const { return strength_; } + + const geometry::Points3& vertices() const { return vertices_; } + + private: + geometry::Points3 vertices_; + TriangleFaces3 faces_; + double strength_; + double range_; +}; + +struct StructuralTrendSamples3 { + geometry::Points3 normals; + VecX ratios; + VecX distances; + std::vector dominant_inputs; + std::vector anisotropies; +}; + +class StructuralDomainBuilder3 { + using Bbox = geometry::Bbox3; + using Point = geometry::Point3; + using Points = geometry::Points3; + + public: + explicit StructuralDomainBuilder3(double domain_size = 0.0, + double overlap = 0.0, + Index min_support_points = 4) + : domain_size_(domain_size), + overlap_(overlap), + min_support_points_(min_support_points) { + if (!(domain_size_ >= 0.0)) { + throw std::invalid_argument("domain_size must be non-negative"); + } + if (!(overlap_ >= 0.0)) { + throw std::invalid_argument("overlap must be non-negative"); + } + if (min_support_points_ <= 0) { + throw std::invalid_argument("min_support_points must be positive"); + } + } + + std::vector build( + const Points& points, const std::vector& inputs, + StructuralTrendType trend_type = + StructuralTrendType::kStrongestAlongInputs, + const std::vector& model_parameters = {}) const { + if (points.rows() == 0) { + throw std::invalid_argument("points must not be empty"); + } + validate_inputs(inputs); + + auto prepared = prepare_inputs(inputs); + + Point points_min = points.colwise().minCoeff(); + Point points_max = points.colwise().maxCoeff(); + auto widths = points_max - points_min; + + auto max_range = 0.0; + for (const auto& input : inputs) { + max_range = std::max(max_range, input.range()); + } + + auto domain_size = domain_size_; + if (!(domain_size > 0.0)) { + // The recovered Leapfrog test is best reproduced by domains whose + // active-axis core width is close to the structural range. 1.25 gives + // enough local orientation detail while the one-range overlap preserves + // continuity between neighbouring local interpolants. + domain_size = 1.25 * max_range; + if (!(domain_size > 0.0)) { + domain_size = std::max({widths(0), widths(1), widths(2)}) / 4.0; + } + if (!(domain_size > 0.0)) { + domain_size = 1.0; + } + } + + auto overlap = overlap_; + if (!(overlap > 0.0)) { + overlap = trend_type == StructuralTrendType::kNonDecaying + ? 0.5 * domain_size + : max_range; + } + + // Meshes are often extruded along one direction. Splitting that invariant + // direction creates duplicate local solves with nearly identical + // anisotropy but different support sets, which can introduce seams. Detect + // axes along which the sign-invariant normal tensor does not vary and keep + // the full point extent in those directions. + auto active_axes = detect_active_axes(inputs, prepared); + + using Cell = std::array; + std::map> cells; + + for (Index i = 0; i < points.rows(); ++i) { + Cell cell{}; + for (Index axis = 0; axis < 3; ++axis) { + if (!active_axes.at(static_cast(axis))) { + cell.at(static_cast(axis)) = 0; + continue; + } + + cell.at(static_cast(axis)) = + static_cast(std::floor( + (points(i, axis) - points_min(axis)) / domain_size)); + } + cells[cell].push_back(i); + } + + std::vector domains; + domains.reserve(cells.size()); + + for (const auto& [cell, core_indices] : cells) { + Point core_min; + Point core_max; + Point sample_point = Point::Zero(); + + for (Index axis = 0; axis < 3; ++axis) { + if (!active_axes.at(static_cast(axis))) { + core_min(axis) = points_min(axis); + core_max(axis) = points_max(axis); + continue; + } + + auto cell_i = + static_cast(cell.at(static_cast(axis))); + core_min(axis) = points_min(axis) + cell_i * domain_size; + core_max(axis) = std::min(core_min(axis) + domain_size, + points_max(axis) + domain_size * 1e-9); + } + + for (auto point_i : core_indices) { + sample_point += points.row(point_i); + } + sample_point /= static_cast(core_indices.size()); + + auto trend = evaluate_one(sample_point, inputs, prepared, trend_type); + + Point bbox_min = core_min.array() - overlap; + Point bbox_max = core_max.array() + overlap; + auto support_indices = points_inside(points, bbox_min, bbox_max); + + if (static_cast(support_indices.size()) < min_support_points_) { + expand_to_minimum_support(points, sample_point, domain_size, + support_indices, bbox_min, bbox_max); + } + + domains.emplace_back(trend.anisotropy, bbox_min, bbox_max, + std::move(support_indices), model_parameters); + } + + return domains; + } + + StructuralTrendSamples3 sample( + const Points& query_points, + const std::vector& inputs, + StructuralTrendType trend_type = + StructuralTrendType::kStrongestAlongInputs) const { + validate_inputs(inputs); + + StructuralTrendSamples3 result; + result.normals.resize(query_points.rows(), 3); + result.ratios.resize(query_points.rows()); + result.distances.resize(query_points.rows()); + result.dominant_inputs.resize( + static_cast(query_points.rows())); + result.anisotropies.reserve( + static_cast(query_points.rows())); + + auto prepared = prepare_inputs(inputs); + + for (Index i = 0; i < query_points.rows(); ++i) { + auto trend = + evaluate_one(query_points.row(i), inputs, prepared, trend_type); + result.normals.row(i) = trend.normal; + result.ratios(i) = trend.ratio; + result.distances(i) = trend.distance; + result.dominant_inputs.at(static_cast(i)) = + trend.dominant_input; + result.anisotropies.push_back(trend.anisotropy); + } + + return result; + } + + double domain_size() const { return domain_size_; } + + Index min_support_points() const { return min_support_points_; } + + double overlap() const { return overlap_; } + + private: + struct PreparedInput { + Points vertex_normals; + std::unique_ptr> tree; + }; + + struct TrendValue { + Point normal; + Mat3 anisotropy; + double ratio; + double distance; + Index dominant_input; + }; + + static Mat3 anisotropy_from_normal(const Point& normal, double ratio) { + Eigen::Vector3d n = normal.transpose(); + n.normalize(); + + auto tangent_scale = std::pow(ratio, -1.0 / 3.0); + auto normal_scale = std::pow(ratio, 2.0 / 3.0); + + Mat3 projector = n * n.transpose(); + return tangent_scale * (Mat3::Identity() - projector) + + normal_scale * projector; + } + + static Points compute_vertex_normals(const StructuralTrendInput3& input) { + Points normals = Points::Zero(input.vertices().rows(), 3); + + for (Index face_i = 0; face_i < input.faces().rows(); ++face_i) { + auto a = input.faces()(face_i, 0); + auto b = input.faces()(face_i, 1); + auto c = input.faces()(face_i, 2); + + Eigen::Vector3d va = input.vertices().row(a).transpose(); + Eigen::Vector3d vb = input.vertices().row(b).transpose(); + Eigen::Vector3d vc = input.vertices().row(c).transpose(); + + Eigen::Vector3d face_normal = (vb - va).cross(vc - va); + auto norm = face_normal.norm(); + if (!(norm > 0.0)) { + continue; + } + face_normal /= norm; + + normals.row(a) += face_normal.transpose(); + normals.row(b) += face_normal.transpose(); + normals.row(c) += face_normal.transpose(); + } + + for (Index vertex_i = 0; vertex_i < normals.rows(); ++vertex_i) { + auto norm = normals.row(vertex_i).norm(); + if (norm > 0.0) { + normals.row(vertex_i) /= norm; + } else { + normals.row(vertex_i) << 0.0, 0.0, 1.0; + } + } + + return normals; + } + + static std::vector prepare_inputs( + const std::vector& inputs) { + std::vector prepared; + prepared.reserve(inputs.size()); + + for (const auto& input : inputs) { + PreparedInput item; + item.vertex_normals = compute_vertex_normals(input); + item.tree = + std::make_unique>(input.vertices()); + prepared.push_back(std::move(item)); + } + + return prepared; + } + + static std::array detect_active_axes( + const std::vector& inputs, + const std::vector& prepared) { + constexpr Index kNumBins = 8; + constexpr double kVariationThreshold = 0.02; + + std::array active_axes{false, false, false}; + + for (Index axis = 0; axis < 3; ++axis) { + auto maximum_variation = 0.0; + + for (Index input_i = 0; + input_i < static_cast(inputs.size()); ++input_i) { + const auto& vertices = + inputs.at(static_cast(input_i)).vertices(); + const auto& normals = + prepared.at(static_cast(input_i)).vertex_normals; + + auto coordinate_min = vertices.col(axis).minCoeff(); + auto coordinate_max = vertices.col(axis).maxCoeff(); + auto coordinate_span = coordinate_max - coordinate_min; + if (!(coordinate_span > 0.0)) { + continue; + } + + Mat3 global_tensor = Mat3::Zero(); + for (Index vertex_i = 0; vertex_i < normals.rows(); ++vertex_i) { + Eigen::Vector3d normal = normals.row(vertex_i).transpose(); + global_tensor += normal * normal.transpose(); + } + global_tensor /= static_cast(normals.rows()); + + std::array bin_tensors; + std::array bin_counts{}; + for (auto& tensor : bin_tensors) { + tensor.setZero(); + } + + for (Index vertex_i = 0; vertex_i < vertices.rows(); ++vertex_i) { + auto normalized_coordinate = + (vertices(vertex_i, axis) - coordinate_min) / coordinate_span; + auto bin_i = static_cast( + std::floor(normalized_coordinate * kNumBins)); + bin_i = std::clamp(bin_i, 0, kNumBins - 1); + + Eigen::Vector3d normal = normals.row(vertex_i).transpose(); + bin_tensors.at(static_cast(bin_i)) += + normal * normal.transpose(); + ++bin_counts.at(static_cast(bin_i)); + } + + for (Index bin_i = 0; bin_i < kNumBins; ++bin_i) { + auto count = bin_counts.at(static_cast(bin_i)); + if (count == 0) { + continue; + } + + auto local_tensor = + bin_tensors.at(static_cast(bin_i)) / + static_cast(count); + maximum_variation = + std::max(maximum_variation, + (local_tensor - global_tensor).norm()); + } + } + + active_axes.at(static_cast(axis)) = + maximum_variation > kVariationThreshold; + } + + return active_axes; + } + + static void validate_inputs( + const std::vector& inputs) { + if (inputs.empty()) { + throw std::invalid_argument("structural trend inputs must not be empty"); + } + } + + TrendValue evaluate_one( + const Point& point, + const std::vector& inputs, + const std::vector& prepared, + StructuralTrendType trend_type) const { + std::vector normals(inputs.size()); + std::vector distances(inputs.size()); + std::vector q_values(inputs.size()); + std::vector contributions(inputs.size()); + + Index dominant_input = 0; + auto best_contribution = -1.0; + auto best_distance = std::numeric_limits::infinity(); + + for (Index input_i = 0; input_i < static_cast(inputs.size()); + ++input_i) { + std::vector indices; + std::vector nearest_distances; + prepared.at(static_cast(input_i)) + .tree->knn_search(point, 1, indices, nearest_distances); + + auto vertex_i = indices.at(0); + auto distance = nearest_distances.at(0); + auto normal = prepared.at(static_cast(input_i)) + .vertex_normals.row(vertex_i); + + auto q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-distance / + inputs.at(static_cast(input_i)) + .range()); + auto contribution = + (inputs.at(static_cast(input_i)).strength() - 1.0) * q; + + normals.at(static_cast(input_i)) = normal; + distances.at(static_cast(input_i)) = distance; + q_values.at(static_cast(input_i)) = q; + contributions.at(static_cast(input_i)) = contribution; + + if (contribution > best_contribution || + (contribution == best_contribution && distance < best_distance)) { + best_contribution = contribution; + best_distance = distance; + dominant_input = input_i; + } + } + + Point normal = normals.at(static_cast(dominant_input)); + auto ratio = + 1.0 + contributions.at(static_cast(dominant_input)); + + if (trend_type == StructuralTrendType::kBlending && inputs.size() > 1) { + Point blended = Point::Zero(); + auto weight_sum = 0.0; + auto q_sum = 0.0; + auto reference = normal; + + for (Index input_i = 0; input_i < static_cast(inputs.size()); + ++input_i) { + auto aligned = normals.at(static_cast(input_i)); + if (aligned.dot(reference) < 0.0) { + aligned *= -1.0; + } + + auto weight = contributions.at(static_cast(input_i)); + if (!(weight > 0.0)) { + weight = q_values.at(static_cast(input_i)); + } + + blended += weight * aligned; + weight_sum += weight; + q_sum += q_values.at(static_cast(input_i)); + } + + if (weight_sum > 0.0 && blended.norm() > 0.0) { + normal = blended / blended.norm(); + } + + auto contribution_sum = 0.0; + for (auto contribution : contributions) { + contribution_sum += contribution; + } + ratio = 1.0 + contribution_sum / std::max(1.0, q_sum); + } + + normal.normalize(); + + auto minimum_distance = + *std::min_element(distances.begin(), distances.end()); + + return TrendValue{normal, + anisotropy_from_normal(normal, ratio), + ratio, + minimum_distance, + dominant_input}; + } + + static std::vector points_inside(const Points& points, + const Point& bbox_min, + const Point& bbox_max) { + std::vector result; + result.reserve(static_cast(points.rows())); + + for (Index i = 0; i < points.rows(); ++i) { + auto inside = true; + for (Index axis = 0; axis < 3; ++axis) { + if (points(i, axis) < bbox_min(axis) || + points(i, axis) > bbox_max(axis)) { + inside = false; + break; + } + } + if (inside) { + result.push_back(i); + } + } + + return result; + } + + void expand_to_minimum_support(const Points& points, + const Point& sample_point, + double domain_size, + std::vector& support_indices, + Point& bbox_min, Point& bbox_max) const { + auto expansion = 0.5 * domain_size; + + for (int attempt = 0; + attempt < 8 && + static_cast(support_indices.size()) < min_support_points_; + ++attempt) { + bbox_min.array() -= expansion; + bbox_max.array() += expansion; + support_indices = points_inside(points, bbox_min, bbox_max); + expansion *= 1.5; + } + + if (static_cast(support_indices.size()) >= min_support_points_) { + return; + } + + std::vector> nearest; + nearest.reserve(static_cast(points.rows())); + for (Index i = 0; i < points.rows(); ++i) { + nearest.emplace_back((points.row(i) - sample_point).squaredNorm(), i); + } + std::sort(nearest.begin(), nearest.end()); + + support_indices.clear(); + auto count = std::min(min_support_points_, points.rows()); + for (Index i = 0; i < count; ++i) { + auto point_i = nearest.at(static_cast(i)).second; + support_indices.push_back(point_i); + bbox_min = bbox_min.cwiseMin(points.row(point_i)); + bbox_max = bbox_max.cwiseMax(points.row(point_i)); + } + } + + double domain_size_; + double overlap_; + Index min_support_points_; +}; + +} // namespace polatory::structural diff --git a/include/polatory/structural/domain_orientation_average.hpp b/include/polatory/structural/domain_orientation_average.hpp new file mode 100644 index 000000000..c3ba5ec5a --- /dev/null +++ b/include/polatory/structural/domain_orientation_average.hpp @@ -0,0 +1,250 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +namespace detail { + +struct NearestTrendVertex3 { + Index index; + double distance; +}; + +inline NearestTrendVertex3 nearest_trend_vertex( + const geometry::Point3& point, const StructuralTrendInput3& input) { + auto best_index = Index{0}; + auto best_squared_distance = std::numeric_limits::infinity(); + + for (Index vertex_i = 0; vertex_i < input.vertices().rows(); ++vertex_i) { + auto squared_distance = + (input.vertices().row(vertex_i) - point).squaredNorm(); + if (squared_distance < best_squared_distance) { + best_squared_distance = squared_distance; + best_index = vertex_i; + } + } + + return NearestTrendVertex3{best_index, std::sqrt(best_squared_distance)}; +} + +inline geometry::Points3 averaged_vertex_normals( + const StructuralTrendInput3& input) { + geometry::Points3 normals = + geometry::Points3::Zero(input.vertices().rows(), 3); + + for (Index face_i = 0; face_i < input.faces().rows(); ++face_i) { + auto a = input.faces()(face_i, 0); + auto b = input.faces()(face_i, 1); + auto c = input.faces()(face_i, 2); + + Eigen::Vector3d va = input.vertices().row(a).transpose(); + Eigen::Vector3d vb = input.vertices().row(b).transpose(); + Eigen::Vector3d vc = input.vertices().row(c).transpose(); + + Eigen::Vector3d face_normal = (vb - va).cross(vc - va); + auto norm = face_normal.norm(); + if (!(norm > 0.0)) { + continue; + } + face_normal /= norm; + + normals.row(a) += face_normal.transpose(); + normals.row(b) += face_normal.transpose(); + normals.row(c) += face_normal.transpose(); + } + + for (Index vertex_i = 0; vertex_i < normals.rows(); ++vertex_i) { + auto norm = normals.row(vertex_i).norm(); + if (norm > 0.0) { + normals.row(vertex_i) /= norm; + } else { + normals.row(vertex_i) << 0.0, 0.0, 1.0; + } + } + + return normals; +} + +inline bool point_inside_box(const geometry::Point3& point, + const geometry::Point3& box_min, + const geometry::Point3& box_max) { + for (Index axis = 0; axis < 3; ++axis) { + if (point(axis) < box_min(axis) || point(axis) > box_max(axis)) { + return false; + } + } + return true; +} + +inline Mat3 anisotropy_from_axis(const geometry::Point3& normal, + double ratio) { + Eigen::Vector3d n = normal.transpose(); + n.normalize(); + + auto tangent_scale = std::pow(ratio, -1.0 / 3.0); + auto normal_scale = std::pow(ratio, 2.0 / 3.0); + Mat3 projector = n * n.transpose(); + + return tangent_scale * (Mat3::Identity() - projector) + + normal_scale * projector; +} + +inline double resolved_domain_size( + double requested_domain_size, + const std::vector& inputs, + const geometry::Points3& points) { + if (requested_domain_size > 0.0) { + return requested_domain_size; + } + + auto max_range = 0.0; + for (const auto& input : inputs) { + max_range = std::max(max_range, input.range()); + } + + auto size = 1.25 * max_range; + if (!(size > 0.0)) { + auto point_min = points.colwise().minCoeff(); + auto point_max = points.colwise().maxCoeff(); + auto widths = point_max - point_min; + size = std::max({widths(0), widths(1), widths(2)}) / 4.0; + } + return size > 0.0 ? size : 1.0; +} + +inline double resolved_overlap( + double requested_overlap, double domain_size, + const std::vector& inputs, + StructuralTrendType trend_type) { + if (requested_overlap > 0.0) { + return requested_overlap; + } + + if (trend_type == StructuralTrendType::kNonDecaying) { + return 0.5 * domain_size; + } + + auto max_range = 0.0; + for (const auto& input : inputs) { + max_range = std::max(max_range, input.range()); + } + return max_range; +} + +} // namespace detail + +// Experimental single-input refinement used to test whether Leapfrog derives one +// representative domain orientation from several local structural samples rather +// than from only the domain centroid. Domain boxes and support memberships remain +// exactly those produced by StructuralDomainBuilder3; only each domain's normal is +// replaced by a sign-invariant, decay-weighted axial average of the nearest mesh +// vertex normals sampled at the interpolation points in that domain core. +inline std::vector build_orientation_averaged_domains( + const geometry::Points3& points, + const std::vector& inputs, + StructuralTrendType trend_type = + StructuralTrendType::kStrongestAlongInputs, + const std::vector& model_parameters = {}, + double domain_size = 0.0, double overlap = 0.0, + Index min_support_points = 4) { + StructuralDomainBuilder3 builder(domain_size, overlap, min_support_points); + auto center_sampled_domains = + builder.build(points, inputs, trend_type, model_parameters); + + // Multiple-input equations are still being reverse engineered. Preserve the + // validated centre-sampled behaviour there and isolate this experiment to the + // current single-reference-mesh benchmark. + if (inputs.size() != 1) { + return center_sampled_domains; + } + + const auto& input = inputs.front(); + auto vertex_normals = detail::averaged_vertex_normals(input); + auto actual_domain_size = + detail::resolved_domain_size(domain_size, inputs, points); + auto actual_overlap = detail::resolved_overlap( + overlap, actual_domain_size, inputs, trend_type); + + std::vector averaged_domains; + averaged_domains.reserve(center_sampled_domains.size()); + + for (const auto& domain : center_sampled_domains) { + geometry::Point3 core_min = + domain.bbox().min().array() + actual_overlap; + geometry::Point3 core_max = + domain.bbox().max().array() - actual_overlap; + + std::vector core_indices; + core_indices.reserve(domain.support_indices().size()); + geometry::Point3 sample_point = geometry::Point3::Zero(); + + for (auto point_i : domain.support_indices()) { + if (detail::point_inside_box(points.row(point_i), core_min, core_max)) { + core_indices.push_back(point_i); + sample_point += points.row(point_i); + } + } + + if (!core_indices.empty()) { + sample_point /= static_cast(core_indices.size()); + } else { + sample_point = 0.5 * (core_min + core_max); + } + + auto centre_nearest = detail::nearest_trend_vertex(sample_point, input); + geometry::Point3 centre_normal = + vertex_normals.row(centre_nearest.index); + + auto centre_q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-centre_nearest.distance / input.range()); + auto ratio = 1.0 + (input.strength() - 1.0) * centre_q; + + Mat3 axial_tensor = Mat3::Zero(); + auto weight_sum = 0.0; + + for (auto point_i : core_indices) { + auto nearest = detail::nearest_trend_vertex(points.row(point_i), input); + Eigen::Vector3d normal = + vertex_normals.row(nearest.index).transpose(); + auto q = trend_type == StructuralTrendType::kNonDecaying + ? 1.0 + : std::exp(-nearest.distance / input.range()); + auto weight = std::max(q, 1e-12); + axial_tensor += weight * normal * normal.transpose(); + weight_sum += weight; + } + + geometry::Point3 representative_normal = centre_normal; + if (weight_sum > 0.0) { + axial_tensor /= weight_sum; + Eigen::SelfAdjointEigenSolver solver(axial_tensor); + if (solver.info() == Eigen::Success) { + Eigen::Vector3d axis = solver.eigenvectors().col(2); + Eigen::Vector3d reference = centre_normal.transpose(); + if (axis.dot(reference) < 0.0) { + axis *= -1.0; + } + representative_normal = axis.transpose(); + } + } + + auto anisotropy = + detail::anisotropy_from_axis(representative_normal, ratio); + averaged_domains.emplace_back( + anisotropy, domain.bbox().min(), domain.bbox().max(), + domain.support_indices(), domain.model_parameters()); + } + + return averaged_domains; +} + +} // namespace polatory::structural diff --git a/include/polatory/structural/domain_spec.hpp b/include/polatory/structural/domain_spec.hpp new file mode 100644 index 000000000..e246f3e8b --- /dev/null +++ b/include/polatory/structural/domain_spec.hpp @@ -0,0 +1,48 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +class DomainSpec3 { + public: + DomainSpec3(const Mat3& anisotropy, const geometry::Point3& bbox_min, + const geometry::Point3& bbox_max, std::vector support_indices, + std::vector model_parameters = {}) + : anisotropy_(anisotropy), + bbox_(bbox_min, bbox_max), + support_indices_(std::move(support_indices)), + model_parameters_(std::move(model_parameters)) { + if (!(anisotropy_.determinant() > 0.0)) { + throw std::invalid_argument("anisotropy must have a positive determinant"); + } + if (bbox_.is_empty()) { + throw std::invalid_argument("domain bbox must not be empty"); + } + if (support_indices_.empty()) { + throw std::invalid_argument("support_indices must not be empty"); + } + } + + const Mat3& anisotropy() const { return anisotropy_; } + + const geometry::Bbox3& bbox() const { return bbox_; } + + const std::vector& model_parameters() const { return model_parameters_; } + + const std::vector& support_indices() const { return support_indices_; } + + private: + Mat3 anisotropy_; + geometry::Bbox3 bbox_; + std::vector support_indices_; + std::vector model_parameters_; +}; + +} // namespace polatory::structural diff --git a/include/polatory/structural/interpolant.hpp b/include/polatory/structural/interpolant.hpp new file mode 100644 index 000000000..7323861e8 --- /dev/null +++ b/include/polatory/structural/interpolant.hpp @@ -0,0 +1,399 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace polatory::structural { + +class StructuralInterpolant3 { + static constexpr double kInfinity = std::numeric_limits::infinity(); + using Bbox = geometry::Bbox3; + using Model = Model<3>; + using Point = geometry::Point3; + using Points = geometry::Points3; + using StandardInterpolant = Interpolant<3>; + + public: + explicit StructuralInterpolant3(const Model& base_model, + double outside_value = -1.0, + double blend_power = 1.0, + double alignment_strength = 0.0) + : base_model_(base_model), + outside_value_(outside_value), + blend_power_(blend_power), + alignment_strength_(alignment_strength) { + if (!(blend_power_ > 0.0)) { + throw std::invalid_argument("blend_power must be positive"); + } + if (!(alignment_strength_ >= 0.0)) { + throw std::invalid_argument("alignment_strength must be non-negative"); + } + } + + const Bbox& bbox() const { + throw_if_not_fitted(); + return bbox_; + } + + double blend_power() const { return blend_power_; } + + double alignment_strength() const { return alignment_strength_; } + + Index num_domains() const { return static_cast(domains_.size()); } + + double outside_value() const { return outside_value_; } + + std::vector domain_offsets() const { + throw_if_not_fitted(); + std::vector offsets; + offsets.reserve(domains_.size()); + for (const auto& domain : domains_) { + offsets.push_back(domain.offset); + } + return offsets; + } + + void fit(const Points& points, const VecX& values, + const std::vector& domain_specs, + double tolerance, int max_iter = 100, + double accuracy = kInfinity) { + if (points.rows() != values.rows()) { + throw std::invalid_argument("values.rows() must equal points.rows()"); + } + if (points.rows() == 0) { + throw std::invalid_argument("points must not be empty"); + } + if (domain_specs.empty()) { + throw std::invalid_argument("domain_specs must not be empty"); + } + if (!(tolerance > 0.0)) { + throw std::invalid_argument("tolerance must be positive"); + } + if (max_iter < 0) { + throw std::invalid_argument("max_iter must be nonnegative"); + } + if (!(accuracy > 0.0)) { + throw std::invalid_argument("accuracy must be positive"); + } + + clear(); + + for (const auto& spec : domain_specs) { + const auto& indices = spec.support_indices(); + Points local_points(static_cast(indices.size()), 3); + VecX local_values(static_cast(indices.size())); + + for (Index local_i = 0; + local_i < static_cast(indices.size()); ++local_i) { + auto global_i = indices.at(static_cast(local_i)); + if (global_i < 0 || global_i >= points.rows()) { + throw std::out_of_range("domain support index is outside points"); + } + local_points.row(local_i) = points.row(global_i); + local_values(local_i) = values(global_i); + } + + Model local_model = base_model_; + if (!spec.model_parameters().empty()) { + local_model.set_parameters(spec.model_parameters()); + } + for (auto& rbf : local_model.rbfs()) { + rbf.set_anisotropy(spec.anisotropy()); + } + + auto local_interpolant = + std::make_unique(local_model); + local_interpolant->fit(local_points, local_values, tolerance, + max_iter, accuracy); + + domains_.push_back( + Domain{spec, std::move(local_interpolant), 0.0}); + bbox_ = bbox_.is_empty() ? spec.bbox() + : bbox_.convex_hull(spec.bbox()); + } + + if (alignment_strength_ > 0.0 && domains_.size() > 1) { + align_domain_offsets(accuracy); + } + + fitted_ = true; + } + + VecX evaluate(const Points& points, double accuracy = kInfinity) { + throw_if_not_fitted(); + if (!(accuracy > 0.0)) { + throw std::invalid_argument("accuracy must be positive"); + } + if (points.rows() == 0) { + return VecX(); + } + + set_evaluation_bbox_impl(Bbox::from_points(points), accuracy); + return evaluate_impl(points); + } + + VecX evaluate_impl(const Points& points) const { + throw_if_not_fitted(); + + VecX numerator = VecX::Zero(points.rows()); + VecX denominator = VecX::Zero(points.rows()); + + for (const auto& domain : domains_) { + std::vector active_indices; + std::vector active_weights; + active_indices.reserve(static_cast(points.rows())); + active_weights.reserve(static_cast(points.rows())); + + for (Index i = 0; i < points.rows(); ++i) { + auto weight = box_weight(points.row(i), domain.spec.bbox()); + if (weight > 0.0) { + active_indices.push_back(i); + active_weights.push_back(weight); + } + } + + if (active_indices.empty()) { + continue; + } + + Points active_points(static_cast(active_indices.size()), 3); + for (Index i = 0; + i < static_cast(active_indices.size()); ++i) { + active_points.row(i) = + points.row(active_indices.at(static_cast(i))); + } + + VecX predictions = domain.interpolant->evaluate_impl(active_points); + predictions.array() += domain.offset; + for (Index i = 0; + i < static_cast(active_indices.size()); ++i) { + auto query_i = active_indices.at(static_cast(i)); + auto weight = active_weights.at(static_cast(i)); + numerator(query_i) += weight * predictions(i); + denominator(query_i) += weight; + } + } + + VecX result = VecX::Constant(points.rows(), outside_value_); + for (Index i = 0; i < points.rows(); ++i) { + if (denominator(i) > 0.0) { + result(i) = numerator(i) / denominator(i); + } + } + return result; + } + + void set_evaluation_bbox_impl(const Bbox& bbox, + double accuracy = kInfinity) { + throw_if_not_fitted(); + if (!(accuracy > 0.0)) { + throw std::invalid_argument("accuracy must be positive"); + } + + for (auto& domain : domains_) { + domain.interpolant->set_evaluation_bbox_impl(bbox, accuracy); + } + } + + private: + struct Domain { + DomainSpec3 spec; + std::unique_ptr interpolant; + double offset; + }; + + double box_weight(const Point& point, const Bbox& bbox) const { + if (!bbox.contains(point)) { + return 0.0; + } + + double weight = 1.0; + auto width = bbox.width(); + for (Index axis = 0; axis < 3; ++axis) { + if (!(width(axis) > 0.0)) { + return 0.0; + } + + auto distance_to_face = + std::min(point(axis) - bbox.min()(axis), + bbox.max()(axis) - point(axis)); + auto u = std::clamp(2.0 * distance_to_face / width(axis), + 0.0, 1.0); + auto smooth = u * u * (3.0 - 2.0 * u); + weight *= smooth; + } + + return std::pow(weight, blend_power_); + } + + static bool overlap_bbox(const Bbox& a, const Bbox& b, + Point& overlap_min, Point& overlap_max) { + overlap_min = a.min().cwiseMax(b.min()); + overlap_max = a.max().cwiseMin(b.max()); + return (overlap_max.array() > overlap_min.array()).all(); + } + + static Points make_overlap_samples(const Point& overlap_min, + const Point& overlap_max) { + constexpr std::array fractions{ + 0.0714285714285714, + 0.2142857142857143, + 0.3571428571428571, + 0.5, + 0.6428571428571429, + 0.7857142857142857, + 0.9285714285714286}; + Points samples(343, 3); + Index sample_i = 0; + for (auto fx : fractions) { + for (auto fy : fractions) { + for (auto fz : fractions) { + samples.row(sample_i) = + overlap_min.array() + + Point(fx, fy, fz).array() * + (overlap_max - overlap_min).array(); + ++sample_i; + } + } + } + return samples; + } + + void align_domain_offsets(double accuracy) { + auto n = static_cast(domains_.size()); + Eigen::MatrixXd system = Eigen::MatrixXd::Zero(n, n); + Eigen::VectorXd rhs = Eigen::VectorXd::Zero(n); + Index edge_count = 0; + + for (Index i = 0; i < n; ++i) { + for (Index j = i + 1; j < n; ++j) { + Point overlap_min; + Point overlap_max; + if (!overlap_bbox(domains_.at(static_cast(i)).spec.bbox(), + domains_.at(static_cast(j)).spec.bbox(), + overlap_min, overlap_max)) { + continue; + } + + auto samples = make_overlap_samples(overlap_min, overlap_max); + auto predictions_i = + domains_.at(static_cast(i)) + .interpolant->evaluate(samples, accuracy); + auto predictions_j = + domains_.at(static_cast(j)) + .interpolant->evaluate(samples, accuracy); + + std::vector level_distances( + static_cast(samples.rows())); + for (Index sample_i = 0; sample_i < samples.rows(); ++sample_i) { + level_distances.at(static_cast(sample_i)) = + std::abs(predictions_i(sample_i)) + + std::abs(predictions_j(sample_i)); + } + auto middle = level_distances.begin() + + static_cast(level_distances.size() / 2); + std::nth_element(level_distances.begin(), middle, + level_distances.end()); + auto level_scale = std::max(*middle, 1e-8); + + auto weighted_difference = 0.0; + auto weight_sum = 0.0; + for (Index sample_i = 0; sample_i < samples.rows(); ++sample_i) { + auto wi = box_weight( + samples.row(sample_i), + domains_.at(static_cast(i)).spec.bbox()); + auto wj = box_weight( + samples.row(sample_i), + domains_.at(static_cast(j)).spec.bbox()); + auto overlap_weight = std::sqrt(wi * wj); + + auto level_distance = + std::abs(predictions_i(sample_i)) + + std::abs(predictions_j(sample_i)); + auto normalized_level = level_distance / level_scale; + auto level_weight = + 1.0 / (1.0 + normalized_level * normalized_level); + auto weight = overlap_weight * level_weight; + weighted_difference += + weight * (predictions_j(sample_i) - + predictions_i(sample_i)); + weight_sum += weight; + } + + if (!(weight_sum > 1e-12)) { + continue; + } + + auto difference = weighted_difference / weight_sum; + auto confidence = weight_sum / + static_cast(samples.rows()); + + system(i, i) += confidence; + system(j, j) += confidence; + system(i, j) -= confidence; + system(j, i) -= confidence; + rhs(i) += confidence * difference; + rhs(j) -= confidence * difference; + ++edge_count; + } + } + + if (edge_count == 0) { + return; + } + + // The pairwise system is invariant to a common global offset. A tiny ridge + // gives disconnected overlap components a stable minimum-norm solution. + system.diagonal().array() += 1e-10; + Eigen::VectorXd offsets = system.ldlt().solve(rhs); + if (!offsets.allFinite()) { + return; + } + + // Preserve the model's global zero reference rather than translating every + // local field together. + offsets.array() -= offsets.mean(); + + for (Index i = 0; i < n; ++i) { + domains_.at(static_cast(i)).offset = + alignment_strength_ * offsets(i); + } + } + + void clear() { + fitted_ = false; + domains_.clear(); + bbox_ = Bbox(); + } + + void throw_if_not_fitted() const { + if (!fitted_) { + throw std::runtime_error( + "structural interpolant has not been fitted"); + } + } + + Model base_model_; + double outside_value_; + double blend_power_; + double alignment_strength_; + bool fitted_{}; + std::vector domains_; + Bbox bbox_; +}; + +} // namespace polatory::structural diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index ac57f7222..2b14d0e67 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -2,14 +2,23 @@ set(TARGET _core) pybind11_add_module(${TARGET} python_binding.cpp) -target_compile_definitions(${TARGET} PRIVATE - -DPOLATORY_VERSION=${POLATORY_VERSION} -) +set(STRUCTURAL_TARGET _structural) +pybind11_add_module(${STRUCTURAL_TARGET} structural_binding.cpp) -target_link_libraries(${TARGET} PRIVATE - polatory -) +foreach(PY_TARGET ${TARGET} ${STRUCTURAL_TARGET}) + target_compile_definitions(${PY_TARGET} PRIVATE + -DPOLATORY_VERSION=${POLATORY_VERSION} + ) -if(MSVC) - polatory_target_contents(${TARGET} ${POLATORY_DLLS}) -endif() + target_link_libraries(${PY_TARGET} PRIVATE + polatory + ) + + if(MSVC) + polatory_target_contents(${PY_TARGET} ${POLATORY_DLLS}) + endif() +endforeach() + +# setup.py builds the _core target explicitly. Make that build the companion +# structural extension into the same Python package output directory as well. +add_dependencies(${TARGET} ${STRUCTURAL_TARGET}) diff --git a/python/src/polatory/__init__.py b/python/src/polatory/__init__.py index 72aa94182..62f001fea 100644 --- a/python/src/polatory/__init__.py +++ b/python/src/polatory/__init__.py @@ -1,2 +1,7 @@ from ._core import * from ._core import __doc__, __version__ +from ._structural import * +from .clustered_domain_builder import ( + ClusteredStructuralDomainBuilder3, + fit_from_meshes_clustered, +) diff --git a/python/src/polatory/clustered_domain_builder.py b/python/src/polatory/clustered_domain_builder.py new file mode 100644 index 000000000..89135c887 --- /dev/null +++ b/python/src/polatory/clustered_domain_builder.py @@ -0,0 +1,440 @@ +"""Deterministic clustered structural-domain builder for local-varying anisotropy. + +The builder samples the structural anisotropy at every interpolation point, +clusters points jointly in active spatial coordinates and standardized +anisotropy-matrix space, selects an actual sampled matrix medoid for each +cluster, and expands each core by the anisotropic local RBF range. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from math import ceil, sqrt +from typing import Sequence + +import numpy as np + +from ._structural import ( + StructuralDomain3, + StructuralDomainBuilder3, + StructuralTrendType, +) + + +@dataclass(frozen=True) +class _KMeansResult: + labels: np.ndarray + centers: np.ndarray + inertia: float + + +class ClusteredStructuralDomainBuilder3: + """Build overlapping structural domains by deterministic matrix-space clustering. + + Parameters + ---------- + domain_count: + Number of domains. ``0`` selects ``ceil(sqrt(n_points) / 3)``. + base_range: + Base RBF range before applying ``local_range_scale``. ``0`` derives it + from the final model parameter when available, otherwise from the + largest structural-trend range. + local_range_scale: + Scale applied to the base RBF range for every local solve. The tested + default is ``0.8``. + restarts: + Number of deterministic k-means++ restarts. + random_seed: + Seed controlling the deterministic restarts. + maximum_iterations: + Maximum Lloyd iterations per restart. + minimum_support_points: + Minimum number of interpolation points in an expanded domain box. + spatial_weight, anisotropy_weight: + Optional feature-group weights. Defaults reproduce the validated + implementation. + """ + + def __init__( + self, + domain_count: int = 0, + base_range: float = 0.0, + local_range_scale: float = 0.8, + restarts: int = 20, + random_seed: int = 42, + maximum_iterations: int = 200, + minimum_support_points: int = 4, + spatial_weight: float = 1.0, + anisotropy_weight: float = 1.0, + ) -> None: + if domain_count < 0: + raise ValueError("domain_count must be non-negative") + if base_range < 0.0: + raise ValueError("base_range must be non-negative") + if not local_range_scale > 0.0: + raise ValueError("local_range_scale must be positive") + if restarts <= 0: + raise ValueError("restarts must be positive") + if maximum_iterations <= 0: + raise ValueError("maximum_iterations must be positive") + if minimum_support_points <= 0: + raise ValueError("minimum_support_points must be positive") + if not spatial_weight > 0.0: + raise ValueError("spatial_weight must be positive") + if not anisotropy_weight > 0.0: + raise ValueError("anisotropy_weight must be positive") + + self.domain_count = int(domain_count) + self.base_range = float(base_range) + self.local_range_scale = float(local_range_scale) + self.restarts = int(restarts) + self.random_seed = int(random_seed) + self.maximum_iterations = int(maximum_iterations) + self.minimum_support_points = int(minimum_support_points) + self.spatial_weight = float(spatial_weight) + self.anisotropy_weight = float(anisotropy_weight) + + self.labels_: np.ndarray | None = None + self.active_axes_: np.ndarray | None = None + self.medoid_indices_: np.ndarray | None = None + self.features_: np.ndarray | None = None + self.domain_count_: int | None = None + self.local_range_: float | None = None + self.inertia_: float | None = None + + @staticmethod + def _matrix_features(anisotropies: np.ndarray) -> np.ndarray: + return np.column_stack( + [ + anisotropies[:, 0, 0], + anisotropies[:, 0, 1], + anisotropies[:, 0, 2], + anisotropies[:, 1, 1], + anisotropies[:, 1, 2], + anisotropies[:, 2, 2], + ] + ) + + @classmethod + def _detect_active_axes( + cls, + points: np.ndarray, + anisotropies: np.ndarray, + number_of_bins: int = 8, + variation_threshold: float = 0.02, + ) -> np.ndarray: + matrix_features = cls._matrix_features(anisotropies) + global_mean = matrix_features.mean(axis=0) + active = np.zeros(3, dtype=bool) + + for axis in range(3): + coordinate_min = float(points[:, axis].min()) + coordinate_max = float(points[:, axis].max()) + span = coordinate_max - coordinate_min + if not span > 0.0: + continue + + normalized = (points[:, axis] - coordinate_min) / span + bins = np.clip( + np.floor(normalized * number_of_bins).astype(np.int64), + 0, + number_of_bins - 1, + ) + + maximum_variation = 0.0 + for bin_id in range(number_of_bins): + indices = bins == bin_id + if not np.any(indices): + continue + local_mean = matrix_features[indices].mean(axis=0) + maximum_variation = max( + maximum_variation, + float(np.linalg.norm(local_mean - global_mean)), + ) + + active[axis] = maximum_variation > variation_threshold + + if not np.any(active): + active[:] = True + return active + + @staticmethod + def _kmeans( + features: np.ndarray, + cluster_count: int, + restarts: int, + random_seed: int, + maximum_iterations: int, + ) -> _KMeansResult: + random = np.random.default_rng(random_seed) + point_count = len(features) + best: _KMeansResult | None = None + + for _ in range(restarts): + centers = np.empty((cluster_count, features.shape[1]), dtype=float) + first_index = int(random.integers(point_count)) + centers[0] = features[first_index] + + minimum_distance_squared = np.sum( + (features - centers[0]) ** 2, + axis=1, + ) + + for center_id in range(1, cluster_count): + total = float(minimum_distance_squared.sum()) + if total <= 0.0: + point_index = int(random.integers(point_count)) + else: + probabilities = minimum_distance_squared / total + point_index = int(random.choice(point_count, p=probabilities)) + + centers[center_id] = features[point_index] + new_distance_squared = np.sum( + (features - centers[center_id]) ** 2, + axis=1, + ) + minimum_distance_squared = np.minimum( + minimum_distance_squared, + new_distance_squared, + ) + + labels: np.ndarray | None = None + + for _ in range(maximum_iterations): + distance_squared = np.sum( + ( + features[:, None, :] + - centers[None, :, :] + ) + ** 2, + axis=2, + ) + new_labels = np.argmin(distance_squared, axis=1) + + if labels is not None and np.array_equal(new_labels, labels): + break + labels = new_labels + + for cluster_id in range(cluster_count): + indices = np.flatnonzero(labels == cluster_id) + if len(indices) > 0: + centers[cluster_id] = features[indices].mean(axis=0) + else: + centers[cluster_id] = features[ + int(random.integers(point_count)) + ] + + assert labels is not None + inertia = float(np.sum((features - centers[labels]) ** 2)) + candidate = _KMeansResult(labels.copy(), centers.copy(), inertia) + if best is None or candidate.inertia < best.inertia: + best = candidate + + assert best is not None + return best + + @staticmethod + def _matrix_medoid( + anisotropies: np.ndarray, + indices: np.ndarray, + ) -> tuple[int, np.ndarray]: + local = anisotropies[indices] + average = local.mean(axis=0) + distances = np.linalg.norm(local - average, axis=(1, 2)) + local_medoid = int(np.argmin(distances)) + global_medoid = int(indices[local_medoid]) + return global_medoid, anisotropies[global_medoid] + + @staticmethod + def _support_indices( + points: np.ndarray, + minimum: np.ndarray, + maximum: np.ndarray, + ) -> np.ndarray: + return np.flatnonzero( + np.all( + (points >= minimum) & (points <= maximum), + axis=1, + ) + ).astype(np.int64) + + def _resolved_domain_count(self, point_count: int) -> int: + if self.domain_count > 0: + return min(self.domain_count, point_count) + return min(point_count, max(2, int(ceil(sqrt(point_count) / 3.0)))) + + def _resolved_base_range( + self, + inputs: Sequence[object], + model_parameters: Sequence[float], + ) -> float: + if self.base_range > 0.0: + return self.base_range + if len(model_parameters) >= 3 and float(model_parameters[-1]) > 0.0: + return float(model_parameters[-1]) + + maximum_input_range = max(float(item.range) for item in inputs) + if not maximum_input_range > 0.0: + raise ValueError("could not resolve a positive base range") + return maximum_input_range + + def build( + self, + points: np.ndarray, + inputs: Sequence[object], + trend_type: object = StructuralTrendType.STRONGEST_ALONG_INPUTS, + model_parameters: Sequence[float] = (), + ) -> list[StructuralDomain3]: + points = np.asarray(points, dtype=float) + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError("points must have shape (n, 3)") + if len(points) == 0: + raise ValueError("points must not be empty") + if not inputs: + raise ValueError("inputs must not be empty") + + samples = StructuralDomainBuilder3().sample( + points, + list(inputs), + trend_type, + ) + anisotropies = np.asarray(samples.anisotropies, dtype=float) + if anisotropies.shape != (len(points), 3, 3): + raise RuntimeError("unexpected sampled anisotropy shape") + + active_axes = self._detect_active_axes(points, anisotropies) + maximum_input_range = max(float(item.range) for item in inputs) + + spatial = ( + points[:, active_axes] + - points[:, active_axes].mean(axis=0) + ) / maximum_input_range + + matrix_features = self._matrix_features(anisotropies) + matrix_mean = matrix_features.mean(axis=0) + matrix_scale = matrix_features.std(axis=0) + matrix_scale[matrix_scale < 1e-12] = 1.0 + matrix_features = (matrix_features - matrix_mean) / matrix_scale + + features = np.column_stack( + [ + self.spatial_weight * spatial, + self.anisotropy_weight * matrix_features, + ] + ) + + cluster_count = self._resolved_domain_count(len(points)) + result = self._kmeans( + features, + cluster_count, + self.restarts, + self.random_seed, + self.maximum_iterations, + ) + + base_range = self._resolved_base_range(inputs, model_parameters) + local_range = self.local_range_scale * base_range + if not local_range > 0.0: + raise ValueError("resolved local range must be positive") + + domains: list[StructuralDomain3] = [] + medoid_indices: list[int] = [] + + for cluster_id in range(cluster_count): + core_indices = np.flatnonzero(result.labels == cluster_id) + if len(core_indices) == 0: + continue + + medoid_index, anisotropy = self._matrix_medoid( + anisotropies, + core_indices, + ) + medoid_indices.append(medoid_index) + + core_minimum = points[core_indices].min(axis=0) + core_maximum = points[core_indices].max(axis=0) + + inverse = np.linalg.inv(anisotropy) + expansion = local_range * np.linalg.norm(inverse, axis=1) + bbox_minimum = core_minimum - expansion + bbox_maximum = core_maximum + expansion + + support_indices = self._support_indices( + points, + bbox_minimum, + bbox_maximum, + ) + + attempts = 0 + while ( + len(support_indices) < self.minimum_support_points + and attempts < 8 + ): + expansion *= 1.25 + bbox_minimum = core_minimum - expansion + bbox_maximum = core_maximum + expansion + support_indices = self._support_indices( + points, + bbox_minimum, + bbox_maximum, + ) + attempts += 1 + + local_parameters = list(float(value) for value in model_parameters) + if local_parameters: + local_parameters[-1] = local_range + + domains.append( + StructuralDomain3( + anisotropy=anisotropy, + bbox_min=bbox_minimum, + bbox_max=bbox_maximum, + support_indices=support_indices.tolist(), + model_parameters=local_parameters, + ) + ) + + self.labels_ = result.labels.copy() + self.active_axes_ = active_axes.copy() + self.medoid_indices_ = np.asarray(medoid_indices, dtype=np.int64) + self.features_ = features.copy() + self.domain_count_ = len(domains) + self.local_range_ = float(local_range) + self.inertia_ = float(result.inertia) + + return domains + + +def fit_from_meshes_clustered( + interpolant: object, + points: np.ndarray, + values: np.ndarray, + inputs: Sequence[object], + tolerance: float, + *, + trend_type: object = StructuralTrendType.STRONGEST_ALONG_INPUTS, + model_parameters: Sequence[float] = (), + max_iter: int = 100, + accuracy: float = float("inf"), + builder: ClusteredStructuralDomainBuilder3 | None = None, +) -> list[StructuralDomain3]: + """Build clustered domains, fit a structural interpolant, and return domains.""" + + if builder is None: + builder = ClusteredStructuralDomainBuilder3() + + domains = builder.build( + points, + inputs, + trend_type, + model_parameters, + ) + interpolant.fit( + points, + values, + domains, + tolerance, + max_iter, + accuracy, + ) + return domains diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp new file mode 100644 index 000000000..6ac54e3ff --- /dev/null +++ b/python/structural_binding.cpp @@ -0,0 +1,269 @@ +#include +#undef _GNU_SOURCE +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; +using namespace py::literals; +using namespace polatory; + +static constexpr double kInfinity = std::numeric_limits::infinity(); + +PYBIND11_MODULE(_structural, m) { + // Register the standard Polatory types (Model<3>, FieldFunction, etc.) first. + py::module_::import("polatory._core"); + + using AdaptiveDomainBuilder = structural::AdaptiveStructuralDomainBuilder3; + using DomainBuilder = structural::StructuralDomainBuilder3; + using DomainSpec = structural::DomainSpec3; + using StructuralInterpolant = structural::StructuralInterpolant3; + using TrendInput = structural::StructuralTrendInput3; + using TrendSamples = structural::StructuralTrendSamples3; + using TrendType = structural::StructuralTrendType; + + py::enum_(m, "StructuralTrendType") + .value("STRONGEST_ALONG_INPUTS", TrendType::kStrongestAlongInputs) + .value("BLENDING", TrendType::kBlending) + .value("NON_DECAYING", TrendType::kNonDecaying) + .export_values(); + + py::class_(m, "StructuralTrendInput3") + .def(py::init(), + "vertices"_a, "faces"_a, "strength"_a, "range"_a) + .def_property_readonly("vertices", &TrendInput::vertices) + .def_property_readonly("faces", &TrendInput::faces) + .def_property_readonly("strength", &TrendInput::strength) + .def_property_readonly("range", &TrendInput::range); + + py::class_(m, "StructuralTrendSamples3") + .def_property_readonly("normals", [](const TrendSamples& samples) { + return samples.normals; + }) + .def_property_readonly("ratios", [](const TrendSamples& samples) { + return samples.ratios; + }) + .def_property_readonly("distances", [](const TrendSamples& samples) { + return samples.distances; + }) + .def_property_readonly("dominant_inputs", [](const TrendSamples& samples) { + return samples.dominant_inputs; + }) + .def_property_readonly("anisotropies", [](const TrendSamples& samples) { + return samples.anisotropies; + }); + + py::class_(m, "StructuralDomain3") + .def(py::init([](const Mat3& anisotropy, const Eigen::Vector3d& bbox_min, + const Eigen::Vector3d& bbox_max, + std::vector support_indices, + std::vector model_parameters) { + geometry::Point3 min_row = bbox_min.transpose(); + geometry::Point3 max_row = bbox_max.transpose(); + return DomainSpec(anisotropy, min_row, max_row, + std::move(support_indices), + std::move(model_parameters)); + }), + "anisotropy"_a, "bbox_min"_a, "bbox_max"_a, + "support_indices"_a, + "model_parameters"_a = std::vector{}) + .def_property_readonly("anisotropy", &DomainSpec::anisotropy) + .def_property_readonly("bbox_min", [](const DomainSpec& spec) { + return Eigen::Vector3d(spec.bbox().min().transpose()); + }) + .def_property_readonly("bbox_max", [](const DomainSpec& spec) { + return Eigen::Vector3d(spec.bbox().max().transpose()); + }) + .def_property_readonly("support_indices", &DomainSpec::support_indices) + .def_property_readonly("model_parameters", &DomainSpec::model_parameters); + + py::class_(m, "StructuralDomainBuilder3") + .def(py::init(), + "domain_size"_a = 0.0, + "overlap"_a = 0.0, + "min_support_points"_a = 4) + .def_property_readonly("domain_size", &DomainBuilder::domain_size) + .def_property_readonly("overlap", &DomainBuilder::overlap) + .def_property_readonly("min_support_points", + &DomainBuilder::min_support_points) + .def("build", &DomainBuilder::build, + "points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "model_parameters"_a = std::vector{}) + .def( + "build_orientation_averaged", + [](const DomainBuilder& builder, + const geometry::Points3& points, + const std::vector& inputs, + TrendType trend_type, + const std::vector& model_parameters) { + return structural::build_orientation_averaged_domains( + points, inputs, trend_type, model_parameters, + builder.domain_size(), builder.overlap(), + builder.min_support_points()); + }, + "points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "model_parameters"_a = std::vector{}) + .def("sample", &DomainBuilder::sample, + "query_points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs); + + py::class_(m, "AdaptiveStructuralDomainBuilder3") + .def(py::init(), + "overlap"_a = 0.0, + "orientation_consistency"_a = 0.97, + "minimum_core_size"_a = 0.0, + "maximum_core_size"_a = 0.0, + "minimum_core_points"_a = 24, + "minimum_support_points"_a = 4, + "maximum_depth"_a = 20) + .def_property_readonly("overlap", + &AdaptiveDomainBuilder::overlap) + .def_property_readonly( + "orientation_consistency", + &AdaptiveDomainBuilder::orientation_consistency) + .def_property_readonly("minimum_core_size", + &AdaptiveDomainBuilder::minimum_core_size) + .def_property_readonly("maximum_core_size", + &AdaptiveDomainBuilder::maximum_core_size) + .def_property_readonly("minimum_core_points", + &AdaptiveDomainBuilder::minimum_core_points) + .def_property_readonly("minimum_support_points", + &AdaptiveDomainBuilder::minimum_support_points) + .def_property_readonly("maximum_depth", + &AdaptiveDomainBuilder::maximum_depth) + .def("build", &AdaptiveDomainBuilder::build, + "points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "model_parameters"_a = std::vector{}) + .def( + "build_orientation_averaged", + [](const AdaptiveDomainBuilder& builder, + const geometry::Points3& points, + const std::vector& inputs, + TrendType trend_type, + const std::vector& model_parameters) { + return structural::build_adaptive_orientation_averaged_domains( + points, inputs, trend_type, model_parameters, + builder.overlap(), builder.orientation_consistency(), + builder.minimum_core_size(), builder.maximum_core_size(), + builder.minimum_core_points(), builder.minimum_support_points(), + builder.maximum_depth()); + }, + "points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "model_parameters"_a = std::vector{}); + + py::class_(m, "StructuralInterpolant3") + .def(py::init&, double, double, double>(), + "base_model"_a, + "outside_value"_a = -1.0, + "blend_power"_a = 1.0, + "alignment_strength"_a = 0.0) + .def_property_readonly("bbox_min", [](const StructuralInterpolant& interpolant) { + return Eigen::Vector3d(interpolant.bbox().min().transpose()); + }) + .def_property_readonly("bbox_max", [](const StructuralInterpolant& interpolant) { + return Eigen::Vector3d(interpolant.bbox().max().transpose()); + }) + .def_property_readonly("blend_power", &StructuralInterpolant::blend_power) + .def_property_readonly("alignment_strength", + &StructuralInterpolant::alignment_strength) + .def_property_readonly("domain_offsets", + &StructuralInterpolant::domain_offsets) + .def_property_readonly("num_domains", &StructuralInterpolant::num_domains) + .def_property_readonly("outside_value", &StructuralInterpolant::outside_value) + .def("fit", &StructuralInterpolant::fit, "points"_a, "values"_a, + "domains"_a, "tolerance"_a, "max_iter"_a = 100, + "accuracy"_a = kInfinity) + .def( + "fit_from_meshes", + [](StructuralInterpolant& interpolant, + const geometry::Points3& points, + const VecX& values, + const std::vector& inputs, + double tolerance, + TrendType trend_type, + int max_iter, + double accuracy, + double domain_size, + double overlap, + Index min_support_points) { + DomainBuilder builder(domain_size, overlap, min_support_points); + auto domains = builder.build(points, inputs, trend_type); + interpolant.fit(points, values, domains, tolerance, max_iter, + accuracy); + return domains; + }, + "points"_a, + "values"_a, + "inputs"_a, + "tolerance"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "max_iter"_a = 100, + "accuracy"_a = kInfinity, + "domain_size"_a = 0.0, + "overlap"_a = 0.0, + "min_support_points"_a = 4) + .def( + "fit_from_meshes_adaptive", + [](StructuralInterpolant& interpolant, + const geometry::Points3& points, + const VecX& values, + const std::vector& inputs, + double tolerance, + TrendType trend_type, + int max_iter, + double accuracy, + double overlap, + double orientation_consistency, + double minimum_core_size, + double maximum_core_size, + Index minimum_core_points, + Index minimum_support_points, + int maximum_depth) { + AdaptiveDomainBuilder builder( + overlap, orientation_consistency, minimum_core_size, + maximum_core_size, minimum_core_points, + minimum_support_points, maximum_depth); + auto domains = builder.build(points, inputs, trend_type); + interpolant.fit(points, values, domains, tolerance, max_iter, + accuracy); + return domains; + }, + "points"_a, + "values"_a, + "inputs"_a, + "tolerance"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "max_iter"_a = 100, + "accuracy"_a = kInfinity, + "overlap"_a = 0.0, + "orientation_consistency"_a = 0.97, + "minimum_core_size"_a = 0.0, + "maximum_core_size"_a = 0.0, + "minimum_core_points"_a = 24, + "minimum_support_points"_a = 4, + "maximum_depth"_a = 20) + .def("evaluate", &StructuralInterpolant::evaluate, "points"_a, + "accuracy"_a = kInfinity); + + py::class_( + m, "StructuralRbfFieldFunction") + .def(py::init(), "interpolant"_a, + "accuracy"_a = kInfinity); +} diff --git a/setup.py b/setup.py index 6a68b8ac6..3d1e521a8 100644 --- a/setup.py +++ b/setup.py @@ -41,7 +41,7 @@ def windows_buildenv() -> dict[str, str]: vcvars64 = (vs_dir / "VC/Auxiliary/Build/vcvars64.bat").resolve() output = subprocess.run( - f'"{vcvars64}" -vcvars_ver=14.39 && set', stdout=subprocess.PIPE, check=True, text=True + f'"{vcvars64}" && set', stdout=subprocess.PIPE, check=True, text=True ).stdout env = {} diff --git a/src/point_cloud/sdf_data_generator.cpp b/src/point_cloud/sdf_data_generator.cpp index af71d04f6..028d271c4 100644 --- a/src/point_cloud/sdf_data_generator.cpp +++ b/src/point_cloud/sdf_data_generator.cpp @@ -40,6 +40,18 @@ SdfDataGenerator::SdfDataGenerator(const geometry::Points3& points, } } +SdfDataGenerator::SdfDataGenerator(const geometry::Points3& points, + const geometry::Vectors3& normals, double min_distance, + double max_distance, const Mat3& aniso) + : SdfDataGenerator(points, normals, max_distance, aniso) { + if (!(min_distance >= 0.0)) { + throw std::invalid_argument("min_distance must be non-negative"); + } + if (!(max_distance >= min_distance)) { + throw std::invalid_argument("max_distance must be greater than or equal to min_distance"); + } +} + std::pair SdfDataGenerator::estimate_impl( const geometry::Points3& points, const geometry::Vectors3& normals, double offset) { KdTree tree(points); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 108671f6d..7d224b019 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -29,6 +29,7 @@ add_executable(${TARGET} preconditioner/test_domain_divider.cpp preconditioner/test_fine_grid.cpp rbf/test_rbf.cpp + structural/test_structural_interpolant.cpp ) target_link_libraries(${TARGET} PRIVATE diff --git a/test/structural/test_structural_interpolant.cpp b/test/structural/test_structural_interpolant.cpp new file mode 100644 index 000000000..ff94684df --- /dev/null +++ b/test/structural/test_structural_interpolant.cpp @@ -0,0 +1,202 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using polatory::Index; +using polatory::Mat3; +using polatory::Model; +using polatory::VecX; +using polatory::geometry::Point3; +using polatory::geometry::Points3; +using polatory::rbf::CovSpheroidal3; +using polatory::structural::DomainSpec3; +using polatory::structural::StructuralDomainBuilder3; +using polatory::structural::StructuralInterpolant3; +using polatory::structural::StructuralTrendInput3; +using polatory::structural::StructuralTrendType; +using polatory::structural::TriangleFaces3; + +StructuralTrendInput3 horizontal_plane(double strength = 5.0, + double range = 10.0, + double z = 0.0) { + Points3 vertices(4, 3); + vertices << -1.0, -1.0, z, // + 1.0, -1.0, z, // + 1.0, 1.0, z, // + -1.0, 1.0, z; + + TriangleFaces3 faces(2, 3); + faces << 0, 1, 2, // + 0, 2, 3; + + return StructuralTrendInput3(vertices, faces, strength, range); +} + +StructuralTrendInput3 vertical_plane(double strength = 3.0, + double range = 10.0) { + Points3 vertices(4, 3); + vertices << 0.0, -1.0, -1.0, // + 0.0, 1.0, -1.0, // + 0.0, 1.0, 1.0, // + 0.0, -1.0, 1.0; + + TriangleFaces3 faces(2, 3); + faces << 0, 1, 2, // + 0, 2, 3; + + return StructuralTrendInput3(vertices, faces, strength, range); +} + +TEST(structural_interpolant, one_domain_matches_standard_interpolant) { + Points3 points(8, 3); + points << 0.0, 0.0, 0.0, // + 1.0, 0.0, 0.0, // + 0.0, 1.0, 0.0, // + 1.0, 1.0, 0.0, // + 0.0, 0.0, 1.0, // + 1.0, 0.0, 1.0, // + 0.0, 1.0, 1.0, // + 1.0, 1.0, 1.0; + + VecX values(8); + values << -1.0, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0; + + Mat3 anisotropy; + anisotropy << 0.8, 0.0, 0.0, // + 0.0, 1.2, 0.0, // + 0.0, 0.0, 1.0; + + CovSpheroidal3<3> standard_rbf({1.0, 4.0}); + standard_rbf.set_anisotropy(anisotropy); + Model<3> standard_model(std::move(standard_rbf), 0); + standard_model.set_nugget(0.01); + + polatory::Interpolant<3> standard(standard_model); + standard.fit(points, values, 1e-8, 100); + + CovSpheroidal3<3> structural_rbf({1.0, 4.0}); + Model<3> structural_model(std::move(structural_rbf), 0); + structural_model.set_nugget(0.01); + + Point3 bbox_min; + bbox_min << -1.0, -1.0, -1.0; + Point3 bbox_max; + bbox_max << 2.0, 2.0, 2.0; + + std::vector support_indices{0, 1, 2, 3, 4, 5, 6, 7}; + DomainSpec3 domain(anisotropy, bbox_min, bbox_max, support_indices); + + StructuralInterpolant3 structural(structural_model); + structural.fit(points, values, {domain}, 1e-8, 100); + + Points3 queries(4, 3); + queries << 0.2, 0.2, 0.2, // + 0.8, 0.3, 0.4, // + 0.4, 0.7, 0.9, // + 0.5, 0.5, 0.5; + + VecX expected = standard.evaluate(queries); + VecX actual = structural.evaluate(queries); + + EXPECT_LT((expected - actual).cwiseAbs().maxCoeff(), 1e-9); + EXPECT_EQ(structural.num_domains(), 1); +} + +TEST(structural_interpolant, returns_outside_value_outside_all_domain_boxes) { + Points3 points(4, 3); + points << 0.0, 0.0, 0.0, // + 1.0, 0.0, 0.0, // + 0.0, 1.0, 0.0, // + 0.0, 0.0, 1.0; + + VecX values(4); + values << -1.0, 0.0, 0.5, 1.0; + + CovSpheroidal3<3> rbf({1.0, 4.0}); + Model<3> model(std::move(rbf), 0); + + Point3 bbox_min; + bbox_min << -0.5, -0.5, -0.5; + Point3 bbox_max; + bbox_max << 1.5, 1.5, 1.5; + + DomainSpec3 domain(Mat3::Identity(), bbox_min, bbox_max, {0, 1, 2, 3}); + + StructuralInterpolant3 structural(model, -7.0); + structural.fit(points, values, {domain}, 1e-6, 100); + + Points3 query(1, 3); + query << 10.0, 10.0, 10.0; + + EXPECT_DOUBLE_EQ(structural.evaluate(query)(0), -7.0); +} + +TEST(structural_domain_builder, planar_input_matches_recovered_decay_formula) { + StructuralDomainBuilder3 builder; + auto input = horizontal_plane(); + + Points3 queries(2, 3); + queries << 1.0, 1.0, 0.0, // exactly on a mesh vertex + 1.0, 1.0, 10.0; // one range away + + auto samples = builder.sample( + queries, {input}, StructuralTrendType::kStrongestAlongInputs); + + EXPECT_NEAR(std::abs(samples.normals(0, 2)), 1.0, 1e-12); + EXPECT_NEAR(samples.ratios(0), 5.0, 1e-12); + EXPECT_NEAR(samples.ratios(1), 1.0 + 4.0 * std::exp(-1.0), 1e-12); + EXPECT_NEAR(samples.anisotropies.at(0).determinant(), 1.0, 1e-12); +} + +TEST(structural_domain_builder, strongest_mode_selects_largest_local_influence) { + StructuralDomainBuilder3 builder; + auto horizontal = horizontal_plane(5.0, 10.0); + auto vertical = vertical_plane(3.0, 10.0); + + Points3 query(1, 3); + query << 1.0, 1.0, 0.0; + + auto samples = builder.sample( + query, {horizontal, vertical}, + StructuralTrendType::kStrongestAlongInputs); + + EXPECT_EQ(samples.dominant_inputs.at(0), 0); + EXPECT_NEAR(std::abs(samples.normals(0, 2)), 1.0, 1e-12); +} + +TEST(structural_domain_builder, creates_overlapping_domains_from_mesh) { + Points3 points(27, 3); + Index row = 0; + for (Index x = 0; x < 3; ++x) { + for (Index y = 0; y < 3; ++y) { + for (Index z = 0; z < 3; ++z) { + points.row(row++) << static_cast(x), + static_cast(y), static_cast(z); + } + } + } + + StructuralDomainBuilder3 builder(1.5, 1.0, 4); + auto domains = builder.build( + points, {horizontal_plane()}, + StructuralTrendType::kStrongestAlongInputs); + + EXPECT_FALSE(domains.empty()); + for (const auto& domain : domains) { + EXPECT_GE(domain.support_indices().size(), 4u); + EXPECT_NEAR(domain.anisotropy().determinant(), 1.0, 1e-10); + } +} + +} // namespace