From 775fc025f378cbd07c171dc11990ed7b2aceac21 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:32:21 -0400 Subject: [PATCH 01/31] Add structural LVA domain specification --- include/polatory/structural/domain_spec.hpp | 48 +++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 include/polatory/structural/domain_spec.hpp 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 From 616f66d072a0357e264ad1a824107d0135d62af2 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:32:56 -0400 Subject: [PATCH 02/31] Add clustered structural LVA interpolant --- include/polatory/structural/interpolant.hpp | 220 ++++++++++++++++++++ 1 file changed, 220 insertions(+) create mode 100644 include/polatory/structural/interpolant.hpp diff --git a/include/polatory/structural/interpolant.hpp b/include/polatory/structural/interpolant.hpp new file mode 100644 index 000000000..261db0a25 --- /dev/null +++ b/include/polatory/structural/interpolant.hpp @@ -0,0 +1,220 @@ +#pragma once + +#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 Points = geometry::Points3; + using StandardInterpolant = Interpolant<3>; + + public: + explicit StructuralInterpolant3(const Model& base_model, double outside_value = -1.0, + double blend_power = 1.0) + : base_model_(base_model), outside_value_(outside_value), blend_power_(blend_power) { + if (!(blend_power_ > 0.0)) { + throw std::invalid_argument("blend_power must be positive"); + } + } + + const Bbox& bbox() const { + throw_if_not_fitted(); + return bbox_; + } + + double blend_power() const { return blend_power_; } + + Index num_domains() const { return static_cast(domains_.size()); } + + double outside_value() const { return outside_value_; } + + 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)}); + bbox_ = bbox_.is_empty() ? spec.bbox() : bbox_.convex_hull(spec.bbox()); + } + + 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); + 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 box_weight(const geometry::Point3& 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_); + } + + 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_; + bool fitted_{}; + std::vector domains_; + Bbox bbox_; +}; + +} // namespace polatory::structural From 5d64ff33fa8789194238b9f2336cb578fd2329b4 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:34:07 -0400 Subject: [PATCH 03/31] Add structural interpolant field-function adapter --- .../structural_rbf_field_function.hpp | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 include/polatory/isosurface/structural_rbf_field_function.hpp 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 From a6eb339c905fd1abed7e8398d516d28f1486fefb Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:34:19 -0400 Subject: [PATCH 04/31] Expose structural LVA interpolation headers --- include/polatory/polatory.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/polatory/polatory.hpp b/include/polatory/polatory.hpp index 4bf3717cd..e21dceae2 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,7 @@ #include #include #include +#include +#include #include #include From 6f6e15e27db7608790f2bd7594a70e785aea0baa Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:35:18 -0400 Subject: [PATCH 05/31] Add Python bindings for structural LVA interpolation --- python/structural_binding.cpp | 67 +++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 python/structural_binding.cpp diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp new file mode 100644 index 000000000..98f7069fa --- /dev/null +++ b/python/structural_binding.cpp @@ -0,0 +1,67 @@ +#include +#undef _GNU_SOURCE +#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 DomainSpec = structural::DomainSpec3; + using StructuralInterpolant = structural::StructuralInterpolant3; + + 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, "StructuralInterpolant3") + .def(py::init&, double, double>(), "base_model"_a, + "outside_value"_a = -1.0, "blend_power"_a = 1.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("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("evaluate", &StructuralInterpolant::evaluate, "points"_a, + "accuracy"_a = kInfinity); + + py::class_( + m, "StructuralRbfFieldFunction") + .def(py::init(), "interpolant"_a, + "accuracy"_a = kInfinity); +} From 8c2d4f4e36390a3e35b0a491050e1b231042e20d Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:35:26 -0400 Subject: [PATCH 06/31] Build structural LVA Python extension --- python/CMakeLists.txt | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) 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}) From 95560421a7e0e989daf4604881ff16e1a5be0ec8 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:35:32 -0400 Subject: [PATCH 07/31] Export structural LVA Python classes --- python/src/polatory/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/src/polatory/__init__.py b/python/src/polatory/__init__.py index 72aa94182..d26b31b86 100644 --- a/python/src/polatory/__init__.py +++ b/python/src/polatory/__init__.py @@ -1,2 +1,3 @@ from ._core import * from ._core import __doc__, __version__ +from ._structural import * From c33c1d3139393b12b570621f6d014cae0bd7068f Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:35:56 -0400 Subject: [PATCH 08/31] Add structural interpolant regression tests --- .../test_structural_interpolant.cpp | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 test/structural/test_structural_interpolant.cpp diff --git a/test/structural/test_structural_interpolant.cpp b/test/structural/test_structural_interpolant.cpp new file mode 100644 index 000000000..b068b9a5b --- /dev/null +++ b/test/structural/test_structural_interpolant.cpp @@ -0,0 +1,108 @@ +#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::StructuralInterpolant3; + +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); +} + +} // namespace From ac08aadb8db1e0aa88768c1a821675761bff798c Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:36:07 -0400 Subject: [PATCH 09/31] Register structural interpolant unit tests --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) 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 From 8017abfaed7be3df6a1fa0f6b822cc28b352257e Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 21:51:23 -0400 Subject: [PATCH 10/31] Use installed MSVC toolset on Windows --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 = {} From c8f37cf8f9f9689c2b2e3e66fb6e5fa0b15359cd Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 22:51:13 -0400 Subject: [PATCH 11/31] Fix SdfDataGenerator Python binding compatibility --- include/polatory/point_cloud/sdf_data_generator.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) 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_; From 4b994d3b50e8837024c768a719970333dc35b464 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Mon, 13 Jul 2026 22:51:32 -0400 Subject: [PATCH 12/31] Implement SdfDataGenerator compatibility overload --- src/point_cloud/sdf_data_generator.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) 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); From 908a702226968b7e63f59402c5fba332e9ad921f Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:19:49 -0400 Subject: [PATCH 13/31] Add automatic structural trend domain builder --- .../polatory/structural/domain_builder.hpp | 484 ++++++++++++++++++ 1 file changed, 484 insertions(+) create mode 100644 include/polatory/structural/domain_builder.hpp diff --git a/include/polatory/structural/domain_builder.hpp b/include/polatory/structural/domain_builder.hpp new file mode 100644 index 000000000..173416b6a --- /dev/null +++ b/include/polatory/structural/domain_builder.hpp @@ -0,0 +1,484 @@ +#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)) { + domain_size = 2.0 * max_range; + if (!(domain_size > 0.0)) { + domain_size = std::max({widths(0), widths(1), widths(2)}) / 3.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; + } + + using Cell = std::array; + std::map> cells; + + for (Index i = 0; i < points.rows(); ++i) { + Cell cell{}; + for (Index axis = 0; axis < 3; ++axis) { + 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) { + 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 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 From b982bb59227085f9680677879e2ddee682f698bf Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:20:03 -0400 Subject: [PATCH 14/31] Expose structural domain builder in public API --- include/polatory/polatory.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/polatory/polatory.hpp b/include/polatory/polatory.hpp index e21dceae2..010ca2223 100644 --- a/include/polatory/polatory.hpp +++ b/include/polatory/polatory.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include From fee006df6fceee74c0aafa2c5c8a842001918a48 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:20:31 -0400 Subject: [PATCH 15/31] Expose automatic mesh-driven structural trend builder --- python/structural_binding.cpp | 82 +++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index 98f7069fa..ead3d23e5 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -19,8 +19,43 @@ PYBIND11_MODULE(_structural, m) { // Register the standard Polatory types (Model<3>, FieldFunction, etc.) first. py::module_::import("polatory._core"); + 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, @@ -43,6 +78,25 @@ PYBIND11_MODULE(_structural, m) { .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("sample", &DomainBuilder::sample, + "query_points"_a, + "inputs"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs); + py::class_(m, "StructuralInterpolant3") .def(py::init&, double, double>(), "base_model"_a, "outside_value"_a = -1.0, "blend_power"_a = 1.0) @@ -57,6 +111,34 @@ PYBIND11_MODULE(_structural, m) { .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, + TrendType trend_type, + double tolerance, + 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, + "trend_type"_a = TrendType::kStrongestAlongInputs, + "tolerance"_a, + "max_iter"_a = 100, + "accuracy"_a = kInfinity, + "domain_size"_a = 0.0, + "overlap"_a = 0.0, + "min_support_points"_a = 4) .def("evaluate", &StructuralInterpolant::evaluate, "points"_a, "accuracy"_a = kInfinity); From a5c1f6861fd7bf502c0830e754f97bcd9ca41f18 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:21:27 -0400 Subject: [PATCH 16/31] Fix fit_from_meshes Python argument ordering --- python/structural_binding.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index ead3d23e5..10d3925a5 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -117,8 +117,8 @@ PYBIND11_MODULE(_structural, m) { const geometry::Points3& points, const VecX& values, const std::vector& inputs, - TrendType trend_type, double tolerance, + TrendType trend_type, int max_iter, double accuracy, double domain_size, @@ -132,8 +132,8 @@ PYBIND11_MODULE(_structural, m) { "points"_a, "values"_a, "inputs"_a, - "trend_type"_a = TrendType::kStrongestAlongInputs, "tolerance"_a, + "trend_type"_a = TrendType::kStrongestAlongInputs, "max_iter"_a = 100, "accuracy"_a = kInfinity, "domain_size"_a = 0.0, From 361ca70b991127206f1688805ee39752c6f3f20c Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:22:02 -0400 Subject: [PATCH 17/31] Test automatic structural trend domain builder --- .../test_structural_interpolant.cpp | 94 +++++++++++++++++++ 1 file changed, 94 insertions(+) diff --git a/test/structural/test_structural_interpolant.cpp b/test/structural/test_structural_interpolant.cpp index b068b9a5b..ff94684df 100644 --- a/test/structural/test_structural_interpolant.cpp +++ b/test/structural/test_structural_interpolant.cpp @@ -1,9 +1,11 @@ #include +#include #include #include #include #include +#include #include #include #include @@ -19,7 +21,42 @@ 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); @@ -105,4 +142,61 @@ TEST(structural_interpolant, returns_outside_value_outside_all_domain_boxes) { 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 From 73202fa99a533d339aff170c6ce7e1b328f3cff2 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:22:45 -0400 Subject: [PATCH 18/31] Document mesh-driven structural LVA API and parity limits --- docs/structural_lva.md | 103 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 docs/structural_lva.md 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. From c2d23c923485a76d0758a644321dfaf68cfecb6f Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 15:44:38 -0400 Subject: [PATCH 19/31] Improve automatic structural domain layout and continuity --- .../polatory/structural/domain_builder.hpp | 117 ++++++++++++++++-- 1 file changed, 106 insertions(+), 11 deletions(-) diff --git a/include/polatory/structural/domain_builder.hpp b/include/polatory/structural/domain_builder.hpp index 173416b6a..52a1a4168 100644 --- a/include/polatory/structural/domain_builder.hpp +++ b/include/polatory/structural/domain_builder.hpp @@ -127,9 +127,13 @@ class StructuralDomainBuilder3 { auto domain_size = domain_size_; if (!(domain_size > 0.0)) { - domain_size = 2.0 * max_range; + // 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)}) / 3.0; + domain_size = std::max({widths(0), widths(1), widths(2)}) / 4.0; } if (!(domain_size > 0.0)) { domain_size = 1.0; @@ -143,12 +147,24 @@ class StructuralDomainBuilder3 { : 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)); @@ -165,8 +181,14 @@ class StructuralDomainBuilder3 { Point sample_point = Point::Zero(); for (Index axis = 0; axis < 3; ++axis) { - auto cell_i = static_cast( - cell.at(static_cast(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); @@ -214,8 +236,8 @@ class StructuralDomainBuilder3 { 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); + 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; @@ -311,6 +333,79 @@ class StructuralDomainBuilder3 { 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()) { @@ -346,9 +441,9 @@ class StructuralDomainBuilder3 { auto q = trend_type == StructuralTrendType::kNonDecaying ? 1.0 - : std::exp(-distance / inputs.at( - static_cast(input_i)) - .range()); + : std::exp(-distance / + inputs.at(static_cast(input_i)) + .range()); auto contribution = (inputs.at(static_cast(input_i)).strength() - 1.0) * q; @@ -366,8 +461,8 @@ class StructuralDomainBuilder3 { } Point normal = normals.at(static_cast(dominant_input)); - auto ratio = 1.0 + - contributions.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(); From ec16a91e7e8c9aa85e74144de21ddc2ff6453187 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:39:55 -0400 Subject: [PATCH 20/31] Add axial domain orientation averaging experiment --- .../structural/domain_orientation_average.hpp | 250 ++++++++++++++++++ 1 file changed, 250 insertions(+) create mode 100644 include/polatory/structural/domain_orientation_average.hpp 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 From 0ecd43076a09acd53026d752ae84ebfbb942f4f0 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:40:37 -0400 Subject: [PATCH 21/31] Use axial orientation averaging for Python structural domains --- python/structural_binding.cpp | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index 10d3925a5..a9f422688 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include namespace py = pybind11; @@ -87,7 +88,23 @@ PYBIND11_MODULE(_structural, m) { .def_property_readonly("overlap", &DomainBuilder::overlap) .def_property_readonly("min_support_points", &DomainBuilder::min_support_points) - .def("build", &DomainBuilder::build, + .def( + "build", + [](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("build_center_sampled", &DomainBuilder::build, "points"_a, "inputs"_a, "trend_type"_a = TrendType::kStrongestAlongInputs, @@ -124,8 +141,9 @@ PYBIND11_MODULE(_structural, m) { 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); + auto domains = structural::build_orientation_averaged_domains( + points, inputs, trend_type, {}, domain_size, overlap, + min_support_points); interpolant.fit(points, values, domains, tolerance, max_iter, accuracy); return domains; }, From 28d3cbfeb6aa5eb265a7421f289c24d028ea85bb Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:52:16 -0400 Subject: [PATCH 22/31] Restore center-sampled domains as default after averaging test --- python/structural_binding.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index a9f422688..58c88e5df 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -88,8 +88,13 @@ PYBIND11_MODULE(_structural, m) { .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", + "build_orientation_averaged", [](const DomainBuilder& builder, const geometry::Points3& points, const std::vector& inputs, @@ -104,11 +109,6 @@ PYBIND11_MODULE(_structural, m) { "inputs"_a, "trend_type"_a = TrendType::kStrongestAlongInputs, "model_parameters"_a = std::vector{}) - .def("build_center_sampled", &DomainBuilder::build, - "points"_a, - "inputs"_a, - "trend_type"_a = TrendType::kStrongestAlongInputs, - "model_parameters"_a = std::vector{}) .def("sample", &DomainBuilder::sample, "query_points"_a, "inputs"_a, @@ -141,9 +141,8 @@ PYBIND11_MODULE(_structural, m) { double domain_size, double overlap, Index min_support_points) { - auto domains = structural::build_orientation_averaged_domains( - points, inputs, trend_type, {}, domain_size, overlap, - 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; }, From a318be5672ea137a370dfa8f3fd598c0399f8ba9 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 19:20:00 -0400 Subject: [PATCH 23/31] Add adaptive orientation-consistency structural domain builder --- .../structural/adaptive_domain_builder.hpp | 650 ++++++++++++++++++ 1 file changed, 650 insertions(+) create mode 100644 include/polatory/structural/adaptive_domain_builder.hpp 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 From c306ec22af401b383226bb34f3a39ca255b29464 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 19:20:52 -0400 Subject: [PATCH 24/31] Expose adaptive structural domain builder to Python --- python/structural_binding.cpp | 89 ++++++++++++++++++++++++++++++++--- 1 file changed, 83 insertions(+), 6 deletions(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index 58c88e5df..5250228e5 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -20,6 +21,7 @@ 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; @@ -60,14 +62,17 @@ PYBIND11_MODULE(_structural, m) { py::class_(m, "StructuralDomain3") .def(py::init([](const Mat3& anisotropy, const Eigen::Vector3d& bbox_min, - const Eigen::Vector3d& bbox_max, std::vector support_indices, + 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), + 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, + "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) { @@ -114,6 +119,36 @@ PYBIND11_MODULE(_structural, m) { "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{}); + py::class_(m, "StructuralInterpolant3") .def(py::init&, double, double>(), "base_model"_a, "outside_value"_a = -1.0, "blend_power"_a = 1.0) @@ -126,8 +161,9 @@ PYBIND11_MODULE(_structural, m) { .def_property_readonly("blend_power", &StructuralInterpolant::blend_power) .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", &StructuralInterpolant::fit, "points"_a, "values"_a, + "domains"_a, "tolerance"_a, "max_iter"_a = 100, + "accuracy"_a = kInfinity) .def( "fit_from_meshes", [](StructuralInterpolant& interpolant, @@ -143,7 +179,8 @@ PYBIND11_MODULE(_structural, m) { 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); + interpolant.fit(points, values, domains, tolerance, max_iter, + accuracy); return domains; }, "points"_a, @@ -156,6 +193,46 @@ PYBIND11_MODULE(_structural, m) { "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); From 2d833fde3819d56809ef4b767f5d471c1055e574 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 20:48:45 -0400 Subject: [PATCH 25/31] Add orientation averaging for adaptive structural domains --- .../adaptive_domain_orientation_average.hpp | 133 ++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 include/polatory/structural/adaptive_domain_orientation_average.hpp 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 From f61fdc7987eedbc39d8cb3e6d1ca592fd83cc3d2 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 20:49:40 -0400 Subject: [PATCH 26/31] Expose orientation-averaged adaptive domain build --- python/structural_binding.cpp | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index 5250228e5..300220232 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -147,7 +148,25 @@ PYBIND11_MODULE(_structural, m) { "points"_a, "inputs"_a, "trend_type"_a = TrendType::kStrongestAlongInputs, - "model_parameters"_a = std::vector{}); + "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>(), "base_model"_a, From 6fe6fbc7e1871c01101aa6f9dac9e621e62854af Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 22:14:07 -0400 Subject: [PATCH 27/31] Add optional overlap offset alignment to structural interpolant --- include/polatory/structural/interpolant.hpp | 194 ++++++++++++++++++-- 1 file changed, 176 insertions(+), 18 deletions(-) diff --git a/include/polatory/structural/interpolant.hpp b/include/polatory/structural/interpolant.hpp index 261db0a25..012ae212a 100644 --- a/include/polatory/structural/interpolant.hpp +++ b/include/polatory/structural/interpolant.hpp @@ -1,6 +1,9 @@ #pragma once +#include +#include #include +#include #include #include #include @@ -20,16 +23,25 @@ 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) - : base_model_(base_model), outside_value_(outside_value), blend_power_(blend_power) { + 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 { @@ -39,12 +51,26 @@ class StructuralInterpolant3 { 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_; } - void fit(const Points& points, const VecX& values, const std::vector& domain_specs, - double tolerance, int max_iter = 100, double accuracy = kInfinity) { + 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()"); } @@ -71,7 +97,8 @@ class StructuralInterpolant3 { 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) { + 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"); @@ -88,11 +115,19 @@ class StructuralInterpolant3 { rbf.set_anisotropy(spec.anisotropy()); } - auto local_interpolant = std::make_unique(local_model); - local_interpolant->fit(local_points, local_values, tolerance, max_iter, accuracy); + 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()); + } - domains_.push_back(Domain{spec, std::move(local_interpolant)}); - 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; @@ -136,12 +171,16 @@ class StructuralInterpolant3 { } 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))); + 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); - for (Index i = 0; i < static_cast(active_indices.size()); ++i) { + 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); @@ -158,7 +197,8 @@ class StructuralInterpolant3 { return result; } - void set_evaluation_bbox_impl(const Bbox& bbox, double accuracy = kInfinity) { + 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"); @@ -173,9 +213,10 @@ class StructuralInterpolant3 { struct Domain { DomainSpec3 spec; std::unique_ptr interpolant; + double offset; }; - double box_weight(const geometry::Point3& point, const Bbox& bbox) const { + double box_weight(const Point& point, const Bbox& bbox) const { if (!bbox.contains(point)) { return 0.0; } @@ -188,8 +229,10 @@ class StructuralInterpolant3 { } 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); + 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; } @@ -197,6 +240,119 @@ class StructuralInterpolant3 { 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.2, 0.5, 0.8}; + Points samples(27, 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); + + 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); + + // Give the zero-level neighbourhood more influence while retaining + // information from the complete overlap volume. + auto level_weight = + 1.0 / (1.0 + std::abs(predictions_i(sample_i)) + + std::abs(predictions_j(sample_i))); + 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(); @@ -205,13 +361,15 @@ class StructuralInterpolant3 { void throw_if_not_fitted() const { if (!fitted_) { - throw std::runtime_error("structural interpolant has not been 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_; From 4fb7a8160f97e912a2bcb475cc118d1d9b58dbf1 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Tue, 14 Jul 2026 22:15:22 -0400 Subject: [PATCH 28/31] Expose structural overlap alignment controls --- python/structural_binding.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/python/structural_binding.cpp b/python/structural_binding.cpp index 300220232..6ac54e3ff 100644 --- a/python/structural_binding.cpp +++ b/python/structural_binding.cpp @@ -169,8 +169,11 @@ PYBIND11_MODULE(_structural, m) { "model_parameters"_a = std::vector{}); py::class_(m, "StructuralInterpolant3") - .def(py::init&, double, double>(), "base_model"_a, - "outside_value"_a = -1.0, "blend_power"_a = 1.0) + .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()); }) @@ -178,6 +181,10 @@ PYBIND11_MODULE(_structural, m) { 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, From 1c731321db716f9b0c9ca8c8057dc00d59b10de4 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Wed, 15 Jul 2026 07:21:41 -0400 Subject: [PATCH 29/31] Refine overlap alignment near zero level --- include/polatory/structural/interpolant.hpp | 33 +++++++++++++++++---- 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/include/polatory/structural/interpolant.hpp b/include/polatory/structural/interpolant.hpp index 012ae212a..7323861e8 100644 --- a/include/polatory/structural/interpolant.hpp +++ b/include/polatory/structural/interpolant.hpp @@ -249,8 +249,15 @@ class StructuralInterpolant3 { static Points make_overlap_samples(const Point& overlap_min, const Point& overlap_max) { - constexpr std::array fractions{0.2, 0.5, 0.8}; - Points samples(27, 3); + 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) { @@ -290,6 +297,19 @@ class StructuralInterpolant3 { 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) { @@ -301,11 +321,12 @@ class StructuralInterpolant3 { domains_.at(static_cast(j)).spec.bbox()); auto overlap_weight = std::sqrt(wi * wj); - // Give the zero-level neighbourhood more influence while retaining - // information from the complete overlap volume. + 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 + std::abs(predictions_i(sample_i)) + - std::abs(predictions_j(sample_i))); + 1.0 / (1.0 + normalized_level * normalized_level); auto weight = overlap_weight * level_weight; weighted_difference += weight * (predictions_j(sample_i) - From 2f1131a51617c73bc890351c5a0ca864a1cb9341 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Wed, 15 Jul 2026 22:43:51 -0400 Subject: [PATCH 30/31] Add validated clustered structural domain builder --- .../src/polatory/clustered_domain_builder.py | 440 ++++++++++++++++++ 1 file changed, 440 insertions(+) create mode 100644 python/src/polatory/clustered_domain_builder.py 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 From f0eee0766130df701a9ef3933cca7ea7ee17b575 Mon Sep 17 00:00:00 2001 From: Mining Geologist <127529341+Mining-Geologist@users.noreply.github.com> Date: Wed, 15 Jul 2026 22:44:09 -0400 Subject: [PATCH 31/31] Expose clustered structural domain builder --- python/src/polatory/__init__.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/python/src/polatory/__init__.py b/python/src/polatory/__init__.py index d26b31b86..62f001fea 100644 --- a/python/src/polatory/__init__.py +++ b/python/src/polatory/__init__.py @@ -1,3 +1,7 @@ from ._core import * from ._core import __doc__, __version__ from ._structural import * +from .clustered_domain_builder import ( + ClusteredStructuralDomainBuilder3, + fit_from_meshes_clustered, +)