From 553f43df06a61c759ea7ac96c5aaa8cf4b66850b Mon Sep 17 00:00:00 2001 From: graeme Date: Wed, 12 Aug 2026 11:46:24 -0400 Subject: [PATCH 1/8] Add ferreus POC test versus existing method --- surface_apps/ferreus/__init__.py | 8 + .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 1215 +++++++++++++++++ .../ferreus_rmt_real_geochem_comparison.ipynb | 707 ++++++++++ .../geochem_pts/geochem_pts_example.geoh5 | 3 + 4 files changed, 1933 insertions(+) create mode 100644 surface_apps/ferreus/__init__.py create mode 100644 surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb create mode 100644 surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb create mode 100644 surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 diff --git a/surface_apps/ferreus/__init__.py b/surface_apps/ferreus/__init__.py new file mode 100644 index 0000000..1b3743e --- /dev/null +++ b/surface_apps/ferreus/__init__.py @@ -0,0 +1,8 @@ +# '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of surface-apps package. ' +# ' +# surface-apps is distributed under the terms and conditions of the MIT License +# (see LICENSE file at the root of this source code package). ' +# '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb new file mode 100644 index 0000000..76de64c --- /dev/null +++ b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -0,0 +1,1215 @@ +{ + "cells": [ + { + "metadata": { + "collapsed": true, + "ExecuteTime": { + "end_time": "2026-08-10T15:19:45.243382800Z", + "start_time": "2026-08-10T15:19:45.191472Z" + } + }, + "cell_type": "markdown", + "source": [ + "# Ferreus RMT to geoh5 proof of concept\n", + "\n", + "Objectives:\n", + "\n", + "Explore if the packages allows for\n", + "\n", + " exclusion constraints\n", + "\n", + " orientation constraints\n", + "\n", + "POC with notebook creating a surface transferred to geoh5" + ], + "id": "52652cf59cffb624" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "", + "id": "de308ccc50aea28d" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "import os\n", + "import sys\n", + "from pathlib import Path\n", + "\n", + "from datetime import datetime\n", + "\n", + "from geoh5py.groups import ContainerGroup\n", + "from geoh5py.objects import Points, Surface\n", + "from geoh5py.workspace import Workspace\n", + "import numpy as np\n", + "\n", + "from ferreus_rmt import (\n", + " BoundaryClosure,\n", + " ClusterMethod,\n", + " build_isosurface,\n", + ")\n", + "from ferreus_rbf import RBFInterpolator\n", + "from ferreus_rbf.interpolant_config import (\n", + " Drift,\n", + " InterpolantSettings,\n", + " RBFKernelType,\n", + ")\n", + "import matplotlib.pyplot as plt" + ], + "id": "c14d7389d994a354", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### define a known sphere", + "id": "a8c8dc3ebfbae12c" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "\n", + "CENTER = np.array([0.0, 0.0, 0.0], dtype=np.float64)\n", + "RADIUS = 50.0\n", + "\n", + "def sphere_values(points: np.ndarray) -> np.ndarray:\n", + " \"\"\"negative inside the sphere, zero on it, positive outside.\"\"\"\n", + "\n", + " offsets = points - CENTER\n", + " return np.linalg.norm(offsets, axis=1) - RADIUS\n", + "\n", + "def sphere_values_and_gradients(\n", + " points: np.ndarray,\n", + ") -> tuple[np.ndarray, np.ndarray]:\n", + "\n", + " \"\"\"return values and analytic gradients.\"\"\"\n", + " offsets = points - CENTER\n", + " distances = np.linalg.norm(offsets, axis=1)\n", + "\n", + " gradients = np.divide(\n", + " offsets,\n", + " distances[:, None],\n", + " out=np.zeros_like(offsets),\n", + " where=distances[:, None] > 0.0,\n", + " )\n", + "\n", + " return distances - RADIUS, gradients" + ], + "id": "9be78b4d9cd5c492", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### extract the surface with rmt", + "id": "bf48bdde9dab3801" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "\n", + "\n", + "\n", + "seed_points = np.array(\n", + " [\n", + " [50.0, 0.0, 0.0],\n", + " [-50.0, 0.0, 0.0],\n", + " ],\n", + " dtype=np.float64,\n", + ")\n", + "\n", + "# Define the axis-aligned bounding box extents to extract the isosurface within:\n", + "# [xmin, ymin, zmin, xmax, ymax, zmax].\n", + "extents = np.array(\n", + " [-75.0, -75.0, -75.0, 75.0, 75.0, 75.0],\n", + " dtype=np.float64,\n", + ")\n", + "\n", + "\n", + "\n", + "resolution = 5.0\n", + "isovalue = 0.0\n", + "\n", + "mesh = build_isosurface(\n", + " seed_points, # approx. locations of the surface to seed the wavefront\n", + " extents, # the axis-aligned bounding box limiting extraction\n", + " resolution, # lattice spacing in coordinate units - smaller means more detail\n", + " isovalue, # the scalar value defining the surface\n", + " surface_fn=sphere_values, # surface_fn (sphere_values) is the scalar function that defines the geometry\n", + " gradient_fn=sphere_values_and_gradients, #\n", + ")\n", + "\n", + "vertices = np.asarray(mesh.vertices)\n", + "cells = np.asarray(mesh.facets, dtype=np.uint32)\n", + "\n", + "print(\"Vertices:\", vertices.shape)\n", + "print(\"Triangles:\", cells.shape)" + ], + "id": "8d2fa33ef1334bfa", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### validate the result", + "id": "4178f0ab56928714" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "\n", + "radii = np.linalg.norm(vertices - CENTER, axis=1)\n", + "radius_errors = radii - RADIUS\n", + "\n", + "print(\"Minimum radius:\", radii.min())\n", + "print(\"Maximum radius:\", radii.max())\n", + "print(\"Mean absolute radius error:\", np.abs(radius_errors).mean())\n", + "print(\"Maximum absolute radius error:\", np.abs(radius_errors).max())\n", + "\n", + "assert cells.ndim == 2 and cells.shape[1] == 3\n", + "assert vertices.ndim == 2 and vertices.shape[1] == 3\n", + "assert np.abs(radius_errors).max() < resolution\n", + "\n", + "# numbers should be close to the radius, but likely not exact because of the resolution of the lattice." + ], + "id": "c953cf37f21ad3a0", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### show the sphere we've constructed", + "id": "3b1c992eb5ef149d" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "\n", + "figure = plt.figure(figsize=(8, 7))\n", + "axis = figure.add_subplot(111, projection=\"3d\")\n", + "\n", + "axis.plot_trisurf(\n", + " vertices[:, 0],\n", + " vertices[:, 1],\n", + " vertices[:, 2],\n", + " triangles=cells,\n", + " linewidth=0.05,\n", + " alpha=0.85,\n", + ")\n", + "\n", + "axis.set_xlabel(\"X\")\n", + "axis.set_ylabel(\"Y\")\n", + "axis.set_zlabel(\"Z\")\n", + "axis.set_box_aspect(np.ptp(vertices, axis=0))\n", + "axis.set_title(\"Analytic sphere extracted with Ferreus RMT\")\n", + "\n", + "plt.show()" + ], + "id": "8c91eaf424ec37fb", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Export to geoh5", + "id": "ba702ddc81228ad6" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "\n", + "\n", + "output_directory = Path.cwd() / \"ferreus_output\"\n", + "output_directory.mkdir(exist_ok=True)\n", + "\n", + "output_path = output_directory / (\n", + " f\"ferreus_rmt_poc_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", + ")\n", + "\n", + "with Workspace.create(output_path) as workspace:\n", + " group = ContainerGroup.create(\n", + " workspace,\n", + " name=\"Ferreus RMT POC\",\n", + " )\n", + "\n", + " Surface.create(\n", + " workspace,\n", + " name=\"RMT analytic sphere\",\n", + " vertices=vertices,\n", + " cells=cells,\n", + " parent=group,\n", + " )\n", + "\n", + "print(\"Created:\", output_path)" + ], + "id": "f39977f0e682fcb9", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Fit a synthetic folded horizon with ferreus_rbf\n", + "\n", + "The analytic sphere proved that `ferreus_rmt` can extract a known scalar\n", + "field and that its mesh can be exported to geoH5.\n", + "\n", + "This section replaces the analytic sphere field with an RBF field fitted\n", + "from scattered scalar observations.\n", + "\n", + "A known analytic geological field is retained as ground truth, allowing\n", + "the interpolation and resulting surface errors to be measured." + ], + "id": "1a8d7c5fad014c7f" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "MODEL_EXTENTS = np.array(\n", + " [-100.0, -100.0, -60.0, 100.0, 100.0, 60.0],\n", + " dtype=np.float64,\n", + ")\n", + "\n", + "\n", + "def true_horizon_elevation(\n", + " x: np.ndarray,\n", + " y: np.ndarray,\n", + ") -> np.ndarray:\n", + " \"\"\"Elevation of the synthetic folded horizon.\"\"\"\n", + " return (\n", + " 20.0 * np.sin(x / 35.0)\n", + " + 10.0 * np.cos(y / 40.0)\n", + " )\n", + "\n", + "\n", + "def true_geological_field(points: np.ndarray) -> np.ndarray:\n", + " \"\"\"Signed vertical separation from the synthetic horizon.\"\"\"\n", + " horizon_z = true_horizon_elevation(\n", + " points[:, 0],\n", + " points[:, 1],\n", + " )\n", + "\n", + " return points[:, 2] - horizon_z" + ], + "id": "e8a9c39ec7391034", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Generate scattered scalar observations", + "id": "55396a8e404a0f23" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "random_generator = np.random.default_rng(42)\n", + "\n", + "number_of_source_points = 400\n", + "\n", + "source_points = random_generator.uniform(\n", + " low=MODEL_EXTENTS[:3],\n", + " high=MODEL_EXTENTS[3:],\n", + " size=(number_of_source_points, 3),\n", + ").astype(np.float64)\n", + "\n", + "source_values = true_geological_field(source_points).astype(\n", + " np.float64\n", + ")\n", + "\n", + "print(\"Source points:\", source_points.shape)\n", + "print(\"Source values:\", source_values.shape)\n", + "print(\n", + " \"Value range:\",\n", + " float(source_values.min()),\n", + " \"to\",\n", + " float(source_values.max()),\n", + ")" + ], + "id": "63f6e9858be0976d", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Fit the RBF field", + "id": "203c7442a4c55fec" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "rbf_settings = InterpolantSettings(\n", + " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", + " drift=Drift.Linear, # adds a broad linear trend, useful here because the field contains a strong linear z component.\n", + " nugget=0.0, # 0.0 requests an unsmoothed fit. Fine for synthetic observations as they contain no noise.\n", + ")\n", + "\n", + "rbf_model = RBFInterpolator(\n", + " points=source_points,\n", + " values=source_values,\n", + " interpolant_settings=rbf_settings,\n", + ")" + ], + "id": "b82e4349d22b0c02", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Check the training residuals", + "id": "6b8e055979af129a" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "source_predictions = np.asarray(\n", + " rbf_model.evaluate_at_source(add_nugget=False)\n", + ").reshape(-1)\n", + "\n", + "source_residuals = source_predictions - source_values\n", + "\n", + "print(\n", + " \"Source RMSE:\",\n", + " float(np.sqrt(np.mean(source_residuals**2))),\n", + ")\n", + "print(\n", + " \"Maximum source error:\",\n", + " float(np.max(np.abs(source_residuals))),\n", + ")\n", + "\n", + "# This demonstrates how closely the solved RBF reproduces the values used to fit it.\n", + "# For a zero-nugget, synthetic problem, the residuals should be very small, on the order of machine precision." + ], + "id": "f1900adf375eec99", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Test at locations not used for fitting", + "id": "e25a69ffedc71d99" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "number_of_validation_points = 200\n", + "\n", + "validation_points = random_generator.uniform(\n", + " low=MODEL_EXTENTS[:3],\n", + " high=MODEL_EXTENTS[3:],\n", + " size=(number_of_validation_points, 3),\n", + ").astype(np.float64)\n", + "\n", + "true_validation_values = true_geological_field(\n", + " validation_points\n", + ")\n", + "\n", + "predicted_validation_values = np.asarray(\n", + " rbf_model.evaluate(validation_points)\n", + ").reshape(-1)\n", + "\n", + "validation_residuals = (\n", + " predicted_validation_values - true_validation_values\n", + ")\n", + "\n", + "validation_rmse = np.sqrt(\n", + " np.mean(validation_residuals**2)\n", + ")\n", + "\n", + "print(\"Validation RMSE:\", float(validation_rmse))\n", + "print(\n", + " \"Maximum validation error:\",\n", + " float(np.max(np.abs(validation_residuals))),\n", + ")\n", + "\n", + "# The source residual measures fitting accuracy, i.e. an RMSE of 0.69 means approx 0.69 vertical units\n", + "# The validation residual measures interpolation accuracy between observations" + ], + "id": "28d26151ceb1bbf9", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "### Prepare the RBF for repeated RMT evaluations", + "id": "adcfe7b8d0d585b4" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "folded_resolution = 5.0\n", + "\n", + "# RMT evaluates a small halo of lattice positions around the extraction\n", + "# box. Give the cached RBF evaluator enough room for those calls.\n", + "evaluator_padding = 5.0 * folded_resolution\n", + "\n", + "evaluator_extents = MODEL_EXTENTS + np.array(\n", + " [\n", + " -evaluator_padding,\n", + " -evaluator_padding,\n", + " -evaluator_padding,\n", + " evaluator_padding,\n", + " evaluator_padding,\n", + " evaluator_padding,\n", + " ],\n", + " dtype=np.float64,\n", + ")\n", + "\n", + "rbf_model.build_evaluator(evaluator_extents)" + ], + "id": "3c481cec6ea9a739", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Choose seed points for the folded horizon\n", + "\n", + "RMT is a surface-following algorithm. This means it does not scan every cell in the model volume. Instead, it starts near the requested isosurface at the seed points and follows the connected surface away from them.\n", + "\n", + "These seeds are algorithm inputs, not extra observations used to fit the RBF. In this synthetic example the true horizon is known, so a small grid of exact surface locations is convenient. In a real geological workflow, observed contact points with scalar value zero would be natural seeds.\n", + "\n", + "Multiple seeds are useful when a model may contain disconnected surface components. A component that has no nearby seed may not be discovered." + ], + "id": "e674a1b0a8c42db5" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "seed_coordinates = np.linspace(-80.0, 80.0, 5)\n", + "seed_x, seed_y = np.meshgrid(\n", + " seed_coordinates,\n", + " seed_coordinates,\n", + " indexing=\"xy\",\n", + ")\n", + "\n", + "folded_seed_points = np.column_stack(\n", + " [\n", + " seed_x.ravel(),\n", + " seed_y.ravel(),\n", + " true_horizon_elevation(\n", + " seed_x.ravel(),\n", + " seed_y.ravel(),\n", + " ),\n", + " ]\n", + ").astype(np.float64)\n", + "\n", + "print(\"RMT seed points shape:\", folded_seed_points.shape)\n", + "print(\"Maximum absolute seed field value:\", float(\n", + " np.max(np.abs(true_geological_field(folded_seed_points)))\n", + "))" + ], + "id": "16cbac4b3b5c9df8", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Extract the fitted zero-isosurface with RMT\n", + "\n", + "The requested surface is where the fitted scalar field equals folded_isovalue, here zero. Because the synthetic field is z - horizon_z, zero represents the folded horizon, negative values lie below it, and positive values lie above it.\n", + "\n", + "folded_resolution is the spacing of RMT's internal sampling lattice in the same coordinate units as X, Y, and Z. 5.0 means samples are approximately five units apart. A smaller value can preserve finer geometry, but requires more field evaluations and usually produces more triangles.\n", + "\n", + "evaluate_targets and evaluate_targets_with_gradients use the cached evaluator prepared above. The gradients help RMT project its working vertices toward the requested scalar value; they are not geological orientation observations." + ], + "id": "f25538fc35ed124b" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "folded_isovalue = 0.0\n", + "\n", + "folded_mesh = build_isosurface(\n", + " seed_points=folded_seed_points,\n", + " extents=MODEL_EXTENTS,\n", + " resolution=folded_resolution,\n", + " isovalue=folded_isovalue,\n", + " surface_fn=rbf_model.evaluate_targets,\n", + " gradient_fn=rbf_model.evaluate_targets_with_gradients,\n", + " cluster_method=ClusterMethod.CurvatureWeighted,\n", + " boundary_closure=BoundaryClosure.None_,\n", + ")\n", + "\n", + "folded_vertices = np.asarray(\n", + " folded_mesh.vertices,\n", + " dtype=np.float64,\n", + ")\n", + "folded_cells = np.asarray(\n", + " folded_mesh.facets,\n", + " dtype=np.uint32,\n", + ")\n", + "\n", + "print(\"Folded surface vertices:\", folded_vertices.shape)\n", + "print(\"Folded surface triangles:\", folded_cells.shape)" + ], + "id": "7a79146bd28e5dd1", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Validate both parts of the workflow\n", + "\n", + "There are two distinct errors worth checking:\n", + "\n", + "1. **RMT extraction error:** evaluate the fitted RBF at every mesh vertex. These values should be close to the requested isovalue of zero. This checks whether RMT followed the field it was given.\n", + "2. **Modelling error:** evaluate the known analytic field at every mesh vertex. Here that value is also the signed vertical distance from the true horizon. This checks how closely the RBF-derived surface matches the ground truth between the observations.\n", + "\n", + "Keeping these checks separate matters as a mesh can follow an imperfect RBF field very accurately." + ], + "id": "e203369f6da5b01f" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "fitted_values_at_vertices = np.asarray(\n", + " rbf_model.evaluate(folded_vertices),\n", + " dtype=np.float64,\n", + ").reshape(-1)\n", + "\n", + "rmt_isovalue_errors = (\n", + " fitted_values_at_vertices - folded_isovalue\n", + ")\n", + "\n", + "# For this particular synthetic field, this is the signed vertical\n", + "# separation between a mesh vertex and the true folded horizon.\n", + "true_vertical_errors = true_geological_field(folded_vertices)\n", + "absolute_vertical_errors = np.abs(true_vertical_errors)\n", + "\n", + "print(\n", + " \"Maximum absolute RMT isovalue error:\",\n", + " float(np.max(np.abs(rmt_isovalue_errors))),\n", + ")\n", + "print(\n", + " \"Mean absolute vertical surface error:\",\n", + " float(np.mean(absolute_vertical_errors)),\n", + ")\n", + "print(\n", + " \"Vertical error percentiles (50%, 90%, 95%, 100%):\",\n", + " np.percentile(absolute_vertical_errors, [50, 90, 95, 100]),\n", + ")\n", + "\n", + "assert folded_vertices.ndim == 2 and folded_vertices.shape[1] == 3\n", + "assert folded_cells.ndim == 2 and folded_cells.shape[1] == 3\n", + "assert np.all(np.isfinite(folded_vertices))\n", + "assert np.all(np.isfinite(fitted_values_at_vertices))" + ], + "id": "798b9ee4fbcc148e", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Inspect the fitted surface and nearby observations\n", + "\n", + "Only observations reasonably close to the horizon are plotted, which keeps the diagnostic readable. Their colours show the signed scalar values supplied to the RBF: blue points are below the horizon, red points are above it, and values near zero lie close to the interpreted contact." + ], + "id": "be32da3dd2532365" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "near_horizon = np.abs(source_values) <= 15.0\n", + "\n", + "figure = plt.figure(figsize=(10, 8))\n", + "axis = figure.add_subplot(111, projection=\"3d\")\n", + "\n", + "axis.plot_trisurf(\n", + " folded_vertices[:, 0],\n", + " folded_vertices[:, 1],\n", + " folded_vertices[:, 2],\n", + " triangles=folded_cells,\n", + " color=\"lightsteelblue\",\n", + " linewidth=0.05,\n", + " alpha=0.70,\n", + ")\n", + "\n", + "observation_plot = axis.scatter(\n", + " source_points[near_horizon, 0],\n", + " source_points[near_horizon, 1],\n", + " source_points[near_horizon, 2],\n", + " c=source_values[near_horizon],\n", + " cmap=\"coolwarm\",\n", + " vmin=-15.0,\n", + " vmax=15.0,\n", + " s=20,\n", + " edgecolors=\"black\",\n", + " linewidths=0.2,\n", + ")\n", + "\n", + "figure.colorbar(\n", + " observation_plot,\n", + " ax=axis,\n", + " shrink=0.65,\n", + " pad=0.10,\n", + " label=\"Observed signed field value\",\n", + ")\n", + "\n", + "axis.set_xlabel(\"X\")\n", + "axis.set_ylabel(\"Y\")\n", + "axis.set_zlabel(\"Z\")\n", + "axis.set_box_aspect(MODEL_EXTENTS[3:] - MODEL_EXTENTS[:3])\n", + "axis.view_init(elev=25, azim=-55)\n", + "axis.set_title(\"RBF folded horizon extracted with Ferreus RMT\")\n", + "\n", + "plt.show()" + ], + "id": "40448b38d7a150c4", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Export the folded surface and its observations to geoh5\n", + "\n", + "The output contains:\n", + "\n", + "- a Points object holding every scalar observation used to fit the RBF;\n", + "- the interpolated Surface object produced by RMT;\n", + "- the fitted RBF value on each surface vertex, useful for checking extraction; and\n", + "- the signed vertical error relative to the known synthetic horizon, useful only because this POC has ground truth.\n", + "\n", + "In Geoscience ANALYST, colour the observation points by **Implicit scalar value** and the surface by **True vertical error** to inspect the result." + ], + "id": "9eab773cbf3c33cd" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "output_directory = Path.cwd() / \"ferreus_output\"\n", + "output_directory.mkdir(exist_ok=True)\n", + "\n", + "folded_output_path = output_directory / (\n", + " f\"ferreus_rbf_folded_horizon_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", + ")\n", + "\n", + "with Workspace.create(folded_output_path) as workspace:\n", + " group = ContainerGroup.create(\n", + " workspace,\n", + " name=\"Ferreus folded horizon POC\",\n", + " )\n", + "\n", + " observations = Points.create(\n", + " workspace,\n", + " name=\"RBF scalar observations\",\n", + " vertices=source_points,\n", + " parent=group,\n", + " )\n", + " observations.add_data(\n", + " {\n", + " \"Implicit scalar value\": {\n", + " \"values\": source_values,\n", + " }\n", + " }\n", + " )\n", + "\n", + " folded_surface = Surface.create(\n", + " workspace,\n", + " name=\"RBF folded horizon - isovalue 0\",\n", + " vertices=folded_vertices,\n", + " cells=folded_cells,\n", + " parent=group,\n", + " )\n", + " folded_surface.add_data(\n", + " {\n", + " \"RBF value on surface\": {\n", + " \"values\": fitted_values_at_vertices,\n", + " },\n", + " \"True vertical error\": {\n", + " \"values\": true_vertical_errors,\n", + " },\n", + " }\n", + " )\n", + "\n", + "print(\"Created:\", folded_output_path)" + ], + "id": "156ec95b0d1e1877", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## What this baseline proves\n", + "\n", + "This section demonstrates the complete technical path: scattered scalar observations -> ferreus_rbf scalar field -> ferreus_rmt zero-isosurface -> geoh5 objects for ANALYST.\n", + "\n", + "It isn't yet a test of geological orientation or exclusion constraints. Every source point above was assigned a full scalar value from the known analytic field. Real geological inputs more commonly include:\n", + "\n", + "- **contact/interface points**, which say f(p) = 0;\n", + "- **orientation measurements**, which say that the gradient of f should point along a measured normal; and\n", + "- **exclusion/inequality observations**, which specify only that a point must be on one side of the surface, such as f(p) > 0 or f(p) < 0.\n", + "\n", + "The public ferreus_rbf Python constructor accepts point coordinates and scalar values; it does not expose separate derivative or inequality observations. The next POC sections can therefore test value-based encodings:\n", + "\n", + "1. approximate an orientation at contact point p with offset samples p + delta times n and p - delta times n, assigned positive and negative values respectively; and\n", + "2. approximate an exclusion constraint with a signed pseudo-value on the required side of the surface.\n", + "\n", + "Those are approximations that should be measured against this baseline." + ], + "id": "f74bf0dad647895f" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Compare the existing marching-cubes method with Ferreus RMT\n", + "\n", + "This is an extraction-algorithm comparison. Both algorithms receive the same fitted RBF field, model bounds, lattice spacing, and isovalue. That keeps the scalar model fixed so differences in the resulting meshes come from marching cubes versus RMT rather than from different geological interpolation methods.\n", + "\n", + "The existing application first samples data onto a complete regular grid and then calls scikit-image marching cubes. The code below performs that same dense-grid extraction. RMT instead evaluates the field while following the surface outward from seed points.\n", + "\n", + "This comparison does not yet reproduce every part of the existing application's preprocessing, such as weighted averaging, maximum distance, or topographic masking. A later end-to-end test can feed both implementations the same real geoh5 object." + ], + "id": "23da01a0b8b33419" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "from time import perf_counter\n", + "\n", + "from skimage.measure import marching_cubes\n", + "\n", + "\n", + "# Construct the same kind of dense regular grid used by the existing\n", + "# iso_surfaces application. The array axes are X, Y, then Z.\n", + "comparison_grid = [\n", + " np.arange(\n", + " MODEL_EXTENTS[axis_index],\n", + " MODEL_EXTENTS[axis_index + 3] + 0.5 * folded_resolution,\n", + " folded_resolution,\n", + " dtype=np.float64,\n", + " )\n", + " for axis_index in range(3)\n", + "]\n", + "\n", + "comparison_x, comparison_y, comparison_z = np.meshgrid(\n", + " comparison_grid[0],\n", + " comparison_grid[1],\n", + " comparison_grid[2],\n", + ")\n", + "comparison_grid_points = np.column_stack(\n", + " [\n", + " comparison_x.ravel(),\n", + " comparison_y.ravel(),\n", + " comparison_z.ravel(),\n", + " ]\n", + ")\n", + "\n", + "marching_cubes_total_start = perf_counter()\n", + "\n", + "grid_evaluation_start = perf_counter()\n", + "comparison_grid_values = np.asarray(\n", + " rbf_model.evaluate_targets(comparison_grid_points),\n", + " dtype=np.float64,\n", + ").reshape(comparison_x.shape)\n", + "grid_evaluation_seconds = perf_counter() - grid_evaluation_start\n", + "\n", + "marching_cubes_start = perf_counter()\n", + "mc_index_vertices, mc_cells, _mc_normals, _mc_values = marching_cubes(\n", + " comparison_grid_values,\n", + " level=folded_isovalue,\n", + ")\n", + "\n", + "# marching_cubes returns positions in array-index coordinates. Map each\n", + "# axis back to the original X, Y, and Z coordinates, matching the logic\n", + "# in surface_apps.iso_surfaces.utils.extract_iso_surfaces.\n", + "mc_vertices = np.column_stack(\n", + " [\n", + " np.interp(\n", + " mc_index_vertices[:, axis_index],\n", + " np.arange(comparison_grid[axis_index].size),\n", + " comparison_grid[axis_index],\n", + " )\n", + " for axis_index in range(3)\n", + " ]\n", + ").astype(np.float64)\n", + "mc_cells = np.asarray(mc_cells, dtype=np.uint32)\n", + "\n", + "marching_cubes_algorithm_seconds = (\n", + " perf_counter() - marching_cubes_start\n", + ")\n", + "marching_cubes_total_seconds = (\n", + " perf_counter() - marching_cubes_total_start\n", + ")\n", + "\n", + "print(\"Regular-grid shape:\", comparison_grid_values.shape)\n", + "print(\"Regular-grid sample count:\", comparison_grid_points.shape[0])\n", + "print(\"Marching-cubes vertices:\", mc_vertices.shape)\n", + "print(\"Marching-cubes triangles:\", mc_cells.shape)\n", + "print(\"Grid evaluation seconds:\", grid_evaluation_seconds)\n", + "print(\"Marching-cubes algorithm seconds:\", marching_cubes_algorithm_seconds)\n", + "print(\"Marching-cubes total seconds:\", marching_cubes_total_seconds)" + ], + "id": "9117b165c6dfa0e7", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Time RMT against the same field\n", + "\n", + "The RMT timing includes its field evaluations and mesh construction. The RBF evaluator is already cached for both methods. These single-run timings are useful for orientation, but they are not a rigorous performance benchmark; repeated warm runs and larger datasets would be needed for that." + ], + "id": "b479c5ef098030d6" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "rmt_total_start = perf_counter()\n", + "\n", + "comparison_rmt_mesh = build_isosurface(\n", + " seed_points=folded_seed_points,\n", + " extents=MODEL_EXTENTS,\n", + " resolution=folded_resolution,\n", + " isovalue=folded_isovalue,\n", + " surface_fn=rbf_model.evaluate_targets,\n", + " gradient_fn=rbf_model.evaluate_targets_with_gradients,\n", + " cluster_method=ClusterMethod.CurvatureWeighted,\n", + " boundary_closure=BoundaryClosure.None_,\n", + ")\n", + "\n", + "rmt_total_seconds = perf_counter() - rmt_total_start\n", + "\n", + "comparison_rmt_vertices = np.asarray(\n", + " comparison_rmt_mesh.vertices,\n", + " dtype=np.float64,\n", + ")\n", + "comparison_rmt_cells = np.asarray(\n", + " comparison_rmt_mesh.facets,\n", + " dtype=np.uint32,\n", + ")\n", + "\n", + "print(\"RMT vertices:\", comparison_rmt_vertices.shape)\n", + "print(\"RMT triangles:\", comparison_rmt_cells.shape)\n", + "print(\"RMT total seconds:\", rmt_total_seconds)" + ], + "id": "c20c7765fc73c434", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Compare geometry and topology\n", + "\n", + "The triangle-quality score below is 1.0 for an equilateral triangle and approaches zero for a very thin or degenerate triangle. The 5th percentile is particularly useful because it describes the poor-quality end of each mesh without depending on only one extreme triangle.\n", + "\n", + "An edge used by more than two triangles is non-manifold. Boundary edges are expected here because the folded horizon reaches the sides of the extraction box and boundary closure was disabled. Therefore, a non-zero boundary-edge count does not by itself indicate a broken mesh." + ], + "id": "22b19f3ec6d5d959" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "def calculate_mesh_statistics(\n", + " mesh_vertices: np.ndarray,\n", + " mesh_cells: np.ndarray,\n", + " elapsed_seconds: float,\n", + ") -> tuple[dict[str, float | int], np.ndarray, np.ndarray]:\n", + " \"\"\"Calculate comparable geometric and topological mesh diagnostics.\"\"\"\n", + " triangles = mesh_vertices[mesh_cells]\n", + "\n", + " edge_01 = triangles[:, 1] - triangles[:, 0]\n", + " edge_12 = triangles[:, 2] - triangles[:, 1]\n", + " edge_20 = triangles[:, 0] - triangles[:, 2]\n", + "\n", + " squared_edge_sum = (\n", + " np.sum(edge_01**2, axis=1)\n", + " + np.sum(edge_12**2, axis=1)\n", + " + np.sum(edge_20**2, axis=1)\n", + " )\n", + " twice_area = np.linalg.norm(\n", + " np.cross(edge_01, -edge_20),\n", + " axis=1,\n", + " )\n", + " triangle_quality = np.divide(\n", + " 2.0 * np.sqrt(3.0) * twice_area,\n", + " squared_edge_sum,\n", + " out=np.zeros_like(twice_area),\n", + " where=squared_edge_sum > 0.0,\n", + " )\n", + "\n", + " mesh_edges = np.vstack(\n", + " [\n", + " mesh_cells[:, [0, 1]],\n", + " mesh_cells[:, [1, 2]],\n", + " mesh_cells[:, [2, 0]],\n", + " ]\n", + " )\n", + " mesh_edges.sort(axis=1)\n", + " _unique_edges, edge_use_counts = np.unique(\n", + " mesh_edges,\n", + " axis=0,\n", + " return_counts=True,\n", + " )\n", + "\n", + " fitted_mesh_values = np.asarray(\n", + " rbf_model.evaluate(mesh_vertices),\n", + " dtype=np.float64,\n", + " ).reshape(-1)\n", + " mesh_true_vertical_errors = true_geological_field(mesh_vertices)\n", + "\n", + " statistics = {\n", + " \"vertices\": int(mesh_vertices.shape[0]),\n", + " \"triangles\": int(mesh_cells.shape[0]),\n", + " \"elapsed_seconds\": float(elapsed_seconds),\n", + " \"quality_mean\": float(np.mean(triangle_quality)),\n", + " \"quality_5th_percentile\": float(\n", + " np.percentile(triangle_quality, 5)\n", + " ),\n", + " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", + " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", + " \"mean_absolute_isovalue_error\": float(\n", + " np.mean(np.abs(fitted_mesh_values - folded_isovalue))\n", + " ),\n", + " \"mean_absolute_vertical_error\": float(\n", + " np.mean(np.abs(mesh_true_vertical_errors))\n", + " ),\n", + " \"maximum_absolute_vertical_error\": float(\n", + " np.max(np.abs(mesh_true_vertical_errors))\n", + " ),\n", + " }\n", + "\n", + " return statistics, fitted_mesh_values, mesh_true_vertical_errors\n", + "\n", + "\n", + "mc_statistics, mc_fitted_values, mc_true_vertical_errors = (\n", + " calculate_mesh_statistics(\n", + " mc_vertices,\n", + " mc_cells,\n", + " marching_cubes_total_seconds,\n", + " )\n", + ")\n", + "rmt_statistics, rmt_fitted_values, rmt_true_vertical_errors = (\n", + " calculate_mesh_statistics(\n", + " comparison_rmt_vertices,\n", + " comparison_rmt_cells,\n", + " rmt_total_seconds,\n", + " )\n", + ")\n", + "\n", + "comparison_rows = [\n", + " (\"Vertices\", \"vertices\", \".0f\"),\n", + " (\"Triangles\", \"triangles\", \".0f\"),\n", + " (\"Total time (seconds)\", \"elapsed_seconds\", \".4f\"),\n", + " (\"Mean triangle quality\", \"quality_mean\", \".4f\"),\n", + " (\"5th-percentile quality\", \"quality_5th_percentile\", \".4f\"),\n", + " (\"Boundary edges\", \"boundary_edges\", \".0f\"),\n", + " (\"Non-manifold edges\", \"non_manifold_edges\", \".0f\"),\n", + " (\"Mean absolute isovalue error\", \"mean_absolute_isovalue_error\", \".4f\"),\n", + " (\"Mean absolute vertical error\", \"mean_absolute_vertical_error\", \".4f\"),\n", + " (\"Maximum absolute vertical error\", \"maximum_absolute_vertical_error\", \".4f\"),\n", + "]\n", + "\n", + "print(f\"{'Metric':38s} {'Marching cubes':>18s} {'Ferreus RMT':>18s}\")\n", + "print(\"-\" * 76)\n", + "for label, key, value_format in comparison_rows:\n", + " mc_text = format(mc_statistics[key], value_format)\n", + " rmt_text = format(rmt_statistics[key], value_format)\n", + " print(f\"{label:38s} {mc_text:>18s} {rmt_text:>18s}\")" + ], + "id": "412f81de43194b5a", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### View the two extraction results\n", + "\n", + "The surfaces should have the same large-scale folded geometry because they represent the same zero-isovalue of the same RBF. Differences to look for are triangulation density, skinny triangles, boundary shape, and small positional deviations introduced by RMT's regularisation and vertex clustering, though those things are difficult to discern when there are so many triangles." + ], + "id": "d85ca22ca792b883" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "figure = plt.figure(figsize=(14, 6))\n", + "\n", + "comparison_meshes = [\n", + " (\"Existing method: marching cubes\", mc_vertices, mc_cells, \"darkorange\"),\n", + " (\"Ferreus RMT\", comparison_rmt_vertices, comparison_rmt_cells, \"steelblue\"),\n", + "]\n", + "\n", + "for plot_index, (title, plot_vertices, plot_cells, colour) in enumerate(\n", + " comparison_meshes,\n", + " start=1,\n", + "):\n", + " axis = figure.add_subplot(1, 2, plot_index, projection=\"3d\")\n", + " axis.plot_trisurf(\n", + " plot_vertices[:, 0],\n", + " plot_vertices[:, 1],\n", + " plot_vertices[:, 2],\n", + " triangles=plot_cells,\n", + " color=colour,\n", + " linewidth=0.05,\n", + " alpha=0.85,\n", + " )\n", + " axis.scatter(\n", + " folded_seed_points[:, 0],\n", + " folded_seed_points[:, 1],\n", + " folded_seed_points[:, 2],\n", + " color=\"black\",\n", + " s=8,\n", + " label=\"RMT seeds\",\n", + " )\n", + " axis.set_xlim(MODEL_EXTENTS[[0, 3]])\n", + " axis.set_ylim(MODEL_EXTENTS[[1, 4]])\n", + " axis.set_zlim(MODEL_EXTENTS[[2, 5]])\n", + " axis.set_box_aspect(MODEL_EXTENTS[3:] - MODEL_EXTENTS[:3])\n", + " axis.set_xlabel(\"X\")\n", + " axis.set_ylabel(\"Y\")\n", + " axis.set_zlabel(\"Z\")\n", + " axis.view_init(elev=25, azim=-55)\n", + " axis.set_title(title)\n", + "\n", + "figure.tight_layout()\n", + "plt.show()" + ], + "id": "c46d0141b8471c12", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### Export both comparison surfaces to geoh5\n", + "\n", + "Both meshes are placed in the same group so they can be toggled, coloured, and compared in Geoscience ANALYST. Each receives the RBF value and known vertical error at its vertices." + ], + "id": "2852865e68b17840" + }, + { + "metadata": {}, + "cell_type": "code", + "source": [ + "comparison_output_path = output_directory / (\n", + " f\"marching_cubes_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", + ")\n", + "\n", + "with Workspace.create(comparison_output_path) as workspace:\n", + " comparison_group = ContainerGroup.create(\n", + " workspace,\n", + " name=\"Marching cubes versus Ferreus RMT\",\n", + " )\n", + "\n", + " mc_surface = Surface.create(\n", + " workspace,\n", + " name=\"Existing method - marching cubes\",\n", + " vertices=mc_vertices,\n", + " cells=mc_cells,\n", + " parent=comparison_group,\n", + " )\n", + " mc_surface.add_data(\n", + " {\n", + " \"RBF value on surface\": {\"values\": mc_fitted_values},\n", + " \"True vertical error\": {\"values\": mc_true_vertical_errors},\n", + " }\n", + " )\n", + "\n", + " rmt_surface = Surface.create(\n", + " workspace,\n", + " name=\"Ferreus RMT\",\n", + " vertices=comparison_rmt_vertices,\n", + " cells=comparison_rmt_cells,\n", + " parent=comparison_group,\n", + " )\n", + " rmt_surface.add_data(\n", + " {\n", + " \"RBF value on surface\": {\"values\": rmt_fitted_values},\n", + " \"True vertical error\": {\"values\": rmt_true_vertical_errors},\n", + " }\n", + " )\n", + "\n", + "print(\"Created:\", comparison_output_path)" + ], + "id": "ccfd27e9a9805429", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "### What this comparison does and does not establish\n", + "\n", + "This establishes whether RMT is a credible alternative mesh extractor for a field that already exists. It can compare mesh density, triangle quality, topology, approximate accuracy, and timing under controlled conditions.\n", + "\n", + "It does not establish support for geological orientation or exclusion constraints. Those affect how the scalar field is fitted from observations, which happens before either marching cubes or RMT is called." + ], + "id": "e66498617a269a8f" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb new file mode 100644 index 0000000..f42388c --- /dev/null +++ b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb @@ -0,0 +1,707 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8cbb863ab5d141a7", + "metadata": {}, + "source": [ + "# Existing iso-surfaces versus Ferreus RMT on real geochemistry data\n", + "\n", + "This notebook compares the existing surface_apps.iso_surfaces marching-cubes extraction with ferreus_rmt using the geochem_pts dataset.\n", + "\n", + "The comparison deliberately holds the interpolated scalar field constant:\n", + "\n", + "1. the existing application interpolates Cu_ppm observations onto a regular grid;\n", + "2. marching cubes extracts a surface from that grid; and\n", + "3. Ferreus RMT extracts the same isovalue from a trilinear function over that same grid.\n", + "\n", + "Consequently, differences measured here are primarily differences between the mesh extractors." + ] + }, + { + "cell_type": "code", + "id": "d0d17fd27a2945db", + "metadata": {}, + "source": [ + "from datetime import datetime\n", + "from pathlib import Path\n", + "from time import perf_counter\n", + "\n", + "from ferreus_rmt import BoundaryClosure, ClusterMethod, build_isosurface\n", + "from geoh5py.groups import ContainerGroup\n", + "from geoh5py.objects import Points, Surface\n", + "from geoh5py.workspace import Workspace\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.interpolate import RegularGridInterpolator\n", + "from scipy.sparse import coo_matrix\n", + "from scipy.sparse.csgraph import connected_components\n", + "from scipy.spatial import cKDTree\n", + "\n", + "from surface_apps.iso_surfaces.utils import (\n", + " entity_to_grid,\n", + " extract_iso_surfaces,\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "d2863f57f31b4888", + "metadata": {}, + "source": [ + "## Configuration\n", + "\n", + "Cu_ppm = 200 was selected as a useful first anomaly boundary: it is above most observations but produces enough geometry for a meaningful comparison. These are experiment parameters, not universally correct geological thresholds.\n", + "\n", + "The 50 m grid resolution and 500 m maximum interpolation distance match the existing application's UI defaults. A grid point farther than the maximum distance from data remains unsupported (NaN)." + ] + }, + { + "cell_type": "code", + "id": "77ded4c6f1424af2", + "metadata": {}, + "source": [ + "input_candidates = [\n", + " Path.cwd() / \"geochem_pts\" / \"geochem_pts_example.geoh5\",\n", + " Path.cwd()\n", + " / \"surface_apps\"\n", + " / \"ferreus\"\n", + " / \"geochem_pts\"\n", + " / \"geochem_pts_example.geoh5\",\n", + "]\n", + "\n", + "INPUT_PATH = next(\n", + " (candidate for candidate in input_candidates if candidate.exists()),\n", + " None,\n", + ")\n", + "if INPUT_PATH is None:\n", + " raise FileNotFoundError(\n", + " \"Could not find geochem_pts/geochem_pts_example.geoh5 \"\n", + " \"from the current working directory.\"\n", + " )\n", + "\n", + "OBJECT_NAME = \"geochem_pts\"\n", + "DATA_NAME = \"Cu_ppm\"\n", + "CONTOUR_LEVEL = 200.0\n", + "RESOLUTION = 50.0\n", + "MAX_DISTANCE = 500.0\n", + "RMT_RESOLUTION = RESOLUTION\n", + "RMT_CLUSTER_METHOD = ClusterMethod.CurvatureWeighted\n", + "\n", + "OUTPUT_DIRECTORY = INPUT_PATH.parent / \"ferreus_output\"\n", + "OUTPUT_DIRECTORY.mkdir(exist_ok=True)\n", + "\n", + "print(\"Input:\", INPUT_PATH)\n", + "print(\"Field and isovalue:\", DATA_NAME, CONTOUR_LEVEL)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "be709f46aca747c0", + "metadata": {}, + "source": [ + "## Load the observations and create the shared field\n", + "\n", + "entity_to_grid is the existing application's preprocessing function. For a points object, it uses distance-weighted averaging to estimate the geochemical value at regular-grid locations.\n", + "\n", + "This conversion is shared by both methods below." + ] + }, + { + "cell_type": "code", + "id": "90a877f44c8b44fb", + "metadata": {}, + "source": [ + "grid_start = perf_counter()\n", + "\n", + "with Workspace(INPUT_PATH) as input_workspace:\n", + " source_object = input_workspace.get_entity(OBJECT_NAME)[0]\n", + " source_data = input_workspace.get_entity(DATA_NAME)[0]\n", + "\n", + " if source_object is None or source_data is None:\n", + " raise ValueError(\n", + " f\"Could not find {OBJECT_NAME!r} and {DATA_NAME!r}.\"\n", + " )\n", + "\n", + " source_locations = np.asarray(\n", + " source_object.locations,\n", + " dtype=np.float64,\n", + " ).copy()\n", + " source_values = np.asarray(\n", + " source_data.values,\n", + " dtype=np.float64,\n", + " ).copy()\n", + "\n", + " comparison_grid, comparison_grid_values = entity_to_grid(\n", + " source_object,\n", + " source_data,\n", + " resolution=RESOLUTION,\n", + " max_distance=MAX_DISTANCE,\n", + " horizon=None,\n", + " )\n", + "\n", + "grid_seconds = perf_counter() - grid_start\n", + "\n", + "finite_source_values = source_values[np.isfinite(source_values)]\n", + "finite_grid = np.isfinite(comparison_grid_values)\n", + "\n", + "print(\"Observation count:\", source_locations.shape[0])\n", + "print(\"XYZ span:\", np.ptp(source_locations, axis=0))\n", + "print(\"Finite Cu observations:\", finite_source_values.size)\n", + "print(\"Cu observation percentiles:\", np.percentile(\n", + " finite_source_values,\n", + " [10, 25, 50, 75, 90, 95, 99],\n", + "))\n", + "print(\"Grid shape:\", comparison_grid_values.shape)\n", + "print(\n", + " \"Supported grid nodes:\",\n", + " int(np.count_nonzero(finite_grid)),\n", + " \"/\",\n", + " comparison_grid_values.size,\n", + ")\n", + "print(\"Shared gridding time (seconds):\", grid_seconds)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "0a47f6b27d024ff5", + "metadata": {}, + "source": [ + "## Existing method: marching cubes\n", + "\n", + "This calls the existing application's extraction helper directly. Its result is the baseline users currently receive after the gridding step." + ] + }, + { + "cell_type": "code", + "id": "a953d300954d4afa", + "metadata": {}, + "source": [ + "marching_cubes_start = perf_counter()\n", + "\n", + "existing_surfaces = extract_iso_surfaces(\n", + " source_object,\n", + " comparison_grid,\n", + " [CONTOUR_LEVEL],\n", + " comparison_grid_values,\n", + ")\n", + "\n", + "marching_cubes_seconds = perf_counter() - marching_cubes_start\n", + "\n", + "mc_vertices = np.asarray(existing_surfaces[0][0], dtype=np.float64)\n", + "mc_cells = np.asarray(existing_surfaces[0][1], dtype=np.uint32)\n", + "\n", + "print(\"Marching-cubes vertices:\", mc_vertices.shape)\n", + "print(\"Marching-cubes triangles:\", mc_cells.shape)\n", + "print(\"Marching-cubes extraction time (seconds):\", marching_cubes_seconds)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "db6dad2bde0944d1", + "metadata": {}, + "source": [ + "## Ferreus RMT on the same regular-grid field\n", + "\n", + "RMT expects a callable scalar function instead of a dense value array. RegularGridInterpolator turns the existing grid into that function using trilinear interpolation. Unsupported regions remain NaN.\n", + "\n", + "As a surface-following method, RMT also requires seed points near each surface component. The helper below independently finds regular-grid cells whose eight finite corner values straddle the requested isovalue, then uses their centres as seeds. It does not use the marching-cubes result. This seed search scans the dense grid." + ] + }, + { + "cell_type": "code", + "id": "da81ec4af37d4cd2", + "metadata": {}, + "source": [ + "grid_interpolator = RegularGridInterpolator(\n", + " tuple(comparison_grid),\n", + " comparison_grid_values,\n", + " method=\"linear\",\n", + " bounds_error=False,\n", + " fill_value=np.nan,\n", + ")\n", + "\n", + "\n", + "def shared_grid_field(targets: np.ndarray) -> np.ndarray:\n", + " \"\"\"Evaluate the same trilinear scalar field represented by the grid.\"\"\"\n", + " return np.asarray(grid_interpolator(targets), dtype=np.float64)\n", + "\n", + "\n", + "def find_crossing_cell_seeds(\n", + " grid: list[np.ndarray],\n", + " grid_values: np.ndarray,\n", + " isovalue: float,\n", + ") -> np.ndarray:\n", + " \"\"\"Return centres of fully supported grid cells crossing an isovalue.\"\"\"\n", + " corner_values = np.stack(\n", + " [\n", + " grid_values[x_slice, y_slice, z_slice]\n", + " for x_slice in (slice(None, -1), slice(1, None))\n", + " for y_slice in (slice(None, -1), slice(1, None))\n", + " for z_slice in (slice(None, -1), slice(1, None))\n", + " ],\n", + " axis=-1,\n", + " )\n", + "\n", + " fully_supported = np.all(np.isfinite(corner_values), axis=-1)\n", + " crosses_isovalue = (\n", + " fully_supported\n", + " & (np.min(corner_values, axis=-1) <= isovalue)\n", + " & (np.max(corner_values, axis=-1) >= isovalue)\n", + " )\n", + "\n", + " crossing_indices = np.argwhere(crosses_isovalue)\n", + " if crossing_indices.size == 0:\n", + " raise ValueError(f\"No fully supported grid cells cross {isovalue}.\")\n", + "\n", + " return np.column_stack(\n", + " [\n", + " 0.5\n", + " * (\n", + " grid[axis][crossing_indices[:, axis]]\n", + " + grid[axis][crossing_indices[:, axis] + 1]\n", + " )\n", + " for axis in range(3)\n", + " ]\n", + " ).astype(np.float64)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "id": "aa3c77a2469f4b56", + "metadata": {}, + "source": [ + "rmt_total_start = perf_counter()\n", + "\n", + "rmt_seed_points = find_crossing_cell_seeds(\n", + " comparison_grid,\n", + " comparison_grid_values,\n", + " CONTOUR_LEVEL,\n", + ")\n", + "seed_detection_seconds = perf_counter() - rmt_total_start\n", + "\n", + "rmt_extents = np.array(\n", + " [\n", + " comparison_grid[0][0],\n", + " comparison_grid[1][0],\n", + " comparison_grid[2][0],\n", + " comparison_grid[0][-1],\n", + " comparison_grid[1][-1],\n", + " comparison_grid[2][-1],\n", + " ],\n", + " dtype=np.float64,\n", + ")\n", + "\n", + "rmt_extraction_start = perf_counter()\n", + "rmt_mesh = build_isosurface(\n", + " seed_points=rmt_seed_points,\n", + " extents=rmt_extents,\n", + " resolution=RMT_RESOLUTION,\n", + " isovalue=CONTOUR_LEVEL,\n", + " surface_fn=shared_grid_field,\n", + " cluster_method=RMT_CLUSTER_METHOD,\n", + " boundary_closure=BoundaryClosure.None_,\n", + ")\n", + "rmt_extraction_seconds = perf_counter() - rmt_extraction_start\n", + "rmt_total_seconds = perf_counter() - rmt_total_start\n", + "\n", + "rmt_vertices = np.asarray(rmt_mesh.vertices, dtype=np.float64)\n", + "rmt_cells = np.asarray(rmt_mesh.facets, dtype=np.uint32)\n", + "\n", + "print(\"Independent RMT seed points:\", rmt_seed_points.shape)\n", + "print(\"RMT vertices:\", rmt_vertices.shape)\n", + "print(\"RMT triangles:\", rmt_cells.shape)\n", + "print(\"Seed detection time (seconds):\", seed_detection_seconds)\n", + "print(\"RMT extraction time (seconds):\", rmt_extraction_seconds)\n", + "print(\"RMT seed plus extraction time (seconds):\", rmt_total_seconds)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "c55a39891ad44956", + "metadata": {}, + "source": [ + "## Compare mesh quality, topology, and field agreement\n", + "\n", + "Triangle quality is 1.0 for an equilateral triangle and approaches zero for increasingly skinny triangles. An edge shared by more than two triangles is non-manifold. Boundary edges can occur where a surface meets an unsupported region or the extraction bounds.\n", + "\n", + "There is no known true copper-anomaly surface in this dataset, so this notebook cannot report geological accuracy. Instead it measures how closely each mesh follows the shared interpolated 200 ppm field and how closely the two meshes agree spatially." + ] + }, + { + "cell_type": "code", + "id": "e14ca8c55c72457b", + "metadata": {}, + "source": [ + "def calculate_mesh_statistics(\n", + " mesh_vertices: np.ndarray,\n", + " mesh_cells: np.ndarray,\n", + " extraction_seconds: float,\n", + ") -> tuple[dict[str, float | int], np.ndarray]:\n", + " \"\"\"Calculate geometric, topological, and scalar-field diagnostics.\"\"\"\n", + " triangles = mesh_vertices[mesh_cells]\n", + " edge_01 = triangles[:, 1] - triangles[:, 0]\n", + " edge_12 = triangles[:, 2] - triangles[:, 1]\n", + " edge_20 = triangles[:, 0] - triangles[:, 2]\n", + "\n", + " twice_area = np.linalg.norm(np.cross(edge_01, -edge_20), axis=1)\n", + " squared_edge_sum = (\n", + " np.sum(edge_01**2, axis=1)\n", + " + np.sum(edge_12**2, axis=1)\n", + " + np.sum(edge_20**2, axis=1)\n", + " )\n", + " triangle_quality = np.divide(\n", + " 2.0 * np.sqrt(3.0) * twice_area,\n", + " squared_edge_sum,\n", + " out=np.zeros_like(twice_area),\n", + " where=squared_edge_sum > 0.0,\n", + " )\n", + "\n", + " mesh_edges = np.vstack(\n", + " [\n", + " mesh_cells[:, [0, 1]],\n", + " mesh_cells[:, [1, 2]],\n", + " mesh_cells[:, [2, 0]],\n", + " ]\n", + " )\n", + " mesh_edges.sort(axis=1)\n", + " unique_edges, edge_use_counts = np.unique(\n", + " mesh_edges,\n", + " axis=0,\n", + " return_counts=True,\n", + " )\n", + "\n", + " graph = coo_matrix(\n", + " (\n", + " np.ones(unique_edges.shape[0] * 2),\n", + " (\n", + " np.r_[unique_edges[:, 0], unique_edges[:, 1]],\n", + " np.r_[unique_edges[:, 1], unique_edges[:, 0]],\n", + " ),\n", + " ),\n", + " shape=(mesh_vertices.shape[0], mesh_vertices.shape[0]),\n", + " )\n", + " component_count = connected_components(\n", + " graph,\n", + " directed=False,\n", + " return_labels=False,\n", + " )\n", + "\n", + " field_values = shared_grid_field(mesh_vertices)\n", + " finite_field = np.isfinite(field_values)\n", + " isovalue_errors = np.abs(field_values[finite_field] - CONTOUR_LEVEL)\n", + "\n", + " statistics = {\n", + " \"vertices\": int(mesh_vertices.shape[0]),\n", + " \"triangles\": int(mesh_cells.shape[0]),\n", + " \"extraction_seconds\": float(extraction_seconds),\n", + " \"surface_area\": float(np.sum(0.5 * twice_area)),\n", + " \"quality_mean\": float(np.mean(triangle_quality)),\n", + " \"quality_5th_percentile\": float(\n", + " np.percentile(triangle_quality, 5)\n", + " ),\n", + " \"components\": int(component_count),\n", + " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", + " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", + " \"finite_field_fraction\": float(np.mean(finite_field)),\n", + " \"mean_absolute_isovalue_error\": float(np.mean(isovalue_errors)),\n", + " \"maximum_absolute_isovalue_error\": float(np.max(isovalue_errors)),\n", + " }\n", + "\n", + " return statistics, field_values\n", + "\n", + "\n", + "mc_statistics, mc_field_values = calculate_mesh_statistics(\n", + " mc_vertices,\n", + " mc_cells,\n", + " marching_cubes_seconds,\n", + ")\n", + "rmt_statistics, rmt_field_values = calculate_mesh_statistics(\n", + " rmt_vertices,\n", + " rmt_cells,\n", + " rmt_extraction_seconds,\n", + ")\n", + "\n", + "comparison_rows = [\n", + " (\"Vertices\", \"vertices\", \".0f\"),\n", + " (\"Triangles\", \"triangles\", \".0f\"),\n", + " (\"Extraction seconds\", \"extraction_seconds\", \".4f\"),\n", + " (\"Surface area\", \"surface_area\", \".1f\"),\n", + " (\"Mean triangle quality\", \"quality_mean\", \".4f\"),\n", + " (\"5th-percentile triangle quality\", \"quality_5th_percentile\", \".4f\"),\n", + " (\"Connected components\", \"components\", \".0f\"),\n", + " (\"Boundary edges\", \"boundary_edges\", \".0f\"),\n", + " (\"Non-manifold edges\", \"non_manifold_edges\", \".0f\"),\n", + " (\"Finite field-value fraction\", \"finite_field_fraction\", \".4f\"),\n", + " (\"Mean absolute isovalue error\", \"mean_absolute_isovalue_error\", \".4f\"),\n", + " (\"Maximum absolute isovalue error\", \"maximum_absolute_isovalue_error\", \".4f\"),\n", + "]\n", + "\n", + "print(f\"{'Metric':40s} {'Marching cubes':>18s} {'Ferreus RMT':>18s}\")\n", + "print(\"-\" * 78)\n", + "for label, key, value_format in comparison_rows:\n", + " mc_text = format(mc_statistics[key], value_format)\n", + " rmt_text = format(rmt_statistics[key], value_format)\n", + " print(f\"{label:40s} {mc_text:>18s} {rmt_text:>18s}\")\n", + "\n", + "\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "450f5356b5304251", + "metadata": {}, + "source": [ + "### Direct positional agreement\n", + "\n", + "Nearest-vertex distances provide an approximate comparison between the meshes. They are not exact point-to-triangle distances, so interpret them relative to the configured grid resolution." + ] + }, + { + "cell_type": "code", + "id": "a85ce67ad448473a", + "metadata": {}, + "source": [ + "mc_to_rmt_distances = cKDTree(rmt_vertices).query(mc_vertices)[0]\n", + "rmt_to_mc_distances = cKDTree(mc_vertices).query(rmt_vertices)[0]\n", + "symmetric_distances = np.r_[mc_to_rmt_distances, rmt_to_mc_distances]\n", + "\n", + "print(\"Marching cubes to RMT mean distance:\", mc_to_rmt_distances.mean())\n", + "print(\"RMT to marching cubes mean distance:\", rmt_to_mc_distances.mean())\n", + "print(\n", + " \"Symmetric distance percentiles (50%, 90%, 95%, 100%):\",\n", + " np.percentile(symmetric_distances, [50, 90, 95, 100]),\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "fa7ef208f1dd4dab", + "metadata": {}, + "source": [ + "### Reading the result\n", + "\n", + "Compare triangle quality, component counts, field error, and positional separation together. A method can produce attractive triangles while missing small components, or reproduce the isovalue accurately while retaining questionable interpolation artefacts. Because this dataset has no known true 200 ppm surface, neither output can be labelled geologically correct from these metrics alone." + ] + }, + { + "cell_type": "markdown", + "id": "2f3b645319f045de", + "metadata": {}, + "source": [ + "## Visual comparison\n", + "\n", + "The source points shown are finite observations at or above 200 ppm. They are not expected to lie exactly on the interpolated surface: the surface encloses/intersects regions where the continuous estimated field equals 200 ppm.\n", + "\n", + "Matplotlib's 3D renderer seems to have issues even when I tried to limit the complexity of what it's displaying, so I couldn't get this to work and it just keeps working for minutes." + ] + }, + { + "cell_type": "code", + "id": "21305e00d59c4e30", + "metadata": {}, + "source": [ + "'''MAX_PREVIEW_TRIANGLES = 2_500\n", + "\n", + "\n", + "def select_preview_cells(\n", + " cells: np.ndarray,\n", + " maximum_triangles: int = MAX_PREVIEW_TRIANGLES,\n", + ") -> np.ndarray:\n", + " \"\"\"Select evenly spaced triangles for responsive notebook rendering.\"\"\"\n", + " if cells.shape[0] <= maximum_triangles:\n", + " return cells\n", + "\n", + " preview_indices = np.linspace(\n", + " 0,\n", + " cells.shape[0] - 1,\n", + " maximum_triangles,\n", + " dtype=int,\n", + " )\n", + " return cells[preview_indices]\n", + "\n", + "\n", + "high_copper = np.isfinite(source_values) & (source_values >= CONTOUR_LEVEL)\n", + "coordinate_span = np.ptp(source_locations, axis=0)\n", + "\n", + "figure = plt.figure(figsize=(15, 7))\n", + "plot_meshes = [\n", + " (\"Existing marching cubes\", mc_vertices, mc_cells, \"darkorange\"),\n", + " (\"Ferreus RMT\", rmt_vertices, rmt_cells, \"steelblue\"),\n", + "]\n", + "\n", + "for plot_index, (title, vertices, full_cells, colour) in enumerate(\n", + " plot_meshes,\n", + " start=1,\n", + "):\n", + " preview_cells = select_preview_cells(full_cells)\n", + " print(\n", + " f\"Rendering {title}: {preview_cells.shape[0]:,} of \"\n", + " f\"{full_cells.shape[0]:,} triangles\"\n", + " )\n", + "\n", + " axis = figure.add_subplot(1, 2, plot_index, projection=\"3d\")\n", + " axis.plot_trisurf(\n", + " vertices[:, 0],\n", + " vertices[:, 1],\n", + " vertices[:, 2],\n", + " triangles=preview_cells,\n", + " color=colour,\n", + " linewidth=0.0,\n", + " antialiased=False,\n", + " shade=False,\n", + " alpha=1.0,\n", + " )\n", + " axis.scatter(\n", + " source_locations[high_copper, 0],\n", + " source_locations[high_copper, 1],\n", + " source_locations[high_copper, 2],\n", + " color=\"black\",\n", + " s=5,\n", + " alpha=0.6,\n", + " label=f\"Observed {DATA_NAME} >= {CONTOUR_LEVEL:g}\",\n", + " )\n", + " axis.set_box_aspect(coordinate_span)\n", + " axis.set_xlabel(\"X\")\n", + " axis.set_ylabel(\"Y\")\n", + " axis.set_zlabel(\"Z\")\n", + " axis.view_init(elev=22, azim=-55)\n", + " axis.set_title(title)\n", + " axis.legend(loc=\"upper right\")\n", + "\n", + "figure.suptitle(f\"{DATA_NAME} = {CONTOUR_LEVEL:g} ppm\")\n", + "figure.subplots_adjust(left=0.02, right=0.98, bottom=0.02, top=0.90)\n", + "plt.show()'''" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "8b950ca8bc474c13", + "metadata": {}, + "source": [ + "## Export both methods to geoh5\n", + "\n", + "The output includes a copy of the geochemistry points and both surfaces. Each mesh has its interpolated field value and nearest-vertex distance to the other method, making it possible to inspect differences in Geoscience ANALYST." + ] + }, + { + "cell_type": "code", + "id": "62451e97bef04a06", + "metadata": {}, + "source": [ + "# Recreate the output directory if previous results were deleted after\n", + "# the configuration cell was run.\n", + "OUTPUT_DIRECTORY.mkdir(parents=True, exist_ok=True)\n", + "\n", + "output_path = OUTPUT_DIRECTORY / (\n", + " f\"real_geochem_mc_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", + ")\n", + "\n", + "with Workspace.create(output_path) as output_workspace:\n", + " result_group = ContainerGroup.create(\n", + " output_workspace,\n", + " name=f\"{DATA_NAME} {CONTOUR_LEVEL:g} - MC versus RMT\",\n", + " )\n", + "\n", + " output_points = Points.create(\n", + " output_workspace,\n", + " name=OBJECT_NAME,\n", + " vertices=source_locations,\n", + " parent=result_group,\n", + " )\n", + " output_points.add_data({DATA_NAME: {\"values\": source_values}})\n", + "\n", + " output_mc_surface = Surface.create(\n", + " output_workspace,\n", + " name=f\"Existing marching cubes - {CONTOUR_LEVEL:g} ppm\",\n", + " vertices=mc_vertices,\n", + " cells=mc_cells,\n", + " parent=result_group,\n", + " )\n", + " output_mc_surface.add_data(\n", + " {\n", + " \"Interpolated Cu ppm\": {\"values\": mc_field_values},\n", + " \"Nearest RMT vertex distance\": {\n", + " \"values\": mc_to_rmt_distances\n", + " },\n", + " }\n", + " )\n", + "\n", + " output_rmt_surface = Surface.create(\n", + " output_workspace,\n", + " name=f\"Ferreus RMT - {CONTOUR_LEVEL:g} ppm\",\n", + " vertices=rmt_vertices,\n", + " cells=rmt_cells,\n", + " parent=result_group,\n", + " )\n", + " output_rmt_surface.add_data(\n", + " {\n", + " \"Interpolated Cu ppm\": {\"values\": rmt_field_values},\n", + " \"Nearest marching-cubes vertex distance\": {\n", + " \"values\": rmt_to_mc_distances\n", + " },\n", + " }\n", + " )\n", + "\n", + "print(\"Created:\", output_path)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "e642a55b31c948eb", + "metadata": {}, + "source": [ + "## Interpretation limits and next tests\n", + "\n", + "A successful run establishes that RMT can extract real, disconnected geochemical anomaly surfaces from the same field used by the current application, including in the presence of unsupported regions.\n", + "\n", + "It does not show that either surface is geological truth. Both inherit assumptions and artefacts from the shared weighted-average grid. It also does not test the RBF package or constraint fitting.\n", + "\n", + "Useful follow-up experiments are:\n", + "\n", + "- repeat at several contour levels and resolutions;\n", + "- compare results at approximately equal triangle count or positional error;\n", + "- test RMT boundary-closure modes on a genuinely enclosed component;\n", + "- replace grid-cell seed detection with seeds derived from selected observations or another production-safe strategy; and\n", + "- separately compare the existing weighted-average field with a ferreus_rbf field fitted from the raw geochemistry points." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 b/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 new file mode 100644 index 0000000..6d6d881 --- /dev/null +++ b/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b1410455e6510fe01bf6f4f32c380eb04aa9aba4dc8007fa6eecdb6bd1f576af +size 1376256 From 6fe24384076c0c266f70a57f35618999e91d3494 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 12 Aug 2026 15:50:29 +0000 Subject: [PATCH 2/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 504 ++++++++---------- .../ferreus_rmt_real_geochem_comparison.ipynb | 95 ++-- 2 files changed, 280 insertions(+), 319 deletions(-) diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb index 76de64c..684a97f 100644 --- a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -1,14 +1,15 @@ { "cells": [ { + "cell_type": "markdown", + "id": "52652cf59cffb624", "metadata": { - "collapsed": true, "ExecuteTime": { "end_time": "2026-08-10T15:19:45.243382800Z", "start_time": "2026-08-10T15:19:45.191472Z" - } + }, + "collapsed": true }, - "cell_type": "markdown", "source": [ "# Ferreus RMT to geoh5 proof of concept\n", "\n", @@ -21,71 +22,71 @@ " orientation constraints\n", "\n", "POC with notebook creating a surface transferred to geoh5" - ], - "id": "52652cf59cffb624" + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "", - "id": "de308ccc50aea28d" + "id": "de308ccc50aea28d", + "metadata": {}, + "source": "" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "c14d7389d994a354", + "metadata": {}, + "outputs": [], "source": [ "import os\n", "import sys\n", - "from pathlib import Path\n", - "\n", "from datetime import datetime\n", + "from pathlib import Path\n", "\n", - "from geoh5py.groups import ContainerGroup\n", - "from geoh5py.objects import Points, Surface\n", - "from geoh5py.workspace import Workspace\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "\n", - "from ferreus_rmt import (\n", - " BoundaryClosure,\n", - " ClusterMethod,\n", - " build_isosurface,\n", - ")\n", "from ferreus_rbf import RBFInterpolator\n", "from ferreus_rbf.interpolant_config import (\n", " Drift,\n", " InterpolantSettings,\n", " RBFKernelType,\n", ")\n", - "import matplotlib.pyplot as plt" - ], - "id": "c14d7389d994a354", - "outputs": [], - "execution_count": null + "from ferreus_rmt import (\n", + " BoundaryClosure,\n", + " ClusterMethod,\n", + " build_isosurface,\n", + ")\n", + "from geoh5py.groups import ContainerGroup\n", + "from geoh5py.objects import Points, Surface\n", + "from geoh5py.workspace import Workspace" + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### define a known sphere", - "id": "a8c8dc3ebfbae12c" + "id": "a8c8dc3ebfbae12c", + "metadata": {}, + "source": "### define a known sphere" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "9be78b4d9cd5c492", + "metadata": {}, + "outputs": [], "source": [ - "\n", "CENTER = np.array([0.0, 0.0, 0.0], dtype=np.float64)\n", "RADIUS = 50.0\n", "\n", + "\n", "def sphere_values(points: np.ndarray) -> np.ndarray:\n", " \"\"\"negative inside the sphere, zero on it, positive outside.\"\"\"\n", "\n", " offsets = points - CENTER\n", " return np.linalg.norm(offsets, axis=1) - RADIUS\n", "\n", + "\n", "def sphere_values_and_gradients(\n", - " points: np.ndarray,\n", + " points: np.ndarray,\n", ") -> tuple[np.ndarray, np.ndarray]:\n", - "\n", " \"\"\"return values and analytic gradients.\"\"\"\n", " offsets = points - CENTER\n", " distances = np.linalg.norm(offsets, axis=1)\n", @@ -98,24 +99,21 @@ " )\n", "\n", " return distances - RADIUS, gradients" - ], - "id": "9be78b4d9cd5c492", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### extract the surface with rmt", - "id": "bf48bdde9dab3801" + "id": "bf48bdde9dab3801", + "metadata": {}, + "source": "### extract the surface with rmt" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "8d2fa33ef1334bfa", + "metadata": {}, + "outputs": [], "source": [ - "\n", - "\n", - "\n", "seed_points = np.array(\n", " [\n", " [50.0, 0.0, 0.0],\n", @@ -132,17 +130,16 @@ ")\n", "\n", "\n", - "\n", "resolution = 5.0\n", "isovalue = 0.0\n", "\n", "mesh = build_isosurface(\n", - " seed_points, # approx. locations of the surface to seed the wavefront\n", - " extents, # the axis-aligned bounding box limiting extraction\n", - " resolution, # lattice spacing in coordinate units - smaller means more detail\n", - " isovalue, # the scalar value defining the surface\n", - " surface_fn=sphere_values, # surface_fn (sphere_values) is the scalar function that defines the geometry\n", - " gradient_fn=sphere_values_and_gradients, #\n", + " seed_points, # approx. locations of the surface to seed the wavefront\n", + " extents, # the axis-aligned bounding box limiting extraction\n", + " resolution, # lattice spacing in coordinate units - smaller means more detail\n", + " isovalue, # the scalar value defining the surface\n", + " surface_fn=sphere_values, # surface_fn (sphere_values) is the scalar function that defines the geometry\n", + " gradient_fn=sphere_values_and_gradients, #\n", ")\n", "\n", "vertices = np.asarray(mesh.vertices)\n", @@ -150,22 +147,21 @@ "\n", "print(\"Vertices:\", vertices.shape)\n", "print(\"Triangles:\", cells.shape)" - ], - "id": "8d2fa33ef1334bfa", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### validate the result", - "id": "4178f0ab56928714" + "id": "4178f0ab56928714", + "metadata": {}, + "source": "### validate the result" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "c953cf37f21ad3a0", + "metadata": {}, + "outputs": [], "source": [ - "\n", "radii = np.linalg.norm(vertices - CENTER, axis=1)\n", "radius_errors = radii - RADIUS\n", "\n", @@ -179,22 +175,21 @@ "assert np.abs(radius_errors).max() < resolution\n", "\n", "# numbers should be close to the radius, but likely not exact because of the resolution of the lattice." - ], - "id": "c953cf37f21ad3a0", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### show the sphere we've constructed", - "id": "3b1c992eb5ef149d" + "id": "3b1c992eb5ef149d", + "metadata": {}, + "source": "### show the sphere we've constructed" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "8c91eaf424ec37fb", + "metadata": {}, + "outputs": [], "source": [ - "\n", "figure = plt.figure(figsize=(8, 7))\n", "axis = figure.add_subplot(111, projection=\"3d\")\n", "\n", @@ -214,23 +209,21 @@ "axis.set_title(\"Analytic sphere extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ], - "id": "8c91eaf424ec37fb", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Export to geoh5", - "id": "ba702ddc81228ad6" + "id": "ba702ddc81228ad6", + "metadata": {}, + "source": "### Export to geoh5" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "f39977f0e682fcb9", + "metadata": {}, + "outputs": [], "source": [ - "\n", - "\n", "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", "\n", @@ -253,14 +246,12 @@ " )\n", "\n", "print(\"Created:\", output_path)" - ], - "id": "f39977f0e682fcb9", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "1a8d7c5fad014c7f", + "metadata": {}, "source": [ "## Fit a synthetic folded horizon with ferreus_rbf\n", "\n", @@ -272,12 +263,14 @@ "\n", "A known analytic geological field is retained as ground truth, allowing\n", "the interpolation and resulting surface errors to be measured." - ], - "id": "1a8d7c5fad014c7f" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "e8a9c39ec7391034", + "metadata": {}, + "outputs": [], "source": [ "MODEL_EXTENTS = np.array(\n", " [-100.0, -100.0, -60.0, 100.0, 100.0, 60.0],\n", @@ -290,10 +283,7 @@ " y: np.ndarray,\n", ") -> np.ndarray:\n", " \"\"\"Elevation of the synthetic folded horizon.\"\"\"\n", - " return (\n", - " 20.0 * np.sin(x / 35.0)\n", - " + 10.0 * np.cos(y / 40.0)\n", - " )\n", + " return 20.0 * np.sin(x / 35.0) + 10.0 * np.cos(y / 40.0)\n", "\n", "\n", "def true_geological_field(points: np.ndarray) -> np.ndarray:\n", @@ -304,20 +294,20 @@ " )\n", "\n", " return points[:, 2] - horizon_z" - ], - "id": "e8a9c39ec7391034", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Generate scattered scalar observations", - "id": "55396a8e404a0f23" + "id": "55396a8e404a0f23", + "metadata": {}, + "source": "### Generate scattered scalar observations" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "63f6e9858be0976d", + "metadata": {}, + "outputs": [], "source": [ "random_generator = np.random.default_rng(42)\n", "\n", @@ -329,9 +319,7 @@ " size=(number_of_source_points, 3),\n", ").astype(np.float64)\n", "\n", - "source_values = true_geological_field(source_points).astype(\n", - " np.float64\n", - ")\n", + "source_values = true_geological_field(source_points).astype(np.float64)\n", "\n", "print(\"Source points:\", source_points.shape)\n", "print(\"Source values:\", source_values.shape)\n", @@ -341,25 +329,25 @@ " \"to\",\n", " float(source_values.max()),\n", ")" - ], - "id": "63f6e9858be0976d", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Fit the RBF field", - "id": "203c7442a4c55fec" + "id": "203c7442a4c55fec", + "metadata": {}, + "source": "### Fit the RBF field" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "b82e4349d22b0c02", + "metadata": {}, + "outputs": [], "source": [ "rbf_settings = InterpolantSettings(\n", - " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", - " drift=Drift.Linear, # adds a broad linear trend, useful here because the field contains a strong linear z component.\n", - " nugget=0.0, # 0.0 requests an unsmoothed fit. Fine for synthetic observations as they contain no noise.\n", + " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", + " drift=Drift.Linear, # adds a broad linear trend, useful here because the field contains a strong linear z component.\n", + " nugget=0.0, # 0.0 requests an unsmoothed fit. Fine for synthetic observations as they contain no noise.\n", ")\n", "\n", "rbf_model = RBFInterpolator(\n", @@ -367,24 +355,24 @@ " values=source_values,\n", " interpolant_settings=rbf_settings,\n", ")" - ], - "id": "b82e4349d22b0c02", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Check the training residuals", - "id": "6b8e055979af129a" + "id": "6b8e055979af129a", + "metadata": {}, + "source": "### Check the training residuals" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "f1900adf375eec99", + "metadata": {}, + "outputs": [], "source": [ - "source_predictions = np.asarray(\n", - " rbf_model.evaluate_at_source(add_nugget=False)\n", - ").reshape(-1)\n", + "source_predictions = np.asarray(rbf_model.evaluate_at_source(add_nugget=False)).reshape(\n", + " -1\n", + ")\n", "\n", "source_residuals = source_predictions - source_values\n", "\n", @@ -399,20 +387,20 @@ "\n", "# This demonstrates how closely the solved RBF reproduces the values used to fit it.\n", "# For a zero-nugget, synthetic problem, the residuals should be very small, on the order of machine precision." - ], - "id": "f1900adf375eec99", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Test at locations not used for fitting", - "id": "e25a69ffedc71d99" + "id": "e25a69ffedc71d99", + "metadata": {}, + "source": "### Test at locations not used for fitting" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "28d26151ceb1bbf9", + "metadata": {}, + "outputs": [], "source": [ "number_of_validation_points = 200\n", "\n", @@ -422,21 +410,15 @@ " size=(number_of_validation_points, 3),\n", ").astype(np.float64)\n", "\n", - "true_validation_values = true_geological_field(\n", - " validation_points\n", - ")\n", + "true_validation_values = true_geological_field(validation_points)\n", "\n", - "predicted_validation_values = np.asarray(\n", - " rbf_model.evaluate(validation_points)\n", - ").reshape(-1)\n", - "\n", - "validation_residuals = (\n", - " predicted_validation_values - true_validation_values\n", + "predicted_validation_values = np.asarray(rbf_model.evaluate(validation_points)).reshape(\n", + " -1\n", ")\n", "\n", - "validation_rmse = np.sqrt(\n", - " np.mean(validation_residuals**2)\n", - ")\n", + "validation_residuals = predicted_validation_values - true_validation_values\n", + "\n", + "validation_rmse = np.sqrt(np.mean(validation_residuals**2))\n", "\n", "print(\"Validation RMSE:\", float(validation_rmse))\n", "print(\n", @@ -446,20 +428,20 @@ "\n", "# The source residual measures fitting accuracy, i.e. an RMSE of 0.69 means approx 0.69 vertical units\n", "# The validation residual measures interpolation accuracy between observations" - ], - "id": "28d26151ceb1bbf9", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "### Prepare the RBF for repeated RMT evaluations", - "id": "adcfe7b8d0d585b4" + "id": "adcfe7b8d0d585b4", + "metadata": {}, + "source": "### Prepare the RBF for repeated RMT evaluations" }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "3c481cec6ea9a739", + "metadata": {}, + "outputs": [], "source": [ "folded_resolution = 5.0\n", "\n", @@ -480,14 +462,12 @@ ")\n", "\n", "rbf_model.build_evaluator(evaluator_extents)" - ], - "id": "3c481cec6ea9a739", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "e674a1b0a8c42db5", + "metadata": {}, "source": [ "### Choose seed points for the folded horizon\n", "\n", @@ -496,12 +476,14 @@ "These seeds are algorithm inputs, not extra observations used to fit the RBF. In this synthetic example the true horizon is known, so a small grid of exact surface locations is convenient. In a real geological workflow, observed contact points with scalar value zero would be natural seeds.\n", "\n", "Multiple seeds are useful when a model may contain disconnected surface components. A component that has no nearby seed may not be discovered." - ], - "id": "e674a1b0a8c42db5" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "16cbac4b3b5c9df8", + "metadata": {}, + "outputs": [], "source": [ "seed_coordinates = np.linspace(-80.0, 80.0, 5)\n", "seed_x, seed_y = np.meshgrid(\n", @@ -522,17 +504,16 @@ ").astype(np.float64)\n", "\n", "print(\"RMT seed points shape:\", folded_seed_points.shape)\n", - "print(\"Maximum absolute seed field value:\", float(\n", - " np.max(np.abs(true_geological_field(folded_seed_points)))\n", - "))" - ], - "id": "16cbac4b3b5c9df8", - "outputs": [], - "execution_count": null + "print(\n", + " \"Maximum absolute seed field value:\",\n", + " float(np.max(np.abs(true_geological_field(folded_seed_points)))),\n", + ")" + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "f25538fc35ed124b", + "metadata": {}, "source": [ "### Extract the fitted zero-isosurface with RMT\n", "\n", @@ -541,12 +522,14 @@ "folded_resolution is the spacing of RMT's internal sampling lattice in the same coordinate units as X, Y, and Z. 5.0 means samples are approximately five units apart. A smaller value can preserve finer geometry, but requires more field evaluations and usually produces more triangles.\n", "\n", "evaluate_targets and evaluate_targets_with_gradients use the cached evaluator prepared above. The gradients help RMT project its working vertices toward the requested scalar value; they are not geological orientation observations." - ], - "id": "f25538fc35ed124b" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "7a79146bd28e5dd1", + "metadata": {}, + "outputs": [], "source": [ "folded_isovalue = 0.0\n", "\n", @@ -572,14 +555,12 @@ "\n", "print(\"Folded surface vertices:\", folded_vertices.shape)\n", "print(\"Folded surface triangles:\", folded_cells.shape)" - ], - "id": "7a79146bd28e5dd1", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "e203369f6da5b01f", + "metadata": {}, "source": [ "### Validate both parts of the workflow\n", "\n", @@ -589,21 +570,21 @@ "2. **Modelling error:** evaluate the known analytic field at every mesh vertex. Here that value is also the signed vertical distance from the true horizon. This checks how closely the RBF-derived surface matches the ground truth between the observations.\n", "\n", "Keeping these checks separate matters as a mesh can follow an imperfect RBF field very accurately." - ], - "id": "e203369f6da5b01f" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "798b9ee4fbcc148e", + "metadata": {}, + "outputs": [], "source": [ "fitted_values_at_vertices = np.asarray(\n", " rbf_model.evaluate(folded_vertices),\n", " dtype=np.float64,\n", ").reshape(-1)\n", "\n", - "rmt_isovalue_errors = (\n", - " fitted_values_at_vertices - folded_isovalue\n", - ")\n", + "rmt_isovalue_errors = fitted_values_at_vertices - folded_isovalue\n", "\n", "# For this particular synthetic field, this is the signed vertical\n", "# separation between a mesh vertex and the true folded horizon.\n", @@ -627,24 +608,24 @@ "assert folded_cells.ndim == 2 and folded_cells.shape[1] == 3\n", "assert np.all(np.isfinite(folded_vertices))\n", "assert np.all(np.isfinite(fitted_values_at_vertices))" - ], - "id": "798b9ee4fbcc148e", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "be32da3dd2532365", + "metadata": {}, "source": [ "### Inspect the fitted surface and nearby observations\n", "\n", "Only observations reasonably close to the horizon are plotted, which keeps the diagnostic readable. Their colours show the signed scalar values supplied to the RBF: blue points are below the horizon, red points are above it, and values near zero lie close to the interpreted contact." - ], - "id": "be32da3dd2532365" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "40448b38d7a150c4", + "metadata": {}, + "outputs": [], "source": [ "near_horizon = np.abs(source_values) <= 15.0\n", "\n", @@ -690,14 +671,12 @@ "axis.set_title(\"RBF folded horizon extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ], - "id": "40448b38d7a150c4", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "9eab773cbf3c33cd", + "metadata": {}, "source": [ "### Export the folded surface and its observations to geoh5\n", "\n", @@ -709,12 +688,14 @@ "- the signed vertical error relative to the known synthetic horizon, useful only because this POC has ground truth.\n", "\n", "In Geoscience ANALYST, colour the observation points by **Implicit scalar value** and the surface by **True vertical error** to inspect the result." - ], - "id": "9eab773cbf3c33cd" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "156ec95b0d1e1877", + "metadata": {}, + "outputs": [], "source": [ "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", @@ -762,14 +743,12 @@ " )\n", "\n", "print(\"Created:\", folded_output_path)" - ], - "id": "156ec95b0d1e1877", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "f74bf0dad647895f", + "metadata": {}, "source": [ "## What this baseline proves\n", "\n", @@ -787,12 +766,12 @@ "2. approximate an exclusion constraint with a signed pseudo-value on the required side of the surface.\n", "\n", "Those are approximations that should be measured against this baseline." - ], - "id": "f74bf0dad647895f" + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "23da01a0b8b33419", + "metadata": {}, "source": [ "## Compare the existing marching-cubes method with Ferreus RMT\n", "\n", @@ -801,12 +780,14 @@ "The existing application first samples data onto a complete regular grid and then calls scikit-image marching cubes. The code below performs that same dense-grid extraction. RMT instead evaluates the field while following the surface outward from seed points.\n", "\n", "This comparison does not yet reproduce every part of the existing application's preprocessing, such as weighted averaging, maximum distance, or topographic masking. A later end-to-end test can feed both implementations the same real geoh5 object." - ], - "id": "23da01a0b8b33419" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "9117b165c6dfa0e7", + "metadata": {}, + "outputs": [], "source": [ "from time import perf_counter\n", "\n", @@ -868,12 +849,8 @@ ").astype(np.float64)\n", "mc_cells = np.asarray(mc_cells, dtype=np.uint32)\n", "\n", - "marching_cubes_algorithm_seconds = (\n", - " perf_counter() - marching_cubes_start\n", - ")\n", - "marching_cubes_total_seconds = (\n", - " perf_counter() - marching_cubes_total_start\n", - ")\n", + "marching_cubes_algorithm_seconds = perf_counter() - marching_cubes_start\n", + "marching_cubes_total_seconds = perf_counter() - marching_cubes_total_start\n", "\n", "print(\"Regular-grid shape:\", comparison_grid_values.shape)\n", "print(\"Regular-grid sample count:\", comparison_grid_points.shape[0])\n", @@ -882,24 +859,24 @@ "print(\"Grid evaluation seconds:\", grid_evaluation_seconds)\n", "print(\"Marching-cubes algorithm seconds:\", marching_cubes_algorithm_seconds)\n", "print(\"Marching-cubes total seconds:\", marching_cubes_total_seconds)" - ], - "id": "9117b165c6dfa0e7", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "b479c5ef098030d6", + "metadata": {}, "source": [ "### Time RMT against the same field\n", "\n", "The RMT timing includes its field evaluations and mesh construction. The RBF evaluator is already cached for both methods. These single-run timings are useful for orientation, but they are not a rigorous performance benchmark; repeated warm runs and larger datasets would be needed for that." - ], - "id": "b479c5ef098030d6" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "c20c7765fc73c434", + "metadata": {}, + "outputs": [], "source": [ "rmt_total_start = perf_counter()\n", "\n", @@ -928,26 +905,26 @@ "print(\"RMT vertices:\", comparison_rmt_vertices.shape)\n", "print(\"RMT triangles:\", comparison_rmt_cells.shape)\n", "print(\"RMT total seconds:\", rmt_total_seconds)" - ], - "id": "c20c7765fc73c434", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "22b19f3ec6d5d959", + "metadata": {}, "source": [ "### Compare geometry and topology\n", "\n", "The triangle-quality score below is 1.0 for an equilateral triangle and approaches zero for a very thin or degenerate triangle. The 5th percentile is particularly useful because it describes the poor-quality end of each mesh without depending on only one extreme triangle.\n", "\n", "An edge used by more than two triangles is non-manifold. Boundary edges are expected here because the folded horizon reaches the sides of the extraction box and boundary closure was disabled. Therefore, a non-zero boundary-edge count does not by itself indicate a broken mesh." - ], - "id": "22b19f3ec6d5d959" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "412f81de43194b5a", + "metadata": {}, + "outputs": [], "source": [ "def calculate_mesh_statistics(\n", " mesh_vertices: np.ndarray,\n", @@ -1002,9 +979,7 @@ " \"triangles\": int(mesh_cells.shape[0]),\n", " \"elapsed_seconds\": float(elapsed_seconds),\n", " \"quality_mean\": float(np.mean(triangle_quality)),\n", - " \"quality_5th_percentile\": float(\n", - " np.percentile(triangle_quality, 5)\n", - " ),\n", + " \"quality_5th_percentile\": float(np.percentile(triangle_quality, 5)),\n", " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", " \"mean_absolute_isovalue_error\": float(\n", @@ -1021,19 +996,15 @@ " return statistics, fitted_mesh_values, mesh_true_vertical_errors\n", "\n", "\n", - "mc_statistics, mc_fitted_values, mc_true_vertical_errors = (\n", - " calculate_mesh_statistics(\n", - " mc_vertices,\n", - " mc_cells,\n", - " marching_cubes_total_seconds,\n", - " )\n", + "mc_statistics, mc_fitted_values, mc_true_vertical_errors = calculate_mesh_statistics(\n", + " mc_vertices,\n", + " mc_cells,\n", + " marching_cubes_total_seconds,\n", ")\n", - "rmt_statistics, rmt_fitted_values, rmt_true_vertical_errors = (\n", - " calculate_mesh_statistics(\n", - " comparison_rmt_vertices,\n", - " comparison_rmt_cells,\n", - " rmt_total_seconds,\n", - " )\n", + "rmt_statistics, rmt_fitted_values, rmt_true_vertical_errors = calculate_mesh_statistics(\n", + " comparison_rmt_vertices,\n", + " comparison_rmt_cells,\n", + " rmt_total_seconds,\n", ")\n", "\n", "comparison_rows = [\n", @@ -1055,24 +1026,24 @@ " mc_text = format(mc_statistics[key], value_format)\n", " rmt_text = format(rmt_statistics[key], value_format)\n", " print(f\"{label:38s} {mc_text:>18s} {rmt_text:>18s}\")" - ], - "id": "412f81de43194b5a", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "d85ca22ca792b883", + "metadata": {}, "source": [ "### View the two extraction results\n", "\n", "The surfaces should have the same large-scale folded geometry because they represent the same zero-isovalue of the same RBF. Differences to look for are triangulation density, skinny triangles, boundary shape, and small positional deviations introduced by RMT's regularisation and vertex clustering, though those things are difficult to discern when there are so many triangles." - ], - "id": "d85ca22ca792b883" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "c46d0141b8471c12", + "metadata": {}, + "outputs": [], "source": [ "figure = plt.figure(figsize=(14, 6))\n", "\n", @@ -1115,24 +1086,24 @@ "\n", "figure.tight_layout()\n", "plt.show()" - ], - "id": "c46d0141b8471c12", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "2852865e68b17840", + "metadata": {}, "source": [ "### Export both comparison surfaces to geoh5\n", "\n", "Both meshes are placed in the same group so they can be toggled, coloured, and compared in Geoscience ANALYST. Each receives the RBF value and known vertical error at its vertices." - ], - "id": "2852865e68b17840" + ] }, { - "metadata": {}, "cell_type": "code", + "execution_count": null, + "id": "ccfd27e9a9805429", + "metadata": {}, + "outputs": [], "source": [ "comparison_output_path = output_directory / (\n", " f\"marching_cubes_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", @@ -1173,22 +1144,19 @@ " )\n", "\n", "print(\"Created:\", comparison_output_path)" - ], - "id": "ccfd27e9a9805429", - "outputs": [], - "execution_count": null + ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "e66498617a269a8f", + "metadata": {}, "source": [ "### What this comparison does and does not establish\n", "\n", "This establishes whether RMT is a credible alternative mesh extractor for a field that already exists. It can compare mesh density, triangle quality, topology, approximate accuracy, and timing under controlled conditions.\n", "\n", "It does not establish support for geological orientation or exclusion constraints. Those affect how the scalar field is fitted from observations, which happens before either marching cubes or RMT is called." - ], - "id": "e66498617a269a8f" + ] } ], "metadata": { diff --git a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb index f42388c..5752b56 100644 --- a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb @@ -20,19 +20,21 @@ }, { "cell_type": "code", + "execution_count": null, "id": "d0d17fd27a2945db", "metadata": {}, + "outputs": [], "source": [ "from datetime import datetime\n", "from pathlib import Path\n", "from time import perf_counter\n", "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "from ferreus_rmt import BoundaryClosure, ClusterMethod, build_isosurface\n", "from geoh5py.groups import ContainerGroup\n", "from geoh5py.objects import Points, Surface\n", "from geoh5py.workspace import Workspace\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", "from scipy.interpolate import RegularGridInterpolator\n", "from scipy.sparse import coo_matrix\n", "from scipy.sparse.csgraph import connected_components\n", @@ -42,9 +44,7 @@ " entity_to_grid,\n", " extract_iso_surfaces,\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -60,8 +60,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "77ded4c6f1424af2", "metadata": {}, + "outputs": [], "source": [ "input_candidates = [\n", " Path.cwd() / \"geochem_pts\" / \"geochem_pts_example.geoh5\",\n", @@ -95,9 +97,7 @@ "\n", "print(\"Input:\", INPUT_PATH)\n", "print(\"Field and isovalue:\", DATA_NAME, CONTOUR_LEVEL)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -113,8 +113,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "90a877f44c8b44fb", "metadata": {}, + "outputs": [], "source": [ "grid_start = perf_counter()\n", "\n", @@ -123,9 +125,7 @@ " source_data = input_workspace.get_entity(DATA_NAME)[0]\n", "\n", " if source_object is None or source_data is None:\n", - " raise ValueError(\n", - " f\"Could not find {OBJECT_NAME!r} and {DATA_NAME!r}.\"\n", - " )\n", + " raise ValueError(f\"Could not find {OBJECT_NAME!r} and {DATA_NAME!r}.\")\n", "\n", " source_locations = np.asarray(\n", " source_object.locations,\n", @@ -152,10 +152,13 @@ "print(\"Observation count:\", source_locations.shape[0])\n", "print(\"XYZ span:\", np.ptp(source_locations, axis=0))\n", "print(\"Finite Cu observations:\", finite_source_values.size)\n", - "print(\"Cu observation percentiles:\", np.percentile(\n", - " finite_source_values,\n", - " [10, 25, 50, 75, 90, 95, 99],\n", - "))\n", + "print(\n", + " \"Cu observation percentiles:\",\n", + " np.percentile(\n", + " finite_source_values,\n", + " [10, 25, 50, 75, 90, 95, 99],\n", + " ),\n", + ")\n", "print(\"Grid shape:\", comparison_grid_values.shape)\n", "print(\n", " \"Supported grid nodes:\",\n", @@ -164,9 +167,7 @@ " comparison_grid_values.size,\n", ")\n", "print(\"Shared gridding time (seconds):\", grid_seconds)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -180,8 +181,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "a953d300954d4afa", "metadata": {}, + "outputs": [], "source": [ "marching_cubes_start = perf_counter()\n", "\n", @@ -200,9 +203,7 @@ "print(\"Marching-cubes vertices:\", mc_vertices.shape)\n", "print(\"Marching-cubes triangles:\", mc_cells.shape)\n", "print(\"Marching-cubes extraction time (seconds):\", marching_cubes_seconds)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -218,8 +219,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "da81ec4af37d4cd2", "metadata": {}, + "outputs": [], "source": [ "grid_interpolator = RegularGridInterpolator(\n", " tuple(comparison_grid),\n", @@ -272,14 +275,14 @@ " for axis in range(3)\n", " ]\n", " ).astype(np.float64)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "id": "aa3c77a2469f4b56", "metadata": {}, + "outputs": [], "source": [ "rmt_total_start = perf_counter()\n", "\n", @@ -324,9 +327,7 @@ "print(\"Seed detection time (seconds):\", seed_detection_seconds)\n", "print(\"RMT extraction time (seconds):\", rmt_extraction_seconds)\n", "print(\"RMT seed plus extraction time (seconds):\", rmt_total_seconds)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -342,8 +343,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "e14ca8c55c72457b", "metadata": {}, + "outputs": [], "source": [ "def calculate_mesh_statistics(\n", " mesh_vertices: np.ndarray,\n", @@ -409,9 +412,7 @@ " \"extraction_seconds\": float(extraction_seconds),\n", " \"surface_area\": float(np.sum(0.5 * twice_area)),\n", " \"quality_mean\": float(np.mean(triangle_quality)),\n", - " \"quality_5th_percentile\": float(\n", - " np.percentile(triangle_quality, 5)\n", - " ),\n", + " \"quality_5th_percentile\": float(np.percentile(triangle_quality, 5)),\n", " \"components\": int(component_count),\n", " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", @@ -454,12 +455,8 @@ "for label, key, value_format in comparison_rows:\n", " mc_text = format(mc_statistics[key], value_format)\n", " rmt_text = format(rmt_statistics[key], value_format)\n", - " print(f\"{label:40s} {mc_text:>18s} {rmt_text:>18s}\")\n", - "\n", - "\n" - ], - "outputs": [], - "execution_count": null + " print(f\"{label:40s} {mc_text:>18s} {rmt_text:>18s}\")" + ] }, { "cell_type": "markdown", @@ -473,8 +470,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "a85ce67ad448473a", "metadata": {}, + "outputs": [], "source": [ "mc_to_rmt_distances = cKDTree(rmt_vertices).query(mc_vertices)[0]\n", "rmt_to_mc_distances = cKDTree(mc_vertices).query(rmt_vertices)[0]\n", @@ -486,9 +485,7 @@ " \"Symmetric distance percentiles (50%, 90%, 95%, 100%):\",\n", " np.percentile(symmetric_distances, [50, 90, 95, 100]),\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -514,8 +511,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "21305e00d59c4e30", "metadata": {}, + "outputs": [], "source": [ "'''MAX_PREVIEW_TRIANGLES = 2_500\n", "\n", @@ -588,9 +587,7 @@ "figure.suptitle(f\"{DATA_NAME} = {CONTOUR_LEVEL:g} ppm\")\n", "figure.subplots_adjust(left=0.02, right=0.98, bottom=0.02, top=0.90)\n", "plt.show()'''" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -604,8 +601,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "62451e97bef04a06", "metadata": {}, + "outputs": [], "source": [ "# Recreate the output directory if previous results were deleted after\n", "# the configuration cell was run.\n", @@ -639,9 +638,7 @@ " output_mc_surface.add_data(\n", " {\n", " \"Interpolated Cu ppm\": {\"values\": mc_field_values},\n", - " \"Nearest RMT vertex distance\": {\n", - " \"values\": mc_to_rmt_distances\n", - " },\n", + " \"Nearest RMT vertex distance\": {\"values\": mc_to_rmt_distances},\n", " }\n", " )\n", "\n", @@ -655,16 +652,12 @@ " output_rmt_surface.add_data(\n", " {\n", " \"Interpolated Cu ppm\": {\"values\": rmt_field_values},\n", - " \"Nearest marching-cubes vertex distance\": {\n", - " \"values\": rmt_to_mc_distances\n", - " },\n", + " \"Nearest marching-cubes vertex distance\": {\"values\": rmt_to_mc_distances},\n", " }\n", " )\n", "\n", "print(\"Created:\", output_path)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", From 4a326c6317738436dab18605d54f712d523fc5b1 Mon Sep 17 00:00:00 2001 From: graeme Date: Thu, 20 Aug 2026 10:50:06 -0400 Subject: [PATCH 3/8] Fixes, add ferreus vs scipy interpolation comparison, add notes about next steps for constraints --- .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 417 +++++++++++++++--- 1 file changed, 351 insertions(+), 66 deletions(-) diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb index 684a97f..98b7e87 100644 --- a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -32,10 +32,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "c14d7389d994a354", "metadata": {}, - "outputs": [], "source": [ "import os\n", "import sys\n", @@ -58,7 +56,9 @@ "from geoh5py.groups import ContainerGroup\n", "from geoh5py.objects import Points, Surface\n", "from geoh5py.workspace import Workspace" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -68,10 +68,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "9be78b4d9cd5c492", "metadata": {}, - "outputs": [], "source": [ "CENTER = np.array([0.0, 0.0, 0.0], dtype=np.float64)\n", "RADIUS = 50.0\n", @@ -99,7 +97,9 @@ " )\n", "\n", " return distances - RADIUS, gradients" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -109,10 +109,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "8d2fa33ef1334bfa", "metadata": {}, - "outputs": [], "source": [ "seed_points = np.array(\n", " [\n", @@ -147,7 +145,9 @@ "\n", "print(\"Vertices:\", vertices.shape)\n", "print(\"Triangles:\", cells.shape)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -157,10 +157,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "c953cf37f21ad3a0", "metadata": {}, - "outputs": [], "source": [ "radii = np.linalg.norm(vertices - CENTER, axis=1)\n", "radius_errors = radii - RADIUS\n", @@ -175,7 +173,9 @@ "assert np.abs(radius_errors).max() < resolution\n", "\n", "# numbers should be close to the radius, but likely not exact because of the resolution of the lattice." - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -185,10 +185,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "8c91eaf424ec37fb", "metadata": {}, - "outputs": [], "source": [ "figure = plt.figure(figsize=(8, 7))\n", "axis = figure.add_subplot(111, projection=\"3d\")\n", @@ -209,7 +207,9 @@ "axis.set_title(\"Analytic sphere extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -219,10 +219,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "f39977f0e682fcb9", "metadata": {}, - "outputs": [], "source": [ "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", @@ -246,7 +244,9 @@ " )\n", "\n", "print(\"Created:\", output_path)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -267,10 +267,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "e8a9c39ec7391034", "metadata": {}, - "outputs": [], "source": [ "MODEL_EXTENTS = np.array(\n", " [-100.0, -100.0, -60.0, 100.0, 100.0, 60.0],\n", @@ -294,7 +292,9 @@ " )\n", "\n", " return points[:, 2] - horizon_z" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -304,10 +304,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "63f6e9858be0976d", "metadata": {}, - "outputs": [], "source": [ "random_generator = np.random.default_rng(42)\n", "\n", @@ -329,7 +327,9 @@ " \"to\",\n", " float(source_values.max()),\n", ")" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -339,10 +339,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "b82e4349d22b0c02", "metadata": {}, - "outputs": [], "source": [ "rbf_settings = InterpolantSettings(\n", " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", @@ -355,7 +353,9 @@ " values=source_values,\n", " interpolant_settings=rbf_settings,\n", ")" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -365,10 +365,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "f1900adf375eec99", "metadata": {}, - "outputs": [], "source": [ "source_predictions = np.asarray(rbf_model.evaluate_at_source(add_nugget=False)).reshape(\n", " -1\n", @@ -387,7 +385,9 @@ "\n", "# This demonstrates how closely the solved RBF reproduces the values used to fit it.\n", "# For a zero-nugget, synthetic problem, the residuals should be very small, on the order of machine precision." - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -397,10 +397,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "28d26151ceb1bbf9", "metadata": {}, - "outputs": [], "source": [ "number_of_validation_points = 200\n", "\n", @@ -428,7 +426,9 @@ "\n", "# The source residual measures fitting accuracy, i.e. an RMSE of 0.69 means approx 0.69 vertical units\n", "# The validation residual measures interpolation accuracy between observations" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -438,10 +438,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "3c481cec6ea9a739", "metadata": {}, - "outputs": [], "source": [ "folded_resolution = 5.0\n", "\n", @@ -462,7 +460,9 @@ ")\n", "\n", "rbf_model.build_evaluator(evaluator_extents)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -480,10 +480,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "16cbac4b3b5c9df8", "metadata": {}, - "outputs": [], "source": [ "seed_coordinates = np.linspace(-80.0, 80.0, 5)\n", "seed_x, seed_y = np.meshgrid(\n", @@ -508,7 +506,9 @@ " \"Maximum absolute seed field value:\",\n", " float(np.max(np.abs(true_geological_field(folded_seed_points)))),\n", ")" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -526,10 +526,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "7a79146bd28e5dd1", "metadata": {}, - "outputs": [], "source": [ "folded_isovalue = 0.0\n", "\n", @@ -555,7 +553,9 @@ "\n", "print(\"Folded surface vertices:\", folded_vertices.shape)\n", "print(\"Folded surface triangles:\", folded_cells.shape)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -574,10 +574,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "798b9ee4fbcc148e", "metadata": {}, - "outputs": [], "source": [ "fitted_values_at_vertices = np.asarray(\n", " rbf_model.evaluate(folded_vertices),\n", @@ -608,7 +606,9 @@ "assert folded_cells.ndim == 2 and folded_cells.shape[1] == 3\n", "assert np.all(np.isfinite(folded_vertices))\n", "assert np.all(np.isfinite(fitted_values_at_vertices))" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -622,10 +622,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "40448b38d7a150c4", "metadata": {}, - "outputs": [], "source": [ "near_horizon = np.abs(source_values) <= 15.0\n", "\n", @@ -671,7 +669,9 @@ "axis.set_title(\"RBF folded horizon extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -692,10 +692,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "156ec95b0d1e1877", "metadata": {}, - "outputs": [], "source": [ "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", @@ -743,7 +741,9 @@ " )\n", "\n", "print(\"Created:\", folded_output_path)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -784,10 +784,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "9117b165c6dfa0e7", "metadata": {}, - "outputs": [], "source": [ "from time import perf_counter\n", "\n", @@ -810,6 +808,23 @@ " comparison_grid[0],\n", " comparison_grid[1],\n", " comparison_grid[2],\n", + " indexing=\"ij\",\n", + ")\n", + "\n", + "# With indexing=\"ij\", array axes 0, 1, and 2 are X, Y, and Z.\n", + "# C-order ravel and reshape are therefore a matched pair. Using\n", + "# order=\"F\" here would reorder otherwise correctly paired values.\n", + "\n", + "expected_grid_shape = tuple(axis.size for axis in comparison_grid)\n", + "assert comparison_x.shape == expected_grid_shape\n", + "sample_index = tuple(min(2, axis.size - 1) for axis in comparison_grid)\n", + "assert np.allclose(\n", + " [\n", + " comparison_x[sample_index],\n", + " comparison_y[sample_index],\n", + " comparison_z[sample_index],\n", + " ],\n", + " [comparison_grid[axis][sample_index[axis]] for axis in range(3)],\n", ")\n", "comparison_grid_points = np.column_stack(\n", " [\n", @@ -822,10 +837,12 @@ "marching_cubes_total_start = perf_counter()\n", "\n", "grid_evaluation_start = perf_counter()\n", + "\n", "comparison_grid_values = np.asarray(\n", " rbf_model.evaluate_targets(comparison_grid_points),\n", " dtype=np.float64,\n", ").reshape(comparison_x.shape)\n", + "\n", "grid_evaluation_seconds = perf_counter() - grid_evaluation_start\n", "\n", "marching_cubes_start = perf_counter()\n", @@ -859,7 +876,9 @@ "print(\"Grid evaluation seconds:\", grid_evaluation_seconds)\n", "print(\"Marching-cubes algorithm seconds:\", marching_cubes_algorithm_seconds)\n", "print(\"Marching-cubes total seconds:\", marching_cubes_total_seconds)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -873,10 +892,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "c20c7765fc73c434", "metadata": {}, - "outputs": [], "source": [ "rmt_total_start = perf_counter()\n", "\n", @@ -905,7 +922,9 @@ "print(\"RMT vertices:\", comparison_rmt_vertices.shape)\n", "print(\"RMT triangles:\", comparison_rmt_cells.shape)\n", "print(\"RMT total seconds:\", rmt_total_seconds)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -921,10 +940,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "412f81de43194b5a", "metadata": {}, - "outputs": [], "source": [ "def calculate_mesh_statistics(\n", " mesh_vertices: np.ndarray,\n", @@ -1026,7 +1043,9 @@ " mc_text = format(mc_statistics[key], value_format)\n", " rmt_text = format(rmt_statistics[key], value_format)\n", " print(f\"{label:38s} {mc_text:>18s} {rmt_text:>18s}\")" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -1040,10 +1059,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "c46d0141b8471c12", "metadata": {}, - "outputs": [], "source": [ "figure = plt.figure(figsize=(14, 6))\n", "\n", @@ -1086,7 +1103,9 @@ "\n", "figure.tight_layout()\n", "plt.show()" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -1100,10 +1119,8 @@ }, { "cell_type": "code", - "execution_count": null, "id": "ccfd27e9a9805429", "metadata": {}, - "outputs": [], "source": [ "comparison_output_path = output_directory / (\n", " f\"marching_cubes_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", @@ -1144,7 +1161,9 @@ " )\n", "\n", "print(\"Created:\", comparison_output_path)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -1157,6 +1176,272 @@ "\n", "It does not establish support for geological orientation or exclusion constraints. Those affect how the scalar field is fitted from observations, which happens before either marching cubes or RMT is called." ] + }, + { + "cell_type": "markdown", + "id": "f91e0a6edbb84baa", + "metadata": {}, + "source": [ + "## Compare Ferreus and SciPy RBF interpolation\n", + "\n", + "The mesh extractor can only be as good as the scalar field it receives. This section therefore compares the Ferreus and SciPy interpolators before doing any isosurfacing. Both models receive exactly the same source observations and validation points.\n", + "\n", + "A fair comparison requires equivalent model choices. Ferreus uses its linear radial kernel, a linear drift, and zero nugget. SciPy uses its linear kernel, a degree-one polynomial, and zero smoothing. These settings are the closest direct match; changing only one model's kernel, polynomial term, or smoothing would compare different mathematical models rather than implementations.\n", + "\n", + "The two source counts below show a small problem and a moderately larger one. Timings are single runs, so treat them as directional rather than as a formal performance benchmark." + ] + }, + { + "cell_type": "code", + "id": "edaa15ec9a814d49", + "metadata": {}, + "source": [ + "from scipy.interpolate import RBFInterpolator as ScipyRBFInterpolator\n", + "from time import perf_counter\n", + "\n", + "# Compare a small interpolation problem with a moderately larger one\n", + "RBF_BENCHMARK_SIZES = (400, 4431)\n", + "RBF_BENCHMARK_VALIDATION_POINTS = 2000\n", + "\n", + "# Fixed seed for reproducibility\n", + "benchmark_rng = np.random.default_rng(1234567)\n", + "\n", + "# Points used to measure accuracy between observations\n", + "benchmark_validation_points = benchmark_rng.uniform(\n", + " low=MODEL_EXTENTS[:3],\n", + " high=MODEL_EXTENTS[3:],\n", + " size=(RBF_BENCHMARK_VALIDATION_POINTS, 3),\n", + ").astype(np.float64)\n", + "benchmark_validation_truth = true_geological_field(\n", + " benchmark_validation_points\n", + ")\n", + "\n", + "# Ferreus prepares a cached evaluator before RMT repeatedly requests field\n", + "# values. Padding includes the halo RMT may sample outside the model box (prevents an RMT query from falling outside the\n", + "# cached evaluator's allowable domain)\n", + "benchmark_evaluator_extents = MODEL_EXTENTS + np.array(\n", + " [-25.0, -25.0, -25.0, 25.0, 25.0, 25.0],\n", + " dtype=np.float64,\n", + ")\n", + "benchmark_rows = []\n", + "\n", + "for benchmark_size in RBF_BENCHMARK_SIZES:\n", + " # Generate one source data set to give the exact same observations to\n", + " # both libraries\n", + " benchmark_source_points = benchmark_rng.uniform(\n", + " low=MODEL_EXTENTS[:3],\n", + " high=MODEL_EXTENTS[3:],\n", + " size=(benchmark_size, 3),\n", + " ).astype(np.float64)\n", + " benchmark_source_values = true_geological_field(\n", + " benchmark_source_points\n", + " ).astype(np.float64)\n", + "\n", + " # Time construction and solution of the Ferreus model. linear kernel,\n", + " # linear drift, and zero nugget are matched to SciPy below\n", + " ferreus_start = perf_counter()\n", + " ferreus_benchmark_model = RBFInterpolator(\n", + " points=benchmark_source_points,\n", + " values=benchmark_source_values,\n", + " interpolant_settings=InterpolantSettings(\n", + " kernel_type=RBFKernelType.Linear,\n", + " drift=Drift.Linear,\n", + " nugget=0.0,\n", + " ),\n", + " )\n", + " ferreus_fit_seconds = perf_counter() - ferreus_start\n", + "\n", + " # Report evaluator preparation separately.\n", + " # RMT can use the cached evaluate_targets path multiple times.\n", + " ferreus_start = perf_counter()\n", + " ferreus_benchmark_model.build_evaluator(\n", + " benchmark_evaluator_extents\n", + " )\n", + " ferreus_setup_seconds = perf_counter() - ferreus_start\n", + "\n", + " # Evaluate unseen locations using the cached path used by RMT.\n", + " ferreus_start = perf_counter()\n", + " ferreus_validation_values = np.asarray(\n", + " ferreus_benchmark_model.evaluate_targets(\n", + " benchmark_validation_points\n", + " ),\n", + " dtype=np.float64,\n", + " ).reshape(-1)\n", + " ferreus_evaluate_seconds = perf_counter() - ferreus_start\n", + "\n", + " # Re-evaluate source points to test how closely the zero-nugget model\n", + " # reproduces the observations used to fit it\n", + " ferreus_source_values = np.asarray(\n", + " ferreus_benchmark_model.evaluate_at_source(add_nugget=False),\n", + " dtype=np.float64,\n", + " ).reshape(-1)\n", + "\n", + " # These are SciPy's closest equivalent settings. Matching the kernel,\n", + " # polynomial trend, and smoothing should help keep this an implementation comparison\n", + " # rather than a comparison of two different interpolation models.\n", + " scipy_start = perf_counter()\n", + " scipy_benchmark_model = ScipyRBFInterpolator(\n", + " benchmark_source_points,\n", + " benchmark_source_values,\n", + " kernel=\"linear\",\n", + " degree=1,\n", + " smoothing=0.0,\n", + " )\n", + " scipy_fit_seconds = perf_counter() - scipy_start\n", + "\n", + " # Evaluate exactly the same unseen validation locations with SciPy.\n", + " scipy_start = perf_counter()\n", + " scipy_validation_values = np.asarray(\n", + " scipy_benchmark_model(benchmark_validation_points),\n", + " dtype=np.float64,\n", + " ).reshape(-1)\n", + " scipy_evaluate_seconds = perf_counter() - scipy_start\n", + "\n", + " # This source-point call is kept outside the validation timing above.\n", + " scipy_source_values = np.asarray(\n", + " scipy_benchmark_model(benchmark_source_points),\n", + " dtype=np.float64,\n", + " ).reshape(-1)\n", + "\n", + " # Store three different accuracy checks:\n", + " # source RMSE measures reproduction of fitted observations;\n", + " # validation RMSE measures accuracy away from observations;\n", + " # between-model RMSE measures direct Ferreus/SciPy agreement.\n", + " benchmark_rows.append(\n", + " {\n", + " \"sources\": benchmark_size,\n", + " \"ferreus_fit_s\": ferreus_fit_seconds,\n", + " \"ferreus_setup_s\": ferreus_setup_seconds,\n", + " \"ferreus_eval_s\": ferreus_evaluate_seconds,\n", + " \"scipy_fit_s\": scipy_fit_seconds,\n", + " \"scipy_eval_s\": scipy_evaluate_seconds,\n", + " \"ferreus_source_rmse\": float(\n", + " np.sqrt(\n", + " np.mean(\n", + " (ferreus_source_values - benchmark_source_values) ** 2\n", + " )\n", + " )\n", + " ),\n", + " \"scipy_source_rmse\": float(\n", + " np.sqrt(\n", + " np.mean(\n", + " (scipy_source_values - benchmark_source_values) ** 2\n", + " )\n", + " )\n", + " ),\n", + " \"ferreus_validation_rmse\": float(\n", + " np.sqrt(\n", + " np.mean(\n", + " (\n", + " ferreus_validation_values\n", + " - benchmark_validation_truth\n", + " )\n", + " ** 2\n", + " )\n", + " )\n", + " ),\n", + " \"scipy_validation_rmse\": float(\n", + " np.sqrt(\n", + " np.mean(\n", + " (\n", + " scipy_validation_values\n", + " - benchmark_validation_truth\n", + " )\n", + " ** 2\n", + " )\n", + " )\n", + " ),\n", + " \"between_models_rmse\": float(\n", + " np.sqrt(\n", + " np.mean(\n", + " (\n", + " ferreus_validation_values\n", + " - scipy_validation_values\n", + " )\n", + " ** 2\n", + " )\n", + " )\n", + " ),\n", + " }\n", + " )\n", + "\n", + "# F and S = Ferreus and SciPy. Timing columns are elapsed seconds;\n", + "# the remaining columns are scalar-field RMSE values.\n", + "print(\n", + " f\"{'N':>6} {'F fit':>9} {'F setup':>9} {'F eval':>9} \"\n", + " f\"{'S fit':>9} {'S eval':>9} {'F val RMSE':>12} \"\n", + " f\"{'S val RMSE':>12} {'F-S RMSE':>12}\"\n", + ")\n", + "for row in benchmark_rows:\n", + " print(\n", + " f\"{row['sources']:6d} \"\n", + " f\"{row['ferreus_fit_s']:9.4f} \"\n", + " f\"{row['ferreus_setup_s']:9.4f} \"\n", + " f\"{row['ferreus_eval_s']:9.4f} \"\n", + " f\"{row['scipy_fit_s']:9.4f} \"\n", + " f\"{row['scipy_eval_s']:9.4f} \"\n", + " f\"{row['ferreus_validation_rmse']:12.6g} \"\n", + " f\"{row['scipy_validation_rmse']:12.6g} \"\n", + " f\"{row['between_models_rmse']:12.6g}\"\n", + " )\n", + "\n", + "# Source errors answer whether the unsmoothed models pass through their\n", + "# observations, which is different from prediction accuracy between them.\n", + "print(\"\\nSource-point RMSE (zero means exact interpolation):\")\n", + "for row in benchmark_rows:\n", + " print(\n", + " f\"N={row['sources']}: \"\n", + " f\"Ferreus={row['ferreus_source_rmse']:.6g}, \"\n", + " f\"SciPy={row['scipy_source_rmse']:.6g}\"\n", + " )\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "3f25676179f04fc1", + "metadata": {}, + "source": [ + "### How to interpret the RBF comparison\n", + "\n", + "- Source-point RMSE checks whether each zero-smoothing model reproduces the observations used to fit it.\n", + "- Validation RMSE compares each model with the known synthetic field at locations that were not used for fitting. This measures interpolation accuracy for this controlled example.\n", + "- F-S RMSE compares the two implementations directly. If it is tiny relative to the validation error, their matched settings are producing essentially the same scalar field.\n", + "- Ferreus setup is reported separately because build_evaluator prepares its cached evaluation path. RMT can then call evaluate_targets repeatedly. Comparing only this cached evaluation time with SciPy's evaluation time would omit that setup cost.\n", + "- These results test ordinary value interpolation only. They do not demonstrate orientation, contact, or inequality/exclusion-constraint support.\n", + "\n", + "SciPy documentation notes that a global RBF solve has memory requirements that grow quadratically with the number of source points. Its neighbors option can reduce that cost by using a local subset, but doing so changes the interpolation method and would no longer be a direct comparison with the global Ferreus model used here.\n", + "\n", + "Overall, the two options seem to produce similar interpolations. The interpolators reproduce their input observations with little difference between them. At a larger input size, ferreus has slightly weaker precision but both are good. Ferreus seems to be slightly faster to run, with better evaluation performance." + ] + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Constraints\n", + "\n", + "The following is my understanding of how different constraints can be handled, whether they are readily available, etc.\n", + "\n", + "### Current conclusion\n", + "\n", + "The public ferreus_rbf Python fitting API accepts coordinates and scalar values. It can calculate gradients from a fitted field, but it cannot accept measured gradients as fitting observations. It also does not expose an inequality-constrained solver. Therefore, contacts are directly supported, while geological orientation and inequality constraints (used for exclusion) are not directly supported. This limitation belongs to the field-fitting stage; ferreus_rmt can still extract a mesh from any scalar field supplied to it.\n", + "\n", + "### Workarounds which may be worth testing\n", + "\n", + "1. Represent an above/below inequality with a small positive or negative pseudo-value. This is only an approximation because it gives the point an exact value when only its sign is really known.\n", + "2. At a contact with a measured normal vector, add two offset points on opposite sides of the contact and assign them small negative and positive values. This encourages the fitted field to have the requested orientation. It is less reliable for orientations measured away from a known contact because the correct scalar level there is unknown.\n", + "\n", + "Both workarounds introduce choices such as the pseudo-value and offset distance. A result is only convincing if it remains reasonably stable when those choices are varied. They will be approximations.\n", + "\n", + "### Potential first steps to investigate\n", + "\n", + "Start with the existing synthetic folded horizon, where the correct answer is known, and compare: contacts only; contacts plus pseudo-inequalities; contacts plus orientation offset points; and all approximate constraints together. For each case, check contact RMSE, the percentage of above/below conditions satisfied, the angle between requested and predicted normals, and error against the known horizon. Export the resulting RMT surfaces to geoh5 for visual comparison.\n", + "\n", + "The workaround would be promising if it keeps contacts accurate, satisfies most inequalities, improves orientation agreement, and produces similar surfaces across reasonable parameter choices. My ballpark estimate of going down this path is somewhere around 1-2 days, depending on how well things work, what quirks or complications come up, things that I have to spend more time working through and reading about so I can properly understand them, etc." + ], + "id": "25c1e859ea5064ee" } ], "metadata": { From 2b15a78881acb26e7068d9a54f5b00a10f7c3247 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 20 Aug 2026 14:52:38 +0000 Subject: [PATCH 4/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 202 ++++++++---------- 1 file changed, 87 insertions(+), 115 deletions(-) diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb index 98b7e87..d74a375 100644 --- a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -32,8 +32,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "c14d7389d994a354", "metadata": {}, + "outputs": [], "source": [ "import os\n", "import sys\n", @@ -56,9 +58,7 @@ "from geoh5py.groups import ContainerGroup\n", "from geoh5py.objects import Points, Surface\n", "from geoh5py.workspace import Workspace" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -68,8 +68,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "9be78b4d9cd5c492", "metadata": {}, + "outputs": [], "source": [ "CENTER = np.array([0.0, 0.0, 0.0], dtype=np.float64)\n", "RADIUS = 50.0\n", @@ -97,9 +99,7 @@ " )\n", "\n", " return distances - RADIUS, gradients" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -109,8 +109,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "8d2fa33ef1334bfa", "metadata": {}, + "outputs": [], "source": [ "seed_points = np.array(\n", " [\n", @@ -145,9 +147,7 @@ "\n", "print(\"Vertices:\", vertices.shape)\n", "print(\"Triangles:\", cells.shape)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -157,8 +157,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "c953cf37f21ad3a0", "metadata": {}, + "outputs": [], "source": [ "radii = np.linalg.norm(vertices - CENTER, axis=1)\n", "radius_errors = radii - RADIUS\n", @@ -173,9 +175,7 @@ "assert np.abs(radius_errors).max() < resolution\n", "\n", "# numbers should be close to the radius, but likely not exact because of the resolution of the lattice." - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -185,8 +185,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "8c91eaf424ec37fb", "metadata": {}, + "outputs": [], "source": [ "figure = plt.figure(figsize=(8, 7))\n", "axis = figure.add_subplot(111, projection=\"3d\")\n", @@ -207,9 +209,7 @@ "axis.set_title(\"Analytic sphere extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -219,8 +219,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "f39977f0e682fcb9", "metadata": {}, + "outputs": [], "source": [ "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", @@ -244,9 +246,7 @@ " )\n", "\n", "print(\"Created:\", output_path)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -267,8 +267,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "e8a9c39ec7391034", "metadata": {}, + "outputs": [], "source": [ "MODEL_EXTENTS = np.array(\n", " [-100.0, -100.0, -60.0, 100.0, 100.0, 60.0],\n", @@ -292,9 +294,7 @@ " )\n", "\n", " return points[:, 2] - horizon_z" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -304,8 +304,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "63f6e9858be0976d", "metadata": {}, + "outputs": [], "source": [ "random_generator = np.random.default_rng(42)\n", "\n", @@ -327,9 +329,7 @@ " \"to\",\n", " float(source_values.max()),\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -339,8 +339,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "b82e4349d22b0c02", "metadata": {}, + "outputs": [], "source": [ "rbf_settings = InterpolantSettings(\n", " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", @@ -353,9 +355,7 @@ " values=source_values,\n", " interpolant_settings=rbf_settings,\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -365,8 +365,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "f1900adf375eec99", "metadata": {}, + "outputs": [], "source": [ "source_predictions = np.asarray(rbf_model.evaluate_at_source(add_nugget=False)).reshape(\n", " -1\n", @@ -385,9 +387,7 @@ "\n", "# This demonstrates how closely the solved RBF reproduces the values used to fit it.\n", "# For a zero-nugget, synthetic problem, the residuals should be very small, on the order of machine precision." - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -397,8 +397,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "28d26151ceb1bbf9", "metadata": {}, + "outputs": [], "source": [ "number_of_validation_points = 200\n", "\n", @@ -426,9 +428,7 @@ "\n", "# The source residual measures fitting accuracy, i.e. an RMSE of 0.69 means approx 0.69 vertical units\n", "# The validation residual measures interpolation accuracy between observations" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -438,8 +438,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "3c481cec6ea9a739", "metadata": {}, + "outputs": [], "source": [ "folded_resolution = 5.0\n", "\n", @@ -460,9 +462,7 @@ ")\n", "\n", "rbf_model.build_evaluator(evaluator_extents)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -480,8 +480,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "16cbac4b3b5c9df8", "metadata": {}, + "outputs": [], "source": [ "seed_coordinates = np.linspace(-80.0, 80.0, 5)\n", "seed_x, seed_y = np.meshgrid(\n", @@ -506,9 +508,7 @@ " \"Maximum absolute seed field value:\",\n", " float(np.max(np.abs(true_geological_field(folded_seed_points)))),\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -526,8 +526,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "7a79146bd28e5dd1", "metadata": {}, + "outputs": [], "source": [ "folded_isovalue = 0.0\n", "\n", @@ -553,9 +555,7 @@ "\n", "print(\"Folded surface vertices:\", folded_vertices.shape)\n", "print(\"Folded surface triangles:\", folded_cells.shape)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -574,8 +574,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "798b9ee4fbcc148e", "metadata": {}, + "outputs": [], "source": [ "fitted_values_at_vertices = np.asarray(\n", " rbf_model.evaluate(folded_vertices),\n", @@ -606,9 +608,7 @@ "assert folded_cells.ndim == 2 and folded_cells.shape[1] == 3\n", "assert np.all(np.isfinite(folded_vertices))\n", "assert np.all(np.isfinite(fitted_values_at_vertices))" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -622,8 +622,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "40448b38d7a150c4", "metadata": {}, + "outputs": [], "source": [ "near_horizon = np.abs(source_values) <= 15.0\n", "\n", @@ -669,9 +671,7 @@ "axis.set_title(\"RBF folded horizon extracted with Ferreus RMT\")\n", "\n", "plt.show()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -692,8 +692,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "156ec95b0d1e1877", "metadata": {}, + "outputs": [], "source": [ "output_directory = Path.cwd() / \"ferreus_output\"\n", "output_directory.mkdir(exist_ok=True)\n", @@ -741,9 +743,7 @@ " )\n", "\n", "print(\"Created:\", folded_output_path)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -784,8 +784,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "9117b165c6dfa0e7", "metadata": {}, + "outputs": [], "source": [ "from time import perf_counter\n", "\n", @@ -876,9 +878,7 @@ "print(\"Grid evaluation seconds:\", grid_evaluation_seconds)\n", "print(\"Marching-cubes algorithm seconds:\", marching_cubes_algorithm_seconds)\n", "print(\"Marching-cubes total seconds:\", marching_cubes_total_seconds)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -892,8 +892,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "c20c7765fc73c434", "metadata": {}, + "outputs": [], "source": [ "rmt_total_start = perf_counter()\n", "\n", @@ -922,9 +924,7 @@ "print(\"RMT vertices:\", comparison_rmt_vertices.shape)\n", "print(\"RMT triangles:\", comparison_rmt_cells.shape)\n", "print(\"RMT total seconds:\", rmt_total_seconds)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -940,8 +940,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "412f81de43194b5a", "metadata": {}, + "outputs": [], "source": [ "def calculate_mesh_statistics(\n", " mesh_vertices: np.ndarray,\n", @@ -1043,9 +1045,7 @@ " mc_text = format(mc_statistics[key], value_format)\n", " rmt_text = format(rmt_statistics[key], value_format)\n", " print(f\"{label:38s} {mc_text:>18s} {rmt_text:>18s}\")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1059,8 +1059,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "c46d0141b8471c12", "metadata": {}, + "outputs": [], "source": [ "figure = plt.figure(figsize=(14, 6))\n", "\n", @@ -1103,9 +1105,7 @@ "\n", "figure.tight_layout()\n", "plt.show()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1119,8 +1119,10 @@ }, { "cell_type": "code", + "execution_count": null, "id": "ccfd27e9a9805429", "metadata": {}, + "outputs": [], "source": [ "comparison_output_path = output_directory / (\n", " f\"marching_cubes_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", @@ -1161,9 +1163,7 @@ " )\n", "\n", "print(\"Created:\", comparison_output_path)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1193,12 +1193,16 @@ }, { "cell_type": "code", + "execution_count": null, "id": "edaa15ec9a814d49", "metadata": {}, + "outputs": [], "source": [ - "from scipy.interpolate import RBFInterpolator as ScipyRBFInterpolator\n", "from time import perf_counter\n", "\n", + "from scipy.interpolate import RBFInterpolator as ScipyRBFInterpolator\n", + "\n", + "\n", "# Compare a small interpolation problem with a moderately larger one\n", "RBF_BENCHMARK_SIZES = (400, 4431)\n", "RBF_BENCHMARK_VALIDATION_POINTS = 2000\n", @@ -1212,9 +1216,7 @@ " high=MODEL_EXTENTS[3:],\n", " size=(RBF_BENCHMARK_VALIDATION_POINTS, 3),\n", ").astype(np.float64)\n", - "benchmark_validation_truth = true_geological_field(\n", - " benchmark_validation_points\n", - ")\n", + "benchmark_validation_truth = true_geological_field(benchmark_validation_points)\n", "\n", "# Ferreus prepares a cached evaluator before RMT repeatedly requests field\n", "# values. Padding includes the halo RMT may sample outside the model box (prevents an RMT query from falling outside the\n", @@ -1233,9 +1235,9 @@ " high=MODEL_EXTENTS[3:],\n", " size=(benchmark_size, 3),\n", " ).astype(np.float64)\n", - " benchmark_source_values = true_geological_field(\n", - " benchmark_source_points\n", - " ).astype(np.float64)\n", + " benchmark_source_values = true_geological_field(benchmark_source_points).astype(\n", + " np.float64\n", + " )\n", "\n", " # Time construction and solution of the Ferreus model. linear kernel,\n", " # linear drift, and zero nugget are matched to SciPy below\n", @@ -1254,17 +1256,13 @@ " # Report evaluator preparation separately.\n", " # RMT can use the cached evaluate_targets path multiple times.\n", " ferreus_start = perf_counter()\n", - " ferreus_benchmark_model.build_evaluator(\n", - " benchmark_evaluator_extents\n", - " )\n", + " ferreus_benchmark_model.build_evaluator(benchmark_evaluator_extents)\n", " ferreus_setup_seconds = perf_counter() - ferreus_start\n", "\n", " # Evaluate unseen locations using the cached path used by RMT.\n", " ferreus_start = perf_counter()\n", " ferreus_validation_values = np.asarray(\n", - " ferreus_benchmark_model.evaluate_targets(\n", - " benchmark_validation_points\n", - " ),\n", + " ferreus_benchmark_model.evaluate_targets(benchmark_validation_points),\n", " dtype=np.float64,\n", " ).reshape(-1)\n", " ferreus_evaluate_seconds = perf_counter() - ferreus_start\n", @@ -1316,50 +1314,26 @@ " \"scipy_fit_s\": scipy_fit_seconds,\n", " \"scipy_eval_s\": scipy_evaluate_seconds,\n", " \"ferreus_source_rmse\": float(\n", - " np.sqrt(\n", - " np.mean(\n", - " (ferreus_source_values - benchmark_source_values) ** 2\n", - " )\n", - " )\n", + " np.sqrt(np.mean((ferreus_source_values - benchmark_source_values) ** 2))\n", " ),\n", " \"scipy_source_rmse\": float(\n", - " np.sqrt(\n", - " np.mean(\n", - " (scipy_source_values - benchmark_source_values) ** 2\n", - " )\n", - " )\n", + " np.sqrt(np.mean((scipy_source_values - benchmark_source_values) ** 2))\n", " ),\n", " \"ferreus_validation_rmse\": float(\n", " np.sqrt(\n", " np.mean(\n", - " (\n", - " ferreus_validation_values\n", - " - benchmark_validation_truth\n", - " )\n", - " ** 2\n", + " (ferreus_validation_values - benchmark_validation_truth) ** 2\n", " )\n", " )\n", " ),\n", " \"scipy_validation_rmse\": float(\n", " np.sqrt(\n", - " np.mean(\n", - " (\n", - " scipy_validation_values\n", - " - benchmark_validation_truth\n", - " )\n", - " ** 2\n", - " )\n", + " np.mean((scipy_validation_values - benchmark_validation_truth) ** 2)\n", " )\n", " ),\n", " \"between_models_rmse\": float(\n", " np.sqrt(\n", - " np.mean(\n", - " (\n", - " ferreus_validation_values\n", - " - scipy_validation_values\n", - " )\n", - " ** 2\n", - " )\n", + " np.mean((ferreus_validation_values - scipy_validation_values) ** 2)\n", " )\n", " ),\n", " }\n", @@ -1393,10 +1367,8 @@ " f\"N={row['sources']}: \"\n", " f\"Ferreus={row['ferreus_source_rmse']:.6g}, \"\n", " f\"SciPy={row['scipy_source_rmse']:.6g}\"\n", - " )\n" - ], - "outputs": [], - "execution_count": null + " )" + ] }, { "cell_type": "markdown", @@ -1417,8 +1389,9 @@ ] }, { - "metadata": {}, "cell_type": "markdown", + "id": "25c1e859ea5064ee", + "metadata": {}, "source": [ "## Constraints\n", "\n", @@ -1440,8 +1413,7 @@ "Start with the existing synthetic folded horizon, where the correct answer is known, and compare: contacts only; contacts plus pseudo-inequalities; contacts plus orientation offset points; and all approximate constraints together. For each case, check contact RMSE, the percentage of above/below conditions satisfied, the angle between requested and predicted normals, and error against the known horizon. Export the resulting RMT surfaces to geoh5 for visual comparison.\n", "\n", "The workaround would be promising if it keeps contacts accurate, satisfies most inequalities, improves orientation agreement, and produces similar surfaces across reasonable parameter choices. My ballpark estimate of going down this path is somewhere around 1-2 days, depending on how well things work, what quirks or complications come up, things that I have to spend more time working through and reading about so I can properly understand them, etc." - ], - "id": "25c1e859ea5064ee" + ] } ], "metadata": { From cb359ca656334c0bbf65bdbf98c0f93a02f8e801 Mon Sep 17 00:00:00 2001 From: domfournier Date: Fri, 21 Aug 2026 10:36:07 -0700 Subject: [PATCH 5/8] Add notebooks to lfs --- .gitattributes | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitattributes b/.gitattributes index 6b1749a..42f9d0d 100644 --- a/.gitattributes +++ b/.gitattributes @@ -2,3 +2,4 @@ *.geoh5 filter=lfs diff=lfs merge=lfs -text *.h5 filter=lfs diff=lfs merge=lfs -text *.min.js filter=lfs diff=lfs merge=lfs -text +*.ipynb filter=lfs diff=lfs merge=lfs -text From 85d6c2d34922332eb683b4eba421b4533a5a1d40 Mon Sep 17 00:00:00 2001 From: domfournier Date: Fri, 21 Aug 2026 10:37:57 -0700 Subject: [PATCH 6/8] Move ipynb to assets --- .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 3 + .../ferreus_rmt_real_geochem_comparison.ipynb | 3 + .../ferreus/geochem_pts_example.geoh5 | 3 + .../ferreus/ferreus_rmt_geoh5_poc.ipynb | 1443 +---------------- .../ferreus_rmt_real_geochem_comparison.ipynb | 703 +------- 5 files changed, 15 insertions(+), 2140 deletions(-) create mode 100644 surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb create mode 100644 surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb create mode 100644 surface_apps-assets/ferreus/geochem_pts_example.geoh5 diff --git a/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb new file mode 100644 index 0000000..8073b52 --- /dev/null +++ b/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:91a597b583533d9228fa34d9b6114163a7c430583f42b38ac0102cd1315f6838 +size 947931 diff --git a/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb new file mode 100644 index 0000000..8d00c73 --- /dev/null +++ b/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0a732fc34d583d75ec4fbb4854e1fa1a269b1762c0cdf9c35b1317d2f303b697 +size 38221 diff --git a/surface_apps-assets/ferreus/geochem_pts_example.geoh5 b/surface_apps-assets/ferreus/geochem_pts_example.geoh5 new file mode 100644 index 0000000..88efc90 --- /dev/null +++ b/surface_apps-assets/ferreus/geochem_pts_example.geoh5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7336127becd1c8fa141a6564049251c6db7dc51a14687eff95e120f5b6593715 +size 1441792 diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb index d74a375..86859aa 100644 --- a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -1,1440 +1,3 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "52652cf59cffb624", - "metadata": { - "ExecuteTime": { - "end_time": "2026-08-10T15:19:45.243382800Z", - "start_time": "2026-08-10T15:19:45.191472Z" - }, - "collapsed": true - }, - "source": [ - "# Ferreus RMT to geoh5 proof of concept\n", - "\n", - "Objectives:\n", - "\n", - "Explore if the packages allows for\n", - "\n", - " exclusion constraints\n", - "\n", - " orientation constraints\n", - "\n", - "POC with notebook creating a surface transferred to geoh5" - ] - }, - { - "cell_type": "markdown", - "id": "de308ccc50aea28d", - "metadata": {}, - "source": "" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c14d7389d994a354", - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import sys\n", - "from datetime import datetime\n", - "from pathlib import Path\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "from ferreus_rbf import RBFInterpolator\n", - "from ferreus_rbf.interpolant_config import (\n", - " Drift,\n", - " InterpolantSettings,\n", - " RBFKernelType,\n", - ")\n", - "from ferreus_rmt import (\n", - " BoundaryClosure,\n", - " ClusterMethod,\n", - " build_isosurface,\n", - ")\n", - "from geoh5py.groups import ContainerGroup\n", - "from geoh5py.objects import Points, Surface\n", - "from geoh5py.workspace import Workspace" - ] - }, - { - "cell_type": "markdown", - "id": "a8c8dc3ebfbae12c", - "metadata": {}, - "source": "### define a known sphere" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9be78b4d9cd5c492", - "metadata": {}, - "outputs": [], - "source": [ - "CENTER = np.array([0.0, 0.0, 0.0], dtype=np.float64)\n", - "RADIUS = 50.0\n", - "\n", - "\n", - "def sphere_values(points: np.ndarray) -> np.ndarray:\n", - " \"\"\"negative inside the sphere, zero on it, positive outside.\"\"\"\n", - "\n", - " offsets = points - CENTER\n", - " return np.linalg.norm(offsets, axis=1) - RADIUS\n", - "\n", - "\n", - "def sphere_values_and_gradients(\n", - " points: np.ndarray,\n", - ") -> tuple[np.ndarray, np.ndarray]:\n", - " \"\"\"return values and analytic gradients.\"\"\"\n", - " offsets = points - CENTER\n", - " distances = np.linalg.norm(offsets, axis=1)\n", - "\n", - " gradients = np.divide(\n", - " offsets,\n", - " distances[:, None],\n", - " out=np.zeros_like(offsets),\n", - " where=distances[:, None] > 0.0,\n", - " )\n", - "\n", - " return distances - RADIUS, gradients" - ] - }, - { - "cell_type": "markdown", - "id": "bf48bdde9dab3801", - "metadata": {}, - "source": "### extract the surface with rmt" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8d2fa33ef1334bfa", - "metadata": {}, - "outputs": [], - "source": [ - "seed_points = np.array(\n", - " [\n", - " [50.0, 0.0, 0.0],\n", - " [-50.0, 0.0, 0.0],\n", - " ],\n", - " dtype=np.float64,\n", - ")\n", - "\n", - "# Define the axis-aligned bounding box extents to extract the isosurface within:\n", - "# [xmin, ymin, zmin, xmax, ymax, zmax].\n", - "extents = np.array(\n", - " [-75.0, -75.0, -75.0, 75.0, 75.0, 75.0],\n", - " dtype=np.float64,\n", - ")\n", - "\n", - "\n", - "resolution = 5.0\n", - "isovalue = 0.0\n", - "\n", - "mesh = build_isosurface(\n", - " seed_points, # approx. locations of the surface to seed the wavefront\n", - " extents, # the axis-aligned bounding box limiting extraction\n", - " resolution, # lattice spacing in coordinate units - smaller means more detail\n", - " isovalue, # the scalar value defining the surface\n", - " surface_fn=sphere_values, # surface_fn (sphere_values) is the scalar function that defines the geometry\n", - " gradient_fn=sphere_values_and_gradients, #\n", - ")\n", - "\n", - "vertices = np.asarray(mesh.vertices)\n", - "cells = np.asarray(mesh.facets, dtype=np.uint32)\n", - "\n", - "print(\"Vertices:\", vertices.shape)\n", - "print(\"Triangles:\", cells.shape)" - ] - }, - { - "cell_type": "markdown", - "id": "4178f0ab56928714", - "metadata": {}, - "source": "### validate the result" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c953cf37f21ad3a0", - "metadata": {}, - "outputs": [], - "source": [ - "radii = np.linalg.norm(vertices - CENTER, axis=1)\n", - "radius_errors = radii - RADIUS\n", - "\n", - "print(\"Minimum radius:\", radii.min())\n", - "print(\"Maximum radius:\", radii.max())\n", - "print(\"Mean absolute radius error:\", np.abs(radius_errors).mean())\n", - "print(\"Maximum absolute radius error:\", np.abs(radius_errors).max())\n", - "\n", - "assert cells.ndim == 2 and cells.shape[1] == 3\n", - "assert vertices.ndim == 2 and vertices.shape[1] == 3\n", - "assert np.abs(radius_errors).max() < resolution\n", - "\n", - "# numbers should be close to the radius, but likely not exact because of the resolution of the lattice." - ] - }, - { - "cell_type": "markdown", - "id": "3b1c992eb5ef149d", - "metadata": {}, - "source": "### show the sphere we've constructed" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8c91eaf424ec37fb", - "metadata": {}, - "outputs": [], - "source": [ - "figure = plt.figure(figsize=(8, 7))\n", - "axis = figure.add_subplot(111, projection=\"3d\")\n", - "\n", - "axis.plot_trisurf(\n", - " vertices[:, 0],\n", - " vertices[:, 1],\n", - " vertices[:, 2],\n", - " triangles=cells,\n", - " linewidth=0.05,\n", - " alpha=0.85,\n", - ")\n", - "\n", - "axis.set_xlabel(\"X\")\n", - "axis.set_ylabel(\"Y\")\n", - "axis.set_zlabel(\"Z\")\n", - "axis.set_box_aspect(np.ptp(vertices, axis=0))\n", - "axis.set_title(\"Analytic sphere extracted with Ferreus RMT\")\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "ba702ddc81228ad6", - "metadata": {}, - "source": "### Export to geoh5" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f39977f0e682fcb9", - "metadata": {}, - "outputs": [], - "source": [ - "output_directory = Path.cwd() / \"ferreus_output\"\n", - "output_directory.mkdir(exist_ok=True)\n", - "\n", - "output_path = output_directory / (\n", - " f\"ferreus_rmt_poc_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", - ")\n", - "\n", - "with Workspace.create(output_path) as workspace:\n", - " group = ContainerGroup.create(\n", - " workspace,\n", - " name=\"Ferreus RMT POC\",\n", - " )\n", - "\n", - " Surface.create(\n", - " workspace,\n", - " name=\"RMT analytic sphere\",\n", - " vertices=vertices,\n", - " cells=cells,\n", - " parent=group,\n", - " )\n", - "\n", - "print(\"Created:\", output_path)" - ] - }, - { - "cell_type": "markdown", - "id": "1a8d7c5fad014c7f", - "metadata": {}, - "source": [ - "## Fit a synthetic folded horizon with ferreus_rbf\n", - "\n", - "The analytic sphere proved that `ferreus_rmt` can extract a known scalar\n", - "field and that its mesh can be exported to geoH5.\n", - "\n", - "This section replaces the analytic sphere field with an RBF field fitted\n", - "from scattered scalar observations.\n", - "\n", - "A known analytic geological field is retained as ground truth, allowing\n", - "the interpolation and resulting surface errors to be measured." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e8a9c39ec7391034", - "metadata": {}, - "outputs": [], - "source": [ - "MODEL_EXTENTS = np.array(\n", - " [-100.0, -100.0, -60.0, 100.0, 100.0, 60.0],\n", - " dtype=np.float64,\n", - ")\n", - "\n", - "\n", - "def true_horizon_elevation(\n", - " x: np.ndarray,\n", - " y: np.ndarray,\n", - ") -> np.ndarray:\n", - " \"\"\"Elevation of the synthetic folded horizon.\"\"\"\n", - " return 20.0 * np.sin(x / 35.0) + 10.0 * np.cos(y / 40.0)\n", - "\n", - "\n", - "def true_geological_field(points: np.ndarray) -> np.ndarray:\n", - " \"\"\"Signed vertical separation from the synthetic horizon.\"\"\"\n", - " horizon_z = true_horizon_elevation(\n", - " points[:, 0],\n", - " points[:, 1],\n", - " )\n", - "\n", - " return points[:, 2] - horizon_z" - ] - }, - { - "cell_type": "markdown", - "id": "55396a8e404a0f23", - "metadata": {}, - "source": "### Generate scattered scalar observations" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "63f6e9858be0976d", - "metadata": {}, - "outputs": [], - "source": [ - "random_generator = np.random.default_rng(42)\n", - "\n", - "number_of_source_points = 400\n", - "\n", - "source_points = random_generator.uniform(\n", - " low=MODEL_EXTENTS[:3],\n", - " high=MODEL_EXTENTS[3:],\n", - " size=(number_of_source_points, 3),\n", - ").astype(np.float64)\n", - "\n", - "source_values = true_geological_field(source_points).astype(np.float64)\n", - "\n", - "print(\"Source points:\", source_points.shape)\n", - "print(\"Source values:\", source_values.shape)\n", - "print(\n", - " \"Value range:\",\n", - " float(source_values.min()),\n", - " \"to\",\n", - " float(source_values.max()),\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "203c7442a4c55fec", - "metadata": {}, - "source": "### Fit the RBF field" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b82e4349d22b0c02", - "metadata": {}, - "outputs": [], - "source": [ - "rbf_settings = InterpolantSettings(\n", - " kernel_type=RBFKernelType.Linear, # selects the radial kernel used to interpolate between observations.\n", - " drift=Drift.Linear, # adds a broad linear trend, useful here because the field contains a strong linear z component.\n", - " nugget=0.0, # 0.0 requests an unsmoothed fit. Fine for synthetic observations as they contain no noise.\n", - ")\n", - "\n", - "rbf_model = RBFInterpolator(\n", - " points=source_points,\n", - " values=source_values,\n", - " interpolant_settings=rbf_settings,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "6b8e055979af129a", - "metadata": {}, - "source": "### Check the training residuals" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f1900adf375eec99", - "metadata": {}, - "outputs": [], - "source": [ - "source_predictions = np.asarray(rbf_model.evaluate_at_source(add_nugget=False)).reshape(\n", - " -1\n", - ")\n", - "\n", - "source_residuals = source_predictions - source_values\n", - "\n", - "print(\n", - " \"Source RMSE:\",\n", - " float(np.sqrt(np.mean(source_residuals**2))),\n", - ")\n", - "print(\n", - " \"Maximum source error:\",\n", - " float(np.max(np.abs(source_residuals))),\n", - ")\n", - "\n", - "# This demonstrates how closely the solved RBF reproduces the values used to fit it.\n", - "# For a zero-nugget, synthetic problem, the residuals should be very small, on the order of machine precision." - ] - }, - { - "cell_type": "markdown", - "id": "e25a69ffedc71d99", - "metadata": {}, - "source": "### Test at locations not used for fitting" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "28d26151ceb1bbf9", - "metadata": {}, - "outputs": [], - "source": [ - "number_of_validation_points = 200\n", - "\n", - "validation_points = random_generator.uniform(\n", - " low=MODEL_EXTENTS[:3],\n", - " high=MODEL_EXTENTS[3:],\n", - " size=(number_of_validation_points, 3),\n", - ").astype(np.float64)\n", - "\n", - "true_validation_values = true_geological_field(validation_points)\n", - "\n", - "predicted_validation_values = np.asarray(rbf_model.evaluate(validation_points)).reshape(\n", - " -1\n", - ")\n", - "\n", - "validation_residuals = predicted_validation_values - true_validation_values\n", - "\n", - "validation_rmse = np.sqrt(np.mean(validation_residuals**2))\n", - "\n", - "print(\"Validation RMSE:\", float(validation_rmse))\n", - "print(\n", - " \"Maximum validation error:\",\n", - " float(np.max(np.abs(validation_residuals))),\n", - ")\n", - "\n", - "# The source residual measures fitting accuracy, i.e. an RMSE of 0.69 means approx 0.69 vertical units\n", - "# The validation residual measures interpolation accuracy between observations" - ] - }, - { - "cell_type": "markdown", - "id": "adcfe7b8d0d585b4", - "metadata": {}, - "source": "### Prepare the RBF for repeated RMT evaluations" - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3c481cec6ea9a739", - "metadata": {}, - "outputs": [], - "source": [ - "folded_resolution = 5.0\n", - "\n", - "# RMT evaluates a small halo of lattice positions around the extraction\n", - "# box. Give the cached RBF evaluator enough room for those calls.\n", - "evaluator_padding = 5.0 * folded_resolution\n", - "\n", - "evaluator_extents = MODEL_EXTENTS + np.array(\n", - " [\n", - " -evaluator_padding,\n", - " -evaluator_padding,\n", - " -evaluator_padding,\n", - " evaluator_padding,\n", - " evaluator_padding,\n", - " evaluator_padding,\n", - " ],\n", - " dtype=np.float64,\n", - ")\n", - "\n", - "rbf_model.build_evaluator(evaluator_extents)" - ] - }, - { - "cell_type": "markdown", - "id": "e674a1b0a8c42db5", - "metadata": {}, - "source": [ - "### Choose seed points for the folded horizon\n", - "\n", - "RMT is a surface-following algorithm. This means it does not scan every cell in the model volume. Instead, it starts near the requested isosurface at the seed points and follows the connected surface away from them.\n", - "\n", - "These seeds are algorithm inputs, not extra observations used to fit the RBF. In this synthetic example the true horizon is known, so a small grid of exact surface locations is convenient. In a real geological workflow, observed contact points with scalar value zero would be natural seeds.\n", - "\n", - "Multiple seeds are useful when a model may contain disconnected surface components. A component that has no nearby seed may not be discovered." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "16cbac4b3b5c9df8", - "metadata": {}, - "outputs": [], - "source": [ - "seed_coordinates = np.linspace(-80.0, 80.0, 5)\n", - "seed_x, seed_y = np.meshgrid(\n", - " seed_coordinates,\n", - " seed_coordinates,\n", - " indexing=\"xy\",\n", - ")\n", - "\n", - "folded_seed_points = np.column_stack(\n", - " [\n", - " seed_x.ravel(),\n", - " seed_y.ravel(),\n", - " true_horizon_elevation(\n", - " seed_x.ravel(),\n", - " seed_y.ravel(),\n", - " ),\n", - " ]\n", - ").astype(np.float64)\n", - "\n", - "print(\"RMT seed points shape:\", folded_seed_points.shape)\n", - "print(\n", - " \"Maximum absolute seed field value:\",\n", - " float(np.max(np.abs(true_geological_field(folded_seed_points)))),\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "f25538fc35ed124b", - "metadata": {}, - "source": [ - "### Extract the fitted zero-isosurface with RMT\n", - "\n", - "The requested surface is where the fitted scalar field equals folded_isovalue, here zero. Because the synthetic field is z - horizon_z, zero represents the folded horizon, negative values lie below it, and positive values lie above it.\n", - "\n", - "folded_resolution is the spacing of RMT's internal sampling lattice in the same coordinate units as X, Y, and Z. 5.0 means samples are approximately five units apart. A smaller value can preserve finer geometry, but requires more field evaluations and usually produces more triangles.\n", - "\n", - "evaluate_targets and evaluate_targets_with_gradients use the cached evaluator prepared above. The gradients help RMT project its working vertices toward the requested scalar value; they are not geological orientation observations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7a79146bd28e5dd1", - "metadata": {}, - "outputs": [], - "source": [ - "folded_isovalue = 0.0\n", - "\n", - "folded_mesh = build_isosurface(\n", - " seed_points=folded_seed_points,\n", - " extents=MODEL_EXTENTS,\n", - " resolution=folded_resolution,\n", - " isovalue=folded_isovalue,\n", - " surface_fn=rbf_model.evaluate_targets,\n", - " gradient_fn=rbf_model.evaluate_targets_with_gradients,\n", - " cluster_method=ClusterMethod.CurvatureWeighted,\n", - " boundary_closure=BoundaryClosure.None_,\n", - ")\n", - "\n", - "folded_vertices = np.asarray(\n", - " folded_mesh.vertices,\n", - " dtype=np.float64,\n", - ")\n", - "folded_cells = np.asarray(\n", - " folded_mesh.facets,\n", - " dtype=np.uint32,\n", - ")\n", - "\n", - "print(\"Folded surface vertices:\", folded_vertices.shape)\n", - "print(\"Folded surface triangles:\", folded_cells.shape)" - ] - }, - { - "cell_type": "markdown", - "id": "e203369f6da5b01f", - "metadata": {}, - "source": [ - "### Validate both parts of the workflow\n", - "\n", - "There are two distinct errors worth checking:\n", - "\n", - "1. **RMT extraction error:** evaluate the fitted RBF at every mesh vertex. These values should be close to the requested isovalue of zero. This checks whether RMT followed the field it was given.\n", - "2. **Modelling error:** evaluate the known analytic field at every mesh vertex. Here that value is also the signed vertical distance from the true horizon. This checks how closely the RBF-derived surface matches the ground truth between the observations.\n", - "\n", - "Keeping these checks separate matters as a mesh can follow an imperfect RBF field very accurately." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "798b9ee4fbcc148e", - "metadata": {}, - "outputs": [], - "source": [ - "fitted_values_at_vertices = np.asarray(\n", - " rbf_model.evaluate(folded_vertices),\n", - " dtype=np.float64,\n", - ").reshape(-1)\n", - "\n", - "rmt_isovalue_errors = fitted_values_at_vertices - folded_isovalue\n", - "\n", - "# For this particular synthetic field, this is the signed vertical\n", - "# separation between a mesh vertex and the true folded horizon.\n", - "true_vertical_errors = true_geological_field(folded_vertices)\n", - "absolute_vertical_errors = np.abs(true_vertical_errors)\n", - "\n", - "print(\n", - " \"Maximum absolute RMT isovalue error:\",\n", - " float(np.max(np.abs(rmt_isovalue_errors))),\n", - ")\n", - "print(\n", - " \"Mean absolute vertical surface error:\",\n", - " float(np.mean(absolute_vertical_errors)),\n", - ")\n", - "print(\n", - " \"Vertical error percentiles (50%, 90%, 95%, 100%):\",\n", - " np.percentile(absolute_vertical_errors, [50, 90, 95, 100]),\n", - ")\n", - "\n", - "assert folded_vertices.ndim == 2 and folded_vertices.shape[1] == 3\n", - "assert folded_cells.ndim == 2 and folded_cells.shape[1] == 3\n", - "assert np.all(np.isfinite(folded_vertices))\n", - "assert np.all(np.isfinite(fitted_values_at_vertices))" - ] - }, - { - "cell_type": "markdown", - "id": "be32da3dd2532365", - "metadata": {}, - "source": [ - "### Inspect the fitted surface and nearby observations\n", - "\n", - "Only observations reasonably close to the horizon are plotted, which keeps the diagnostic readable. Their colours show the signed scalar values supplied to the RBF: blue points are below the horizon, red points are above it, and values near zero lie close to the interpreted contact." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "40448b38d7a150c4", - "metadata": {}, - "outputs": [], - "source": [ - "near_horizon = np.abs(source_values) <= 15.0\n", - "\n", - "figure = plt.figure(figsize=(10, 8))\n", - "axis = figure.add_subplot(111, projection=\"3d\")\n", - "\n", - "axis.plot_trisurf(\n", - " folded_vertices[:, 0],\n", - " folded_vertices[:, 1],\n", - " folded_vertices[:, 2],\n", - " triangles=folded_cells,\n", - " color=\"lightsteelblue\",\n", - " linewidth=0.05,\n", - " alpha=0.70,\n", - ")\n", - "\n", - "observation_plot = axis.scatter(\n", - " source_points[near_horizon, 0],\n", - " source_points[near_horizon, 1],\n", - " source_points[near_horizon, 2],\n", - " c=source_values[near_horizon],\n", - " cmap=\"coolwarm\",\n", - " vmin=-15.0,\n", - " vmax=15.0,\n", - " s=20,\n", - " edgecolors=\"black\",\n", - " linewidths=0.2,\n", - ")\n", - "\n", - "figure.colorbar(\n", - " observation_plot,\n", - " ax=axis,\n", - " shrink=0.65,\n", - " pad=0.10,\n", - " label=\"Observed signed field value\",\n", - ")\n", - "\n", - "axis.set_xlabel(\"X\")\n", - "axis.set_ylabel(\"Y\")\n", - "axis.set_zlabel(\"Z\")\n", - "axis.set_box_aspect(MODEL_EXTENTS[3:] - MODEL_EXTENTS[:3])\n", - "axis.view_init(elev=25, azim=-55)\n", - "axis.set_title(\"RBF folded horizon extracted with Ferreus RMT\")\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "9eab773cbf3c33cd", - "metadata": {}, - "source": [ - "### Export the folded surface and its observations to geoh5\n", - "\n", - "The output contains:\n", - "\n", - "- a Points object holding every scalar observation used to fit the RBF;\n", - "- the interpolated Surface object produced by RMT;\n", - "- the fitted RBF value on each surface vertex, useful for checking extraction; and\n", - "- the signed vertical error relative to the known synthetic horizon, useful only because this POC has ground truth.\n", - "\n", - "In Geoscience ANALYST, colour the observation points by **Implicit scalar value** and the surface by **True vertical error** to inspect the result." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "156ec95b0d1e1877", - "metadata": {}, - "outputs": [], - "source": [ - "output_directory = Path.cwd() / \"ferreus_output\"\n", - "output_directory.mkdir(exist_ok=True)\n", - "\n", - "folded_output_path = output_directory / (\n", - " f\"ferreus_rbf_folded_horizon_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", - ")\n", - "\n", - "with Workspace.create(folded_output_path) as workspace:\n", - " group = ContainerGroup.create(\n", - " workspace,\n", - " name=\"Ferreus folded horizon POC\",\n", - " )\n", - "\n", - " observations = Points.create(\n", - " workspace,\n", - " name=\"RBF scalar observations\",\n", - " vertices=source_points,\n", - " parent=group,\n", - " )\n", - " observations.add_data(\n", - " {\n", - " \"Implicit scalar value\": {\n", - " \"values\": source_values,\n", - " }\n", - " }\n", - " )\n", - "\n", - " folded_surface = Surface.create(\n", - " workspace,\n", - " name=\"RBF folded horizon - isovalue 0\",\n", - " vertices=folded_vertices,\n", - " cells=folded_cells,\n", - " parent=group,\n", - " )\n", - " folded_surface.add_data(\n", - " {\n", - " \"RBF value on surface\": {\n", - " \"values\": fitted_values_at_vertices,\n", - " },\n", - " \"True vertical error\": {\n", - " \"values\": true_vertical_errors,\n", - " },\n", - " }\n", - " )\n", - "\n", - "print(\"Created:\", folded_output_path)" - ] - }, - { - "cell_type": "markdown", - "id": "f74bf0dad647895f", - "metadata": {}, - "source": [ - "## What this baseline proves\n", - "\n", - "This section demonstrates the complete technical path: scattered scalar observations -> ferreus_rbf scalar field -> ferreus_rmt zero-isosurface -> geoh5 objects for ANALYST.\n", - "\n", - "It isn't yet a test of geological orientation or exclusion constraints. Every source point above was assigned a full scalar value from the known analytic field. Real geological inputs more commonly include:\n", - "\n", - "- **contact/interface points**, which say f(p) = 0;\n", - "- **orientation measurements**, which say that the gradient of f should point along a measured normal; and\n", - "- **exclusion/inequality observations**, which specify only that a point must be on one side of the surface, such as f(p) > 0 or f(p) < 0.\n", - "\n", - "The public ferreus_rbf Python constructor accepts point coordinates and scalar values; it does not expose separate derivative or inequality observations. The next POC sections can therefore test value-based encodings:\n", - "\n", - "1. approximate an orientation at contact point p with offset samples p + delta times n and p - delta times n, assigned positive and negative values respectively; and\n", - "2. approximate an exclusion constraint with a signed pseudo-value on the required side of the surface.\n", - "\n", - "Those are approximations that should be measured against this baseline." - ] - }, - { - "cell_type": "markdown", - "id": "23da01a0b8b33419", - "metadata": {}, - "source": [ - "## Compare the existing marching-cubes method with Ferreus RMT\n", - "\n", - "This is an extraction-algorithm comparison. Both algorithms receive the same fitted RBF field, model bounds, lattice spacing, and isovalue. That keeps the scalar model fixed so differences in the resulting meshes come from marching cubes versus RMT rather than from different geological interpolation methods.\n", - "\n", - "The existing application first samples data onto a complete regular grid and then calls scikit-image marching cubes. The code below performs that same dense-grid extraction. RMT instead evaluates the field while following the surface outward from seed points.\n", - "\n", - "This comparison does not yet reproduce every part of the existing application's preprocessing, such as weighted averaging, maximum distance, or topographic masking. A later end-to-end test can feed both implementations the same real geoh5 object." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9117b165c6dfa0e7", - "metadata": {}, - "outputs": [], - "source": [ - "from time import perf_counter\n", - "\n", - "from skimage.measure import marching_cubes\n", - "\n", - "\n", - "# Construct the same kind of dense regular grid used by the existing\n", - "# iso_surfaces application. The array axes are X, Y, then Z.\n", - "comparison_grid = [\n", - " np.arange(\n", - " MODEL_EXTENTS[axis_index],\n", - " MODEL_EXTENTS[axis_index + 3] + 0.5 * folded_resolution,\n", - " folded_resolution,\n", - " dtype=np.float64,\n", - " )\n", - " for axis_index in range(3)\n", - "]\n", - "\n", - "comparison_x, comparison_y, comparison_z = np.meshgrid(\n", - " comparison_grid[0],\n", - " comparison_grid[1],\n", - " comparison_grid[2],\n", - " indexing=\"ij\",\n", - ")\n", - "\n", - "# With indexing=\"ij\", array axes 0, 1, and 2 are X, Y, and Z.\n", - "# C-order ravel and reshape are therefore a matched pair. Using\n", - "# order=\"F\" here would reorder otherwise correctly paired values.\n", - "\n", - "expected_grid_shape = tuple(axis.size for axis in comparison_grid)\n", - "assert comparison_x.shape == expected_grid_shape\n", - "sample_index = tuple(min(2, axis.size - 1) for axis in comparison_grid)\n", - "assert np.allclose(\n", - " [\n", - " comparison_x[sample_index],\n", - " comparison_y[sample_index],\n", - " comparison_z[sample_index],\n", - " ],\n", - " [comparison_grid[axis][sample_index[axis]] for axis in range(3)],\n", - ")\n", - "comparison_grid_points = np.column_stack(\n", - " [\n", - " comparison_x.ravel(),\n", - " comparison_y.ravel(),\n", - " comparison_z.ravel(),\n", - " ]\n", - ")\n", - "\n", - "marching_cubes_total_start = perf_counter()\n", - "\n", - "grid_evaluation_start = perf_counter()\n", - "\n", - "comparison_grid_values = np.asarray(\n", - " rbf_model.evaluate_targets(comparison_grid_points),\n", - " dtype=np.float64,\n", - ").reshape(comparison_x.shape)\n", - "\n", - "grid_evaluation_seconds = perf_counter() - grid_evaluation_start\n", - "\n", - "marching_cubes_start = perf_counter()\n", - "mc_index_vertices, mc_cells, _mc_normals, _mc_values = marching_cubes(\n", - " comparison_grid_values,\n", - " level=folded_isovalue,\n", - ")\n", - "\n", - "# marching_cubes returns positions in array-index coordinates. Map each\n", - "# axis back to the original X, Y, and Z coordinates, matching the logic\n", - "# in surface_apps.iso_surfaces.utils.extract_iso_surfaces.\n", - "mc_vertices = np.column_stack(\n", - " [\n", - " np.interp(\n", - " mc_index_vertices[:, axis_index],\n", - " np.arange(comparison_grid[axis_index].size),\n", - " comparison_grid[axis_index],\n", - " )\n", - " for axis_index in range(3)\n", - " ]\n", - ").astype(np.float64)\n", - "mc_cells = np.asarray(mc_cells, dtype=np.uint32)\n", - "\n", - "marching_cubes_algorithm_seconds = perf_counter() - marching_cubes_start\n", - "marching_cubes_total_seconds = perf_counter() - marching_cubes_total_start\n", - "\n", - "print(\"Regular-grid shape:\", comparison_grid_values.shape)\n", - "print(\"Regular-grid sample count:\", comparison_grid_points.shape[0])\n", - "print(\"Marching-cubes vertices:\", mc_vertices.shape)\n", - "print(\"Marching-cubes triangles:\", mc_cells.shape)\n", - "print(\"Grid evaluation seconds:\", grid_evaluation_seconds)\n", - "print(\"Marching-cubes algorithm seconds:\", marching_cubes_algorithm_seconds)\n", - "print(\"Marching-cubes total seconds:\", marching_cubes_total_seconds)" - ] - }, - { - "cell_type": "markdown", - "id": "b479c5ef098030d6", - "metadata": {}, - "source": [ - "### Time RMT against the same field\n", - "\n", - "The RMT timing includes its field evaluations and mesh construction. The RBF evaluator is already cached for both methods. These single-run timings are useful for orientation, but they are not a rigorous performance benchmark; repeated warm runs and larger datasets would be needed for that." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c20c7765fc73c434", - "metadata": {}, - "outputs": [], - "source": [ - "rmt_total_start = perf_counter()\n", - "\n", - "comparison_rmt_mesh = build_isosurface(\n", - " seed_points=folded_seed_points,\n", - " extents=MODEL_EXTENTS,\n", - " resolution=folded_resolution,\n", - " isovalue=folded_isovalue,\n", - " surface_fn=rbf_model.evaluate_targets,\n", - " gradient_fn=rbf_model.evaluate_targets_with_gradients,\n", - " cluster_method=ClusterMethod.CurvatureWeighted,\n", - " boundary_closure=BoundaryClosure.None_,\n", - ")\n", - "\n", - "rmt_total_seconds = perf_counter() - rmt_total_start\n", - "\n", - "comparison_rmt_vertices = np.asarray(\n", - " comparison_rmt_mesh.vertices,\n", - " dtype=np.float64,\n", - ")\n", - "comparison_rmt_cells = np.asarray(\n", - " comparison_rmt_mesh.facets,\n", - " dtype=np.uint32,\n", - ")\n", - "\n", - "print(\"RMT vertices:\", comparison_rmt_vertices.shape)\n", - "print(\"RMT triangles:\", comparison_rmt_cells.shape)\n", - "print(\"RMT total seconds:\", rmt_total_seconds)" - ] - }, - { - "cell_type": "markdown", - "id": "22b19f3ec6d5d959", - "metadata": {}, - "source": [ - "### Compare geometry and topology\n", - "\n", - "The triangle-quality score below is 1.0 for an equilateral triangle and approaches zero for a very thin or degenerate triangle. The 5th percentile is particularly useful because it describes the poor-quality end of each mesh without depending on only one extreme triangle.\n", - "\n", - "An edge used by more than two triangles is non-manifold. Boundary edges are expected here because the folded horizon reaches the sides of the extraction box and boundary closure was disabled. Therefore, a non-zero boundary-edge count does not by itself indicate a broken mesh." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "412f81de43194b5a", - "metadata": {}, - "outputs": [], - "source": [ - "def calculate_mesh_statistics(\n", - " mesh_vertices: np.ndarray,\n", - " mesh_cells: np.ndarray,\n", - " elapsed_seconds: float,\n", - ") -> tuple[dict[str, float | int], np.ndarray, np.ndarray]:\n", - " \"\"\"Calculate comparable geometric and topological mesh diagnostics.\"\"\"\n", - " triangles = mesh_vertices[mesh_cells]\n", - "\n", - " edge_01 = triangles[:, 1] - triangles[:, 0]\n", - " edge_12 = triangles[:, 2] - triangles[:, 1]\n", - " edge_20 = triangles[:, 0] - triangles[:, 2]\n", - "\n", - " squared_edge_sum = (\n", - " np.sum(edge_01**2, axis=1)\n", - " + np.sum(edge_12**2, axis=1)\n", - " + np.sum(edge_20**2, axis=1)\n", - " )\n", - " twice_area = np.linalg.norm(\n", - " np.cross(edge_01, -edge_20),\n", - " axis=1,\n", - " )\n", - " triangle_quality = np.divide(\n", - " 2.0 * np.sqrt(3.0) * twice_area,\n", - " squared_edge_sum,\n", - " out=np.zeros_like(twice_area),\n", - " where=squared_edge_sum > 0.0,\n", - " )\n", - "\n", - " mesh_edges = np.vstack(\n", - " [\n", - " mesh_cells[:, [0, 1]],\n", - " mesh_cells[:, [1, 2]],\n", - " mesh_cells[:, [2, 0]],\n", - " ]\n", - " )\n", - " mesh_edges.sort(axis=1)\n", - " _unique_edges, edge_use_counts = np.unique(\n", - " mesh_edges,\n", - " axis=0,\n", - " return_counts=True,\n", - " )\n", - "\n", - " fitted_mesh_values = np.asarray(\n", - " rbf_model.evaluate(mesh_vertices),\n", - " dtype=np.float64,\n", - " ).reshape(-1)\n", - " mesh_true_vertical_errors = true_geological_field(mesh_vertices)\n", - "\n", - " statistics = {\n", - " \"vertices\": int(mesh_vertices.shape[0]),\n", - " \"triangles\": int(mesh_cells.shape[0]),\n", - " \"elapsed_seconds\": float(elapsed_seconds),\n", - " \"quality_mean\": float(np.mean(triangle_quality)),\n", - " \"quality_5th_percentile\": float(np.percentile(triangle_quality, 5)),\n", - " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", - " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", - " \"mean_absolute_isovalue_error\": float(\n", - " np.mean(np.abs(fitted_mesh_values - folded_isovalue))\n", - " ),\n", - " \"mean_absolute_vertical_error\": float(\n", - " np.mean(np.abs(mesh_true_vertical_errors))\n", - " ),\n", - " \"maximum_absolute_vertical_error\": float(\n", - " np.max(np.abs(mesh_true_vertical_errors))\n", - " ),\n", - " }\n", - "\n", - " return statistics, fitted_mesh_values, mesh_true_vertical_errors\n", - "\n", - "\n", - "mc_statistics, mc_fitted_values, mc_true_vertical_errors = calculate_mesh_statistics(\n", - " mc_vertices,\n", - " mc_cells,\n", - " marching_cubes_total_seconds,\n", - ")\n", - "rmt_statistics, rmt_fitted_values, rmt_true_vertical_errors = calculate_mesh_statistics(\n", - " comparison_rmt_vertices,\n", - " comparison_rmt_cells,\n", - " rmt_total_seconds,\n", - ")\n", - "\n", - "comparison_rows = [\n", - " (\"Vertices\", \"vertices\", \".0f\"),\n", - " (\"Triangles\", \"triangles\", \".0f\"),\n", - " (\"Total time (seconds)\", \"elapsed_seconds\", \".4f\"),\n", - " (\"Mean triangle quality\", \"quality_mean\", \".4f\"),\n", - " (\"5th-percentile quality\", \"quality_5th_percentile\", \".4f\"),\n", - " (\"Boundary edges\", \"boundary_edges\", \".0f\"),\n", - " (\"Non-manifold edges\", \"non_manifold_edges\", \".0f\"),\n", - " (\"Mean absolute isovalue error\", \"mean_absolute_isovalue_error\", \".4f\"),\n", - " (\"Mean absolute vertical error\", \"mean_absolute_vertical_error\", \".4f\"),\n", - " (\"Maximum absolute vertical error\", \"maximum_absolute_vertical_error\", \".4f\"),\n", - "]\n", - "\n", - "print(f\"{'Metric':38s} {'Marching cubes':>18s} {'Ferreus RMT':>18s}\")\n", - "print(\"-\" * 76)\n", - "for label, key, value_format in comparison_rows:\n", - " mc_text = format(mc_statistics[key], value_format)\n", - " rmt_text = format(rmt_statistics[key], value_format)\n", - " print(f\"{label:38s} {mc_text:>18s} {rmt_text:>18s}\")" - ] - }, - { - "cell_type": "markdown", - "id": "d85ca22ca792b883", - "metadata": {}, - "source": [ - "### View the two extraction results\n", - "\n", - "The surfaces should have the same large-scale folded geometry because they represent the same zero-isovalue of the same RBF. Differences to look for are triangulation density, skinny triangles, boundary shape, and small positional deviations introduced by RMT's regularisation and vertex clustering, though those things are difficult to discern when there are so many triangles." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c46d0141b8471c12", - "metadata": {}, - "outputs": [], - "source": [ - "figure = plt.figure(figsize=(14, 6))\n", - "\n", - "comparison_meshes = [\n", - " (\"Existing method: marching cubes\", mc_vertices, mc_cells, \"darkorange\"),\n", - " (\"Ferreus RMT\", comparison_rmt_vertices, comparison_rmt_cells, \"steelblue\"),\n", - "]\n", - "\n", - "for plot_index, (title, plot_vertices, plot_cells, colour) in enumerate(\n", - " comparison_meshes,\n", - " start=1,\n", - "):\n", - " axis = figure.add_subplot(1, 2, plot_index, projection=\"3d\")\n", - " axis.plot_trisurf(\n", - " plot_vertices[:, 0],\n", - " plot_vertices[:, 1],\n", - " plot_vertices[:, 2],\n", - " triangles=plot_cells,\n", - " color=colour,\n", - " linewidth=0.05,\n", - " alpha=0.85,\n", - " )\n", - " axis.scatter(\n", - " folded_seed_points[:, 0],\n", - " folded_seed_points[:, 1],\n", - " folded_seed_points[:, 2],\n", - " color=\"black\",\n", - " s=8,\n", - " label=\"RMT seeds\",\n", - " )\n", - " axis.set_xlim(MODEL_EXTENTS[[0, 3]])\n", - " axis.set_ylim(MODEL_EXTENTS[[1, 4]])\n", - " axis.set_zlim(MODEL_EXTENTS[[2, 5]])\n", - " axis.set_box_aspect(MODEL_EXTENTS[3:] - MODEL_EXTENTS[:3])\n", - " axis.set_xlabel(\"X\")\n", - " axis.set_ylabel(\"Y\")\n", - " axis.set_zlabel(\"Z\")\n", - " axis.view_init(elev=25, azim=-55)\n", - " axis.set_title(title)\n", - "\n", - "figure.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "2852865e68b17840", - "metadata": {}, - "source": [ - "### Export both comparison surfaces to geoh5\n", - "\n", - "Both meshes are placed in the same group so they can be toggled, coloured, and compared in Geoscience ANALYST. Each receives the RBF value and known vertical error at its vertices." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ccfd27e9a9805429", - "metadata": {}, - "outputs": [], - "source": [ - "comparison_output_path = output_directory / (\n", - " f\"marching_cubes_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", - ")\n", - "\n", - "with Workspace.create(comparison_output_path) as workspace:\n", - " comparison_group = ContainerGroup.create(\n", - " workspace,\n", - " name=\"Marching cubes versus Ferreus RMT\",\n", - " )\n", - "\n", - " mc_surface = Surface.create(\n", - " workspace,\n", - " name=\"Existing method - marching cubes\",\n", - " vertices=mc_vertices,\n", - " cells=mc_cells,\n", - " parent=comparison_group,\n", - " )\n", - " mc_surface.add_data(\n", - " {\n", - " \"RBF value on surface\": {\"values\": mc_fitted_values},\n", - " \"True vertical error\": {\"values\": mc_true_vertical_errors},\n", - " }\n", - " )\n", - "\n", - " rmt_surface = Surface.create(\n", - " workspace,\n", - " name=\"Ferreus RMT\",\n", - " vertices=comparison_rmt_vertices,\n", - " cells=comparison_rmt_cells,\n", - " parent=comparison_group,\n", - " )\n", - " rmt_surface.add_data(\n", - " {\n", - " \"RBF value on surface\": {\"values\": rmt_fitted_values},\n", - " \"True vertical error\": {\"values\": rmt_true_vertical_errors},\n", - " }\n", - " )\n", - "\n", - "print(\"Created:\", comparison_output_path)" - ] - }, - { - "cell_type": "markdown", - "id": "e66498617a269a8f", - "metadata": {}, - "source": [ - "### What this comparison does and does not establish\n", - "\n", - "This establishes whether RMT is a credible alternative mesh extractor for a field that already exists. It can compare mesh density, triangle quality, topology, approximate accuracy, and timing under controlled conditions.\n", - "\n", - "It does not establish support for geological orientation or exclusion constraints. Those affect how the scalar field is fitted from observations, which happens before either marching cubes or RMT is called." - ] - }, - { - "cell_type": "markdown", - "id": "f91e0a6edbb84baa", - "metadata": {}, - "source": [ - "## Compare Ferreus and SciPy RBF interpolation\n", - "\n", - "The mesh extractor can only be as good as the scalar field it receives. This section therefore compares the Ferreus and SciPy interpolators before doing any isosurfacing. Both models receive exactly the same source observations and validation points.\n", - "\n", - "A fair comparison requires equivalent model choices. Ferreus uses its linear radial kernel, a linear drift, and zero nugget. SciPy uses its linear kernel, a degree-one polynomial, and zero smoothing. These settings are the closest direct match; changing only one model's kernel, polynomial term, or smoothing would compare different mathematical models rather than implementations.\n", - "\n", - "The two source counts below show a small problem and a moderately larger one. Timings are single runs, so treat them as directional rather than as a formal performance benchmark." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "edaa15ec9a814d49", - "metadata": {}, - "outputs": [], - "source": [ - "from time import perf_counter\n", - "\n", - "from scipy.interpolate import RBFInterpolator as ScipyRBFInterpolator\n", - "\n", - "\n", - "# Compare a small interpolation problem with a moderately larger one\n", - "RBF_BENCHMARK_SIZES = (400, 4431)\n", - "RBF_BENCHMARK_VALIDATION_POINTS = 2000\n", - "\n", - "# Fixed seed for reproducibility\n", - "benchmark_rng = np.random.default_rng(1234567)\n", - "\n", - "# Points used to measure accuracy between observations\n", - "benchmark_validation_points = benchmark_rng.uniform(\n", - " low=MODEL_EXTENTS[:3],\n", - " high=MODEL_EXTENTS[3:],\n", - " size=(RBF_BENCHMARK_VALIDATION_POINTS, 3),\n", - ").astype(np.float64)\n", - "benchmark_validation_truth = true_geological_field(benchmark_validation_points)\n", - "\n", - "# Ferreus prepares a cached evaluator before RMT repeatedly requests field\n", - "# values. Padding includes the halo RMT may sample outside the model box (prevents an RMT query from falling outside the\n", - "# cached evaluator's allowable domain)\n", - "benchmark_evaluator_extents = MODEL_EXTENTS + np.array(\n", - " [-25.0, -25.0, -25.0, 25.0, 25.0, 25.0],\n", - " dtype=np.float64,\n", - ")\n", - "benchmark_rows = []\n", - "\n", - "for benchmark_size in RBF_BENCHMARK_SIZES:\n", - " # Generate one source data set to give the exact same observations to\n", - " # both libraries\n", - " benchmark_source_points = benchmark_rng.uniform(\n", - " low=MODEL_EXTENTS[:3],\n", - " high=MODEL_EXTENTS[3:],\n", - " size=(benchmark_size, 3),\n", - " ).astype(np.float64)\n", - " benchmark_source_values = true_geological_field(benchmark_source_points).astype(\n", - " np.float64\n", - " )\n", - "\n", - " # Time construction and solution of the Ferreus model. linear kernel,\n", - " # linear drift, and zero nugget are matched to SciPy below\n", - " ferreus_start = perf_counter()\n", - " ferreus_benchmark_model = RBFInterpolator(\n", - " points=benchmark_source_points,\n", - " values=benchmark_source_values,\n", - " interpolant_settings=InterpolantSettings(\n", - " kernel_type=RBFKernelType.Linear,\n", - " drift=Drift.Linear,\n", - " nugget=0.0,\n", - " ),\n", - " )\n", - " ferreus_fit_seconds = perf_counter() - ferreus_start\n", - "\n", - " # Report evaluator preparation separately.\n", - " # RMT can use the cached evaluate_targets path multiple times.\n", - " ferreus_start = perf_counter()\n", - " ferreus_benchmark_model.build_evaluator(benchmark_evaluator_extents)\n", - " ferreus_setup_seconds = perf_counter() - ferreus_start\n", - "\n", - " # Evaluate unseen locations using the cached path used by RMT.\n", - " ferreus_start = perf_counter()\n", - " ferreus_validation_values = np.asarray(\n", - " ferreus_benchmark_model.evaluate_targets(benchmark_validation_points),\n", - " dtype=np.float64,\n", - " ).reshape(-1)\n", - " ferreus_evaluate_seconds = perf_counter() - ferreus_start\n", - "\n", - " # Re-evaluate source points to test how closely the zero-nugget model\n", - " # reproduces the observations used to fit it\n", - " ferreus_source_values = np.asarray(\n", - " ferreus_benchmark_model.evaluate_at_source(add_nugget=False),\n", - " dtype=np.float64,\n", - " ).reshape(-1)\n", - "\n", - " # These are SciPy's closest equivalent settings. Matching the kernel,\n", - " # polynomial trend, and smoothing should help keep this an implementation comparison\n", - " # rather than a comparison of two different interpolation models.\n", - " scipy_start = perf_counter()\n", - " scipy_benchmark_model = ScipyRBFInterpolator(\n", - " benchmark_source_points,\n", - " benchmark_source_values,\n", - " kernel=\"linear\",\n", - " degree=1,\n", - " smoothing=0.0,\n", - " )\n", - " scipy_fit_seconds = perf_counter() - scipy_start\n", - "\n", - " # Evaluate exactly the same unseen validation locations with SciPy.\n", - " scipy_start = perf_counter()\n", - " scipy_validation_values = np.asarray(\n", - " scipy_benchmark_model(benchmark_validation_points),\n", - " dtype=np.float64,\n", - " ).reshape(-1)\n", - " scipy_evaluate_seconds = perf_counter() - scipy_start\n", - "\n", - " # This source-point call is kept outside the validation timing above.\n", - " scipy_source_values = np.asarray(\n", - " scipy_benchmark_model(benchmark_source_points),\n", - " dtype=np.float64,\n", - " ).reshape(-1)\n", - "\n", - " # Store three different accuracy checks:\n", - " # source RMSE measures reproduction of fitted observations;\n", - " # validation RMSE measures accuracy away from observations;\n", - " # between-model RMSE measures direct Ferreus/SciPy agreement.\n", - " benchmark_rows.append(\n", - " {\n", - " \"sources\": benchmark_size,\n", - " \"ferreus_fit_s\": ferreus_fit_seconds,\n", - " \"ferreus_setup_s\": ferreus_setup_seconds,\n", - " \"ferreus_eval_s\": ferreus_evaluate_seconds,\n", - " \"scipy_fit_s\": scipy_fit_seconds,\n", - " \"scipy_eval_s\": scipy_evaluate_seconds,\n", - " \"ferreus_source_rmse\": float(\n", - " np.sqrt(np.mean((ferreus_source_values - benchmark_source_values) ** 2))\n", - " ),\n", - " \"scipy_source_rmse\": float(\n", - " np.sqrt(np.mean((scipy_source_values - benchmark_source_values) ** 2))\n", - " ),\n", - " \"ferreus_validation_rmse\": float(\n", - " np.sqrt(\n", - " np.mean(\n", - " (ferreus_validation_values - benchmark_validation_truth) ** 2\n", - " )\n", - " )\n", - " ),\n", - " \"scipy_validation_rmse\": float(\n", - " np.sqrt(\n", - " np.mean((scipy_validation_values - benchmark_validation_truth) ** 2)\n", - " )\n", - " ),\n", - " \"between_models_rmse\": float(\n", - " np.sqrt(\n", - " np.mean((ferreus_validation_values - scipy_validation_values) ** 2)\n", - " )\n", - " ),\n", - " }\n", - " )\n", - "\n", - "# F and S = Ferreus and SciPy. Timing columns are elapsed seconds;\n", - "# the remaining columns are scalar-field RMSE values.\n", - "print(\n", - " f\"{'N':>6} {'F fit':>9} {'F setup':>9} {'F eval':>9} \"\n", - " f\"{'S fit':>9} {'S eval':>9} {'F val RMSE':>12} \"\n", - " f\"{'S val RMSE':>12} {'F-S RMSE':>12}\"\n", - ")\n", - "for row in benchmark_rows:\n", - " print(\n", - " f\"{row['sources']:6d} \"\n", - " f\"{row['ferreus_fit_s']:9.4f} \"\n", - " f\"{row['ferreus_setup_s']:9.4f} \"\n", - " f\"{row['ferreus_eval_s']:9.4f} \"\n", - " f\"{row['scipy_fit_s']:9.4f} \"\n", - " f\"{row['scipy_eval_s']:9.4f} \"\n", - " f\"{row['ferreus_validation_rmse']:12.6g} \"\n", - " f\"{row['scipy_validation_rmse']:12.6g} \"\n", - " f\"{row['between_models_rmse']:12.6g}\"\n", - " )\n", - "\n", - "# Source errors answer whether the unsmoothed models pass through their\n", - "# observations, which is different from prediction accuracy between them.\n", - "print(\"\\nSource-point RMSE (zero means exact interpolation):\")\n", - "for row in benchmark_rows:\n", - " print(\n", - " f\"N={row['sources']}: \"\n", - " f\"Ferreus={row['ferreus_source_rmse']:.6g}, \"\n", - " f\"SciPy={row['scipy_source_rmse']:.6g}\"\n", - " )" - ] - }, - { - "cell_type": "markdown", - "id": "3f25676179f04fc1", - "metadata": {}, - "source": [ - "### How to interpret the RBF comparison\n", - "\n", - "- Source-point RMSE checks whether each zero-smoothing model reproduces the observations used to fit it.\n", - "- Validation RMSE compares each model with the known synthetic field at locations that were not used for fitting. This measures interpolation accuracy for this controlled example.\n", - "- F-S RMSE compares the two implementations directly. If it is tiny relative to the validation error, their matched settings are producing essentially the same scalar field.\n", - "- Ferreus setup is reported separately because build_evaluator prepares its cached evaluation path. RMT can then call evaluate_targets repeatedly. Comparing only this cached evaluation time with SciPy's evaluation time would omit that setup cost.\n", - "- These results test ordinary value interpolation only. They do not demonstrate orientation, contact, or inequality/exclusion-constraint support.\n", - "\n", - "SciPy documentation notes that a global RBF solve has memory requirements that grow quadratically with the number of source points. Its neighbors option can reduce that cost by using a local subset, but doing so changes the interpolation method and would no longer be a direct comparison with the global Ferreus model used here.\n", - "\n", - "Overall, the two options seem to produce similar interpolations. The interpolators reproduce their input observations with little difference between them. At a larger input size, ferreus has slightly weaker precision but both are good. Ferreus seems to be slightly faster to run, with better evaluation performance." - ] - }, - { - "cell_type": "markdown", - "id": "25c1e859ea5064ee", - "metadata": {}, - "source": [ - "## Constraints\n", - "\n", - "The following is my understanding of how different constraints can be handled, whether they are readily available, etc.\n", - "\n", - "### Current conclusion\n", - "\n", - "The public ferreus_rbf Python fitting API accepts coordinates and scalar values. It can calculate gradients from a fitted field, but it cannot accept measured gradients as fitting observations. It also does not expose an inequality-constrained solver. Therefore, contacts are directly supported, while geological orientation and inequality constraints (used for exclusion) are not directly supported. This limitation belongs to the field-fitting stage; ferreus_rmt can still extract a mesh from any scalar field supplied to it.\n", - "\n", - "### Workarounds which may be worth testing\n", - "\n", - "1. Represent an above/below inequality with a small positive or negative pseudo-value. This is only an approximation because it gives the point an exact value when only its sign is really known.\n", - "2. At a contact with a measured normal vector, add two offset points on opposite sides of the contact and assign them small negative and positive values. This encourages the fitted field to have the requested orientation. It is less reliable for orientations measured away from a known contact because the correct scalar level there is unknown.\n", - "\n", - "Both workarounds introduce choices such as the pseudo-value and offset distance. A result is only convincing if it remains reasonably stable when those choices are varied. They will be approximations.\n", - "\n", - "### Potential first steps to investigate\n", - "\n", - "Start with the existing synthetic folded horizon, where the correct answer is known, and compare: contacts only; contacts plus pseudo-inequalities; contacts plus orientation offset points; and all approximate constraints together. For each case, check contact RMSE, the percentage of above/below conditions satisfied, the angle between requested and predicted normals, and error against the known horizon. Export the resulting RMT surfaces to geoh5 for visual comparison.\n", - "\n", - "The workaround would be promising if it keeps contacts accurate, satisfies most inequalities, improves orientation agreement, and produces similar surfaces across reasonable parameter choices. My ballpark estimate of going down this path is somewhere around 1-2 days, depending on how well things work, what quirks or complications come up, things that I have to spend more time working through and reading about so I can properly understand them, etc." - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} +version https://git-lfs.github.com/spec/v1 +oid sha256:d53a061318b0fd704d6b4d654cf4f41b96d9eba63ad08ea62f7cfa9851d8728b +size 56551 diff --git a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb index 5752b56..48b6a00 100644 --- a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb +++ b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb @@ -1,700 +1,3 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "8cbb863ab5d141a7", - "metadata": {}, - "source": [ - "# Existing iso-surfaces versus Ferreus RMT on real geochemistry data\n", - "\n", - "This notebook compares the existing surface_apps.iso_surfaces marching-cubes extraction with ferreus_rmt using the geochem_pts dataset.\n", - "\n", - "The comparison deliberately holds the interpolated scalar field constant:\n", - "\n", - "1. the existing application interpolates Cu_ppm observations onto a regular grid;\n", - "2. marching cubes extracts a surface from that grid; and\n", - "3. Ferreus RMT extracts the same isovalue from a trilinear function over that same grid.\n", - "\n", - "Consequently, differences measured here are primarily differences between the mesh extractors." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d0d17fd27a2945db", - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import datetime\n", - "from pathlib import Path\n", - "from time import perf_counter\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "from ferreus_rmt import BoundaryClosure, ClusterMethod, build_isosurface\n", - "from geoh5py.groups import ContainerGroup\n", - "from geoh5py.objects import Points, Surface\n", - "from geoh5py.workspace import Workspace\n", - "from scipy.interpolate import RegularGridInterpolator\n", - "from scipy.sparse import coo_matrix\n", - "from scipy.sparse.csgraph import connected_components\n", - "from scipy.spatial import cKDTree\n", - "\n", - "from surface_apps.iso_surfaces.utils import (\n", - " entity_to_grid,\n", - " extract_iso_surfaces,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "d2863f57f31b4888", - "metadata": {}, - "source": [ - "## Configuration\n", - "\n", - "Cu_ppm = 200 was selected as a useful first anomaly boundary: it is above most observations but produces enough geometry for a meaningful comparison. These are experiment parameters, not universally correct geological thresholds.\n", - "\n", - "The 50 m grid resolution and 500 m maximum interpolation distance match the existing application's UI defaults. A grid point farther than the maximum distance from data remains unsupported (NaN)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "77ded4c6f1424af2", - "metadata": {}, - "outputs": [], - "source": [ - "input_candidates = [\n", - " Path.cwd() / \"geochem_pts\" / \"geochem_pts_example.geoh5\",\n", - " Path.cwd()\n", - " / \"surface_apps\"\n", - " / \"ferreus\"\n", - " / \"geochem_pts\"\n", - " / \"geochem_pts_example.geoh5\",\n", - "]\n", - "\n", - "INPUT_PATH = next(\n", - " (candidate for candidate in input_candidates if candidate.exists()),\n", - " None,\n", - ")\n", - "if INPUT_PATH is None:\n", - " raise FileNotFoundError(\n", - " \"Could not find geochem_pts/geochem_pts_example.geoh5 \"\n", - " \"from the current working directory.\"\n", - " )\n", - "\n", - "OBJECT_NAME = \"geochem_pts\"\n", - "DATA_NAME = \"Cu_ppm\"\n", - "CONTOUR_LEVEL = 200.0\n", - "RESOLUTION = 50.0\n", - "MAX_DISTANCE = 500.0\n", - "RMT_RESOLUTION = RESOLUTION\n", - "RMT_CLUSTER_METHOD = ClusterMethod.CurvatureWeighted\n", - "\n", - "OUTPUT_DIRECTORY = INPUT_PATH.parent / \"ferreus_output\"\n", - "OUTPUT_DIRECTORY.mkdir(exist_ok=True)\n", - "\n", - "print(\"Input:\", INPUT_PATH)\n", - "print(\"Field and isovalue:\", DATA_NAME, CONTOUR_LEVEL)" - ] - }, - { - "cell_type": "markdown", - "id": "be709f46aca747c0", - "metadata": {}, - "source": [ - "## Load the observations and create the shared field\n", - "\n", - "entity_to_grid is the existing application's preprocessing function. For a points object, it uses distance-weighted averaging to estimate the geochemical value at regular-grid locations.\n", - "\n", - "This conversion is shared by both methods below." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90a877f44c8b44fb", - "metadata": {}, - "outputs": [], - "source": [ - "grid_start = perf_counter()\n", - "\n", - "with Workspace(INPUT_PATH) as input_workspace:\n", - " source_object = input_workspace.get_entity(OBJECT_NAME)[0]\n", - " source_data = input_workspace.get_entity(DATA_NAME)[0]\n", - "\n", - " if source_object is None or source_data is None:\n", - " raise ValueError(f\"Could not find {OBJECT_NAME!r} and {DATA_NAME!r}.\")\n", - "\n", - " source_locations = np.asarray(\n", - " source_object.locations,\n", - " dtype=np.float64,\n", - " ).copy()\n", - " source_values = np.asarray(\n", - " source_data.values,\n", - " dtype=np.float64,\n", - " ).copy()\n", - "\n", - " comparison_grid, comparison_grid_values = entity_to_grid(\n", - " source_object,\n", - " source_data,\n", - " resolution=RESOLUTION,\n", - " max_distance=MAX_DISTANCE,\n", - " horizon=None,\n", - " )\n", - "\n", - "grid_seconds = perf_counter() - grid_start\n", - "\n", - "finite_source_values = source_values[np.isfinite(source_values)]\n", - "finite_grid = np.isfinite(comparison_grid_values)\n", - "\n", - "print(\"Observation count:\", source_locations.shape[0])\n", - "print(\"XYZ span:\", np.ptp(source_locations, axis=0))\n", - "print(\"Finite Cu observations:\", finite_source_values.size)\n", - "print(\n", - " \"Cu observation percentiles:\",\n", - " np.percentile(\n", - " finite_source_values,\n", - " [10, 25, 50, 75, 90, 95, 99],\n", - " ),\n", - ")\n", - "print(\"Grid shape:\", comparison_grid_values.shape)\n", - "print(\n", - " \"Supported grid nodes:\",\n", - " int(np.count_nonzero(finite_grid)),\n", - " \"/\",\n", - " comparison_grid_values.size,\n", - ")\n", - "print(\"Shared gridding time (seconds):\", grid_seconds)" - ] - }, - { - "cell_type": "markdown", - "id": "0a47f6b27d024ff5", - "metadata": {}, - "source": [ - "## Existing method: marching cubes\n", - "\n", - "This calls the existing application's extraction helper directly. Its result is the baseline users currently receive after the gridding step." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a953d300954d4afa", - "metadata": {}, - "outputs": [], - "source": [ - "marching_cubes_start = perf_counter()\n", - "\n", - "existing_surfaces = extract_iso_surfaces(\n", - " source_object,\n", - " comparison_grid,\n", - " [CONTOUR_LEVEL],\n", - " comparison_grid_values,\n", - ")\n", - "\n", - "marching_cubes_seconds = perf_counter() - marching_cubes_start\n", - "\n", - "mc_vertices = np.asarray(existing_surfaces[0][0], dtype=np.float64)\n", - "mc_cells = np.asarray(existing_surfaces[0][1], dtype=np.uint32)\n", - "\n", - "print(\"Marching-cubes vertices:\", mc_vertices.shape)\n", - "print(\"Marching-cubes triangles:\", mc_cells.shape)\n", - "print(\"Marching-cubes extraction time (seconds):\", marching_cubes_seconds)" - ] - }, - { - "cell_type": "markdown", - "id": "db6dad2bde0944d1", - "metadata": {}, - "source": [ - "## Ferreus RMT on the same regular-grid field\n", - "\n", - "RMT expects a callable scalar function instead of a dense value array. RegularGridInterpolator turns the existing grid into that function using trilinear interpolation. Unsupported regions remain NaN.\n", - "\n", - "As a surface-following method, RMT also requires seed points near each surface component. The helper below independently finds regular-grid cells whose eight finite corner values straddle the requested isovalue, then uses their centres as seeds. It does not use the marching-cubes result. This seed search scans the dense grid." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "da81ec4af37d4cd2", - "metadata": {}, - "outputs": [], - "source": [ - "grid_interpolator = RegularGridInterpolator(\n", - " tuple(comparison_grid),\n", - " comparison_grid_values,\n", - " method=\"linear\",\n", - " bounds_error=False,\n", - " fill_value=np.nan,\n", - ")\n", - "\n", - "\n", - "def shared_grid_field(targets: np.ndarray) -> np.ndarray:\n", - " \"\"\"Evaluate the same trilinear scalar field represented by the grid.\"\"\"\n", - " return np.asarray(grid_interpolator(targets), dtype=np.float64)\n", - "\n", - "\n", - "def find_crossing_cell_seeds(\n", - " grid: list[np.ndarray],\n", - " grid_values: np.ndarray,\n", - " isovalue: float,\n", - ") -> np.ndarray:\n", - " \"\"\"Return centres of fully supported grid cells crossing an isovalue.\"\"\"\n", - " corner_values = np.stack(\n", - " [\n", - " grid_values[x_slice, y_slice, z_slice]\n", - " for x_slice in (slice(None, -1), slice(1, None))\n", - " for y_slice in (slice(None, -1), slice(1, None))\n", - " for z_slice in (slice(None, -1), slice(1, None))\n", - " ],\n", - " axis=-1,\n", - " )\n", - "\n", - " fully_supported = np.all(np.isfinite(corner_values), axis=-1)\n", - " crosses_isovalue = (\n", - " fully_supported\n", - " & (np.min(corner_values, axis=-1) <= isovalue)\n", - " & (np.max(corner_values, axis=-1) >= isovalue)\n", - " )\n", - "\n", - " crossing_indices = np.argwhere(crosses_isovalue)\n", - " if crossing_indices.size == 0:\n", - " raise ValueError(f\"No fully supported grid cells cross {isovalue}.\")\n", - "\n", - " return np.column_stack(\n", - " [\n", - " 0.5\n", - " * (\n", - " grid[axis][crossing_indices[:, axis]]\n", - " + grid[axis][crossing_indices[:, axis] + 1]\n", - " )\n", - " for axis in range(3)\n", - " ]\n", - " ).astype(np.float64)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "aa3c77a2469f4b56", - "metadata": {}, - "outputs": [], - "source": [ - "rmt_total_start = perf_counter()\n", - "\n", - "rmt_seed_points = find_crossing_cell_seeds(\n", - " comparison_grid,\n", - " comparison_grid_values,\n", - " CONTOUR_LEVEL,\n", - ")\n", - "seed_detection_seconds = perf_counter() - rmt_total_start\n", - "\n", - "rmt_extents = np.array(\n", - " [\n", - " comparison_grid[0][0],\n", - " comparison_grid[1][0],\n", - " comparison_grid[2][0],\n", - " comparison_grid[0][-1],\n", - " comparison_grid[1][-1],\n", - " comparison_grid[2][-1],\n", - " ],\n", - " dtype=np.float64,\n", - ")\n", - "\n", - "rmt_extraction_start = perf_counter()\n", - "rmt_mesh = build_isosurface(\n", - " seed_points=rmt_seed_points,\n", - " extents=rmt_extents,\n", - " resolution=RMT_RESOLUTION,\n", - " isovalue=CONTOUR_LEVEL,\n", - " surface_fn=shared_grid_field,\n", - " cluster_method=RMT_CLUSTER_METHOD,\n", - " boundary_closure=BoundaryClosure.None_,\n", - ")\n", - "rmt_extraction_seconds = perf_counter() - rmt_extraction_start\n", - "rmt_total_seconds = perf_counter() - rmt_total_start\n", - "\n", - "rmt_vertices = np.asarray(rmt_mesh.vertices, dtype=np.float64)\n", - "rmt_cells = np.asarray(rmt_mesh.facets, dtype=np.uint32)\n", - "\n", - "print(\"Independent RMT seed points:\", rmt_seed_points.shape)\n", - "print(\"RMT vertices:\", rmt_vertices.shape)\n", - "print(\"RMT triangles:\", rmt_cells.shape)\n", - "print(\"Seed detection time (seconds):\", seed_detection_seconds)\n", - "print(\"RMT extraction time (seconds):\", rmt_extraction_seconds)\n", - "print(\"RMT seed plus extraction time (seconds):\", rmt_total_seconds)" - ] - }, - { - "cell_type": "markdown", - "id": "c55a39891ad44956", - "metadata": {}, - "source": [ - "## Compare mesh quality, topology, and field agreement\n", - "\n", - "Triangle quality is 1.0 for an equilateral triangle and approaches zero for increasingly skinny triangles. An edge shared by more than two triangles is non-manifold. Boundary edges can occur where a surface meets an unsupported region or the extraction bounds.\n", - "\n", - "There is no known true copper-anomaly surface in this dataset, so this notebook cannot report geological accuracy. Instead it measures how closely each mesh follows the shared interpolated 200 ppm field and how closely the two meshes agree spatially." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e14ca8c55c72457b", - "metadata": {}, - "outputs": [], - "source": [ - "def calculate_mesh_statistics(\n", - " mesh_vertices: np.ndarray,\n", - " mesh_cells: np.ndarray,\n", - " extraction_seconds: float,\n", - ") -> tuple[dict[str, float | int], np.ndarray]:\n", - " \"\"\"Calculate geometric, topological, and scalar-field diagnostics.\"\"\"\n", - " triangles = mesh_vertices[mesh_cells]\n", - " edge_01 = triangles[:, 1] - triangles[:, 0]\n", - " edge_12 = triangles[:, 2] - triangles[:, 1]\n", - " edge_20 = triangles[:, 0] - triangles[:, 2]\n", - "\n", - " twice_area = np.linalg.norm(np.cross(edge_01, -edge_20), axis=1)\n", - " squared_edge_sum = (\n", - " np.sum(edge_01**2, axis=1)\n", - " + np.sum(edge_12**2, axis=1)\n", - " + np.sum(edge_20**2, axis=1)\n", - " )\n", - " triangle_quality = np.divide(\n", - " 2.0 * np.sqrt(3.0) * twice_area,\n", - " squared_edge_sum,\n", - " out=np.zeros_like(twice_area),\n", - " where=squared_edge_sum > 0.0,\n", - " )\n", - "\n", - " mesh_edges = np.vstack(\n", - " [\n", - " mesh_cells[:, [0, 1]],\n", - " mesh_cells[:, [1, 2]],\n", - " mesh_cells[:, [2, 0]],\n", - " ]\n", - " )\n", - " mesh_edges.sort(axis=1)\n", - " unique_edges, edge_use_counts = np.unique(\n", - " mesh_edges,\n", - " axis=0,\n", - " return_counts=True,\n", - " )\n", - "\n", - " graph = coo_matrix(\n", - " (\n", - " np.ones(unique_edges.shape[0] * 2),\n", - " (\n", - " np.r_[unique_edges[:, 0], unique_edges[:, 1]],\n", - " np.r_[unique_edges[:, 1], unique_edges[:, 0]],\n", - " ),\n", - " ),\n", - " shape=(mesh_vertices.shape[0], mesh_vertices.shape[0]),\n", - " )\n", - " component_count = connected_components(\n", - " graph,\n", - " directed=False,\n", - " return_labels=False,\n", - " )\n", - "\n", - " field_values = shared_grid_field(mesh_vertices)\n", - " finite_field = np.isfinite(field_values)\n", - " isovalue_errors = np.abs(field_values[finite_field] - CONTOUR_LEVEL)\n", - "\n", - " statistics = {\n", - " \"vertices\": int(mesh_vertices.shape[0]),\n", - " \"triangles\": int(mesh_cells.shape[0]),\n", - " \"extraction_seconds\": float(extraction_seconds),\n", - " \"surface_area\": float(np.sum(0.5 * twice_area)),\n", - " \"quality_mean\": float(np.mean(triangle_quality)),\n", - " \"quality_5th_percentile\": float(np.percentile(triangle_quality, 5)),\n", - " \"components\": int(component_count),\n", - " \"boundary_edges\": int(np.count_nonzero(edge_use_counts == 1)),\n", - " \"non_manifold_edges\": int(np.count_nonzero(edge_use_counts > 2)),\n", - " \"finite_field_fraction\": float(np.mean(finite_field)),\n", - " \"mean_absolute_isovalue_error\": float(np.mean(isovalue_errors)),\n", - " \"maximum_absolute_isovalue_error\": float(np.max(isovalue_errors)),\n", - " }\n", - "\n", - " return statistics, field_values\n", - "\n", - "\n", - "mc_statistics, mc_field_values = calculate_mesh_statistics(\n", - " mc_vertices,\n", - " mc_cells,\n", - " marching_cubes_seconds,\n", - ")\n", - "rmt_statistics, rmt_field_values = calculate_mesh_statistics(\n", - " rmt_vertices,\n", - " rmt_cells,\n", - " rmt_extraction_seconds,\n", - ")\n", - "\n", - "comparison_rows = [\n", - " (\"Vertices\", \"vertices\", \".0f\"),\n", - " (\"Triangles\", \"triangles\", \".0f\"),\n", - " (\"Extraction seconds\", \"extraction_seconds\", \".4f\"),\n", - " (\"Surface area\", \"surface_area\", \".1f\"),\n", - " (\"Mean triangle quality\", \"quality_mean\", \".4f\"),\n", - " (\"5th-percentile triangle quality\", \"quality_5th_percentile\", \".4f\"),\n", - " (\"Connected components\", \"components\", \".0f\"),\n", - " (\"Boundary edges\", \"boundary_edges\", \".0f\"),\n", - " (\"Non-manifold edges\", \"non_manifold_edges\", \".0f\"),\n", - " (\"Finite field-value fraction\", \"finite_field_fraction\", \".4f\"),\n", - " (\"Mean absolute isovalue error\", \"mean_absolute_isovalue_error\", \".4f\"),\n", - " (\"Maximum absolute isovalue error\", \"maximum_absolute_isovalue_error\", \".4f\"),\n", - "]\n", - "\n", - "print(f\"{'Metric':40s} {'Marching cubes':>18s} {'Ferreus RMT':>18s}\")\n", - "print(\"-\" * 78)\n", - "for label, key, value_format in comparison_rows:\n", - " mc_text = format(mc_statistics[key], value_format)\n", - " rmt_text = format(rmt_statistics[key], value_format)\n", - " print(f\"{label:40s} {mc_text:>18s} {rmt_text:>18s}\")" - ] - }, - { - "cell_type": "markdown", - "id": "450f5356b5304251", - "metadata": {}, - "source": [ - "### Direct positional agreement\n", - "\n", - "Nearest-vertex distances provide an approximate comparison between the meshes. They are not exact point-to-triangle distances, so interpret them relative to the configured grid resolution." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a85ce67ad448473a", - "metadata": {}, - "outputs": [], - "source": [ - "mc_to_rmt_distances = cKDTree(rmt_vertices).query(mc_vertices)[0]\n", - "rmt_to_mc_distances = cKDTree(mc_vertices).query(rmt_vertices)[0]\n", - "symmetric_distances = np.r_[mc_to_rmt_distances, rmt_to_mc_distances]\n", - "\n", - "print(\"Marching cubes to RMT mean distance:\", mc_to_rmt_distances.mean())\n", - "print(\"RMT to marching cubes mean distance:\", rmt_to_mc_distances.mean())\n", - "print(\n", - " \"Symmetric distance percentiles (50%, 90%, 95%, 100%):\",\n", - " np.percentile(symmetric_distances, [50, 90, 95, 100]),\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "fa7ef208f1dd4dab", - "metadata": {}, - "source": [ - "### Reading the result\n", - "\n", - "Compare triangle quality, component counts, field error, and positional separation together. A method can produce attractive triangles while missing small components, or reproduce the isovalue accurately while retaining questionable interpolation artefacts. Because this dataset has no known true 200 ppm surface, neither output can be labelled geologically correct from these metrics alone." - ] - }, - { - "cell_type": "markdown", - "id": "2f3b645319f045de", - "metadata": {}, - "source": [ - "## Visual comparison\n", - "\n", - "The source points shown are finite observations at or above 200 ppm. They are not expected to lie exactly on the interpolated surface: the surface encloses/intersects regions where the continuous estimated field equals 200 ppm.\n", - "\n", - "Matplotlib's 3D renderer seems to have issues even when I tried to limit the complexity of what it's displaying, so I couldn't get this to work and it just keeps working for minutes." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "21305e00d59c4e30", - "metadata": {}, - "outputs": [], - "source": [ - "'''MAX_PREVIEW_TRIANGLES = 2_500\n", - "\n", - "\n", - "def select_preview_cells(\n", - " cells: np.ndarray,\n", - " maximum_triangles: int = MAX_PREVIEW_TRIANGLES,\n", - ") -> np.ndarray:\n", - " \"\"\"Select evenly spaced triangles for responsive notebook rendering.\"\"\"\n", - " if cells.shape[0] <= maximum_triangles:\n", - " return cells\n", - "\n", - " preview_indices = np.linspace(\n", - " 0,\n", - " cells.shape[0] - 1,\n", - " maximum_triangles,\n", - " dtype=int,\n", - " )\n", - " return cells[preview_indices]\n", - "\n", - "\n", - "high_copper = np.isfinite(source_values) & (source_values >= CONTOUR_LEVEL)\n", - "coordinate_span = np.ptp(source_locations, axis=0)\n", - "\n", - "figure = plt.figure(figsize=(15, 7))\n", - "plot_meshes = [\n", - " (\"Existing marching cubes\", mc_vertices, mc_cells, \"darkorange\"),\n", - " (\"Ferreus RMT\", rmt_vertices, rmt_cells, \"steelblue\"),\n", - "]\n", - "\n", - "for plot_index, (title, vertices, full_cells, colour) in enumerate(\n", - " plot_meshes,\n", - " start=1,\n", - "):\n", - " preview_cells = select_preview_cells(full_cells)\n", - " print(\n", - " f\"Rendering {title}: {preview_cells.shape[0]:,} of \"\n", - " f\"{full_cells.shape[0]:,} triangles\"\n", - " )\n", - "\n", - " axis = figure.add_subplot(1, 2, plot_index, projection=\"3d\")\n", - " axis.plot_trisurf(\n", - " vertices[:, 0],\n", - " vertices[:, 1],\n", - " vertices[:, 2],\n", - " triangles=preview_cells,\n", - " color=colour,\n", - " linewidth=0.0,\n", - " antialiased=False,\n", - " shade=False,\n", - " alpha=1.0,\n", - " )\n", - " axis.scatter(\n", - " source_locations[high_copper, 0],\n", - " source_locations[high_copper, 1],\n", - " source_locations[high_copper, 2],\n", - " color=\"black\",\n", - " s=5,\n", - " alpha=0.6,\n", - " label=f\"Observed {DATA_NAME} >= {CONTOUR_LEVEL:g}\",\n", - " )\n", - " axis.set_box_aspect(coordinate_span)\n", - " axis.set_xlabel(\"X\")\n", - " axis.set_ylabel(\"Y\")\n", - " axis.set_zlabel(\"Z\")\n", - " axis.view_init(elev=22, azim=-55)\n", - " axis.set_title(title)\n", - " axis.legend(loc=\"upper right\")\n", - "\n", - "figure.suptitle(f\"{DATA_NAME} = {CONTOUR_LEVEL:g} ppm\")\n", - "figure.subplots_adjust(left=0.02, right=0.98, bottom=0.02, top=0.90)\n", - "plt.show()'''" - ] - }, - { - "cell_type": "markdown", - "id": "8b950ca8bc474c13", - "metadata": {}, - "source": [ - "## Export both methods to geoh5\n", - "\n", - "The output includes a copy of the geochemistry points and both surfaces. Each mesh has its interpolated field value and nearest-vertex distance to the other method, making it possible to inspect differences in Geoscience ANALYST." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "62451e97bef04a06", - "metadata": {}, - "outputs": [], - "source": [ - "# Recreate the output directory if previous results were deleted after\n", - "# the configuration cell was run.\n", - "OUTPUT_DIRECTORY.mkdir(parents=True, exist_ok=True)\n", - "\n", - "output_path = OUTPUT_DIRECTORY / (\n", - " f\"real_geochem_mc_vs_rmt_{datetime.now():%Y%m%d_%H%M%S}.geoh5\"\n", - ")\n", - "\n", - "with Workspace.create(output_path) as output_workspace:\n", - " result_group = ContainerGroup.create(\n", - " output_workspace,\n", - " name=f\"{DATA_NAME} {CONTOUR_LEVEL:g} - MC versus RMT\",\n", - " )\n", - "\n", - " output_points = Points.create(\n", - " output_workspace,\n", - " name=OBJECT_NAME,\n", - " vertices=source_locations,\n", - " parent=result_group,\n", - " )\n", - " output_points.add_data({DATA_NAME: {\"values\": source_values}})\n", - "\n", - " output_mc_surface = Surface.create(\n", - " output_workspace,\n", - " name=f\"Existing marching cubes - {CONTOUR_LEVEL:g} ppm\",\n", - " vertices=mc_vertices,\n", - " cells=mc_cells,\n", - " parent=result_group,\n", - " )\n", - " output_mc_surface.add_data(\n", - " {\n", - " \"Interpolated Cu ppm\": {\"values\": mc_field_values},\n", - " \"Nearest RMT vertex distance\": {\"values\": mc_to_rmt_distances},\n", - " }\n", - " )\n", - "\n", - " output_rmt_surface = Surface.create(\n", - " output_workspace,\n", - " name=f\"Ferreus RMT - {CONTOUR_LEVEL:g} ppm\",\n", - " vertices=rmt_vertices,\n", - " cells=rmt_cells,\n", - " parent=result_group,\n", - " )\n", - " output_rmt_surface.add_data(\n", - " {\n", - " \"Interpolated Cu ppm\": {\"values\": rmt_field_values},\n", - " \"Nearest marching-cubes vertex distance\": {\"values\": rmt_to_mc_distances},\n", - " }\n", - " )\n", - "\n", - "print(\"Created:\", output_path)" - ] - }, - { - "cell_type": "markdown", - "id": "e642a55b31c948eb", - "metadata": {}, - "source": [ - "## Interpretation limits and next tests\n", - "\n", - "A successful run establishes that RMT can extract real, disconnected geochemical anomaly surfaces from the same field used by the current application, including in the presence of unsupported regions.\n", - "\n", - "It does not show that either surface is geological truth. Both inherit assumptions and artefacts from the shared weighted-average grid. It also does not test the RBF package or constraint fitting.\n", - "\n", - "Useful follow-up experiments are:\n", - "\n", - "- repeat at several contour levels and resolutions;\n", - "- compare results at approximately equal triangle count or positional error;\n", - "- test RMT boundary-closure modes on a genuinely enclosed component;\n", - "- replace grid-cell seed detection with seeds derived from selected observations or another production-safe strategy; and\n", - "- separately compare the existing weighted-average field with a ferreus_rbf field fitted from the raw geochemistry points." - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} +version https://git-lfs.github.com/spec/v1 +oid sha256:88f4ec321f6001b6892c2cdc7fa566c85c4ff2d4e2470cbb9fdceaeb75b5e610 +size 26859 From c7b4a5272eb26f2ff0975b841dfcf184942b45a0 Mon Sep 17 00:00:00 2001 From: domfournier Date: Fri, 21 Aug 2026 10:38:31 -0700 Subject: [PATCH 7/8] Remove duplicates --- surface_apps/ferreus/__init__.py | 8 -------- surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb | 3 --- .../ferreus/ferreus_rmt_real_geochem_comparison.ipynb | 3 --- .../ferreus/geochem_pts/geochem_pts_example.geoh5 | 3 --- 4 files changed, 17 deletions(-) delete mode 100644 surface_apps/ferreus/__init__.py delete mode 100644 surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb delete mode 100644 surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb delete mode 100644 surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 diff --git a/surface_apps/ferreus/__init__.py b/surface_apps/ferreus/__init__.py deleted file mode 100644 index 1b3743e..0000000 --- a/surface_apps/ferreus/__init__.py +++ /dev/null @@ -1,8 +0,0 @@ -# '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2026 Mira Geoscience Ltd. ' -# ' -# This file is part of surface-apps package. ' -# ' -# surface-apps is distributed under the terms and conditions of the MIT License -# (see LICENSE file at the root of this source code package). ' -# '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb deleted file mode 100644 index 86859aa..0000000 --- a/surface_apps/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:d53a061318b0fd704d6b4d654cf4f41b96d9eba63ad08ea62f7cfa9851d8728b -size 56551 diff --git a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb deleted file mode 100644 index 48b6a00..0000000 --- a/surface_apps/ferreus/ferreus_rmt_real_geochem_comparison.ipynb +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:88f4ec321f6001b6892c2cdc7fa566c85c4ff2d4e2470cbb9fdceaeb75b5e610 -size 26859 diff --git a/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 b/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 deleted file mode 100644 index 6d6d881..0000000 --- a/surface_apps/ferreus/geochem_pts/geochem_pts_example.geoh5 +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:b1410455e6510fe01bf6f4f32c380eb04aa9aba4dc8007fa6eecdb6bd1f576af -size 1376256 From 06d647a52542e31e924468ac3a69b0c14ef7a9cd Mon Sep 17 00:00:00 2001 From: domfournier Date: Fri, 21 Aug 2026 10:46:12 -0700 Subject: [PATCH 8/8] Update paths in notebook --- surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb | 4 ++-- .../ferreus/ferreus_rmt_real_geochem_comparison.ipynb | 4 ++-- surface_apps-assets/ferreus/geochem_pts_example.geoh5 | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb b/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb index 8073b52..19ee96d 100644 --- a/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb +++ b/surface_apps-assets/ferreus/ferreus_rmt_geoh5_poc.ipynb @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:91a597b583533d9228fa34d9b6114163a7c430583f42b38ac0102cd1315f6838 -size 947931 +oid sha256:4643cbf88711cda598ec36e191a2284466a3ee6fe5e07668253bd2ba383709fd +size 947791 diff --git a/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb b/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb index 8d00c73..5d0c1b4 100644 --- a/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb +++ b/surface_apps-assets/ferreus/ferreus_rmt_real_geochem_comparison.ipynb @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:0a732fc34d583d75ec4fbb4854e1fa1a269b1762c0cdf9c35b1317d2f303b697 -size 38221 +oid sha256:cbf1eeb37496c29e7f9c2c494f6905036da3ff77e4a0e4e0e6a7d7fd406a1224 +size 38170 diff --git a/surface_apps-assets/ferreus/geochem_pts_example.geoh5 b/surface_apps-assets/ferreus/geochem_pts_example.geoh5 index 88efc90..655a358 100644 --- a/surface_apps-assets/ferreus/geochem_pts_example.geoh5 +++ b/surface_apps-assets/ferreus/geochem_pts_example.geoh5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:7336127becd1c8fa141a6564049251c6db7dc51a14687eff95e120f5b6593715 +oid sha256:63aeb5733f7676170965f003baae8e9ffcd143e710d24535ba28793b218be869 size 1441792