From ef6804dcfc651420ed78d20288355cc2fe40b2cc Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Wed, 7 Oct 2026 11:29:07 -0300 Subject: [PATCH 1/3] feat: choose the number of interpolation elements from the data Add suggest_nelements and summarise_data, an Automatic setting, and one model manager method for the interpolator arguments that reads the saved settings instead of the class defaults. The feature panel shows the number and keeps a number that the user sets. --- docs/development/usability-plan.md | 106 +++++++++++--- docs/usage/interface.md | 7 + loopstructural/gui/dlg_settings.py | 7 + loopstructural/gui/dlg_settings.ui | 10 ++ .../feature_details_panel/_base.py | 50 ++++++- loopstructural/main/data_manager.py | 2 + loopstructural/main/interpolation_size.py | 132 ++++++++++++++++++ loopstructural/main/model_manager.py | 81 ++++++++--- loopstructural/toolbelt/preferences.py | 1 + tests/unit/test_interpolation_size.py | 107 ++++++++++++++ 10 files changed, 464 insertions(+), 39 deletions(-) create mode 100644 loopstructural/main/interpolation_size.py create mode 100644 tests/unit/test_interpolation_size.py diff --git a/docs/development/usability-plan.md b/docs/development/usability-plan.md index e14da21..4cbbc62 100644 --- a/docs/development/usability-plan.md +++ b/docs/development/usability-plan.md @@ -527,12 +527,12 @@ Design: Tasks: -- [ ] `main/interpolation_size.py` with `summarise_data` and `suggest_nelements`. -- [ ] Setting, settings page and preference test. -- [ ] Model manager: one method for the interpolator arguments, used in all +- [x] `main/interpolation_size.py` with `summarise_data` and `suggest_nelements`. +- [x] Setting, settings page and preference test. +- [x] Model manager: one method for the interpolator arguments, used in all build paths. Fix the use of the class defaults. -- [ ] Feature panel: show the number, and an "Automatic" check box. -- [ ] Docs: say how the number is chosen, and what the user can change. +- [x] Feature panel: show the number, and an "Automatic" check box. +- [x] Docs: say how the number is chosen, and what the user can change. Acceptance: a model with 20 contact points and a model with 5 000 points get different numbers of elements, both inside the limits. A saved fixed number is @@ -636,13 +636,34 @@ user can select any combination that is in the table below. frame). This adds the axial surface constraint (`fold_orientation`). 2. **Fold axis.** The user selects the source of the fold axis: - constant: plunge and azimuth (as now); - - average intersection lineation of the folded feature and the axial - surface (`av_fold_axis`, as now); - - lineation data: a layer of fold axis or intersection lineation - measurements. The plugin fits the fold axis rotation angle to - coordinate 1 (the axis S-plot). + - lineations. The user selects one of two lineation sources: + 1. **Lineation layer**: a point layer (for example a shapefile) with + measured fold axes or intersection lineations. The user selects the + plunge field and the trend (plunge direction) field. + 2. **Calculated intersection lineations**: the plugin calculates a + lineation at each orientation point of the folded feature. The + lineation is the intersection of the folded foliation and the axial + foliation (coordinate 0 of the fold frame) at that point + (`FoldFrame.calculate_intersection_lineation`). This source needs an + axial surface. + + Then the user selects how the plugin uses the lineations: + - **average**: the fold axis is the mean of the lineations, and it is + constant in the model (`av_fold_axis` for the calculated lineations, + as now); + - **fit**: the plugin fits the fold axis rotation angle of the + lineations to coordinate 1 (the axis S-plot), so that the fold axis + can change in the model. This needs an axial surface. This adds the fold axis constraint (`fold_axis_w`). + + | Lineation source | Average | Fit (axis S-plot) | Needs an axial surface | + |---|---|---|---| + | Lineation layer | yes | yes | only for "fit" | + | Calculated intersection lineations | yes | yes | yes | + + The UI shows the lineations of the selected source in the 3D view and in + the axis S-plot, so that the user can compare the two sources. 3. **S-plot.** The user controls the rotation angle profiles: the profile type (Fourier series, trigonometric), the wavelength, and fixed values for the profile parameters. Without this control, the plugin fits the @@ -652,11 +673,11 @@ user can select any combination that is in the table below. | Axial surface | Fold axis | S-plot | Result | |---|---|---|---| | - | - | - | A standard foliation. No fold constraint. | -| x | - | - | Fold frame and DFI. Axial surface constraint on. The fold axis is the average intersection lineation, but the fold axis constraint is off (`fold_axis_w = None`). Automatic profiles. | -| - | x | - | No fold frame. The plugin adds the fold axis as tangent constraints on a regular grid in the bounding box (gradient . axis = 0). The interpolator of the feature does not change. Only a constant fold axis is possible. | +| x | - | - | Fold frame and DFI. Axial surface constraint on. The fold axis is the average of the calculated intersection lineations, but the fold axis constraint is off (`fold_axis_w = None`). Automatic profiles. | +| - | x | - | No fold frame. The plugin adds the fold axis as tangent constraints on a regular grid in the bounding box (gradient . axis = 0). The interpolator of the feature does not change. Only a constant fold axis is possible: plunge and azimuth, or the average of a lineation layer. | | x | x | - | Fold frame and DFI. Both constraints on. Automatic profiles. | | x | - | x | As "axial surface only", but with the limb profile of the user. | -| x | x | x | All constraints on. The limb profile of the user. With "lineation data", also the axis profile of the user. | +| x | x | x | All constraints on. The limb profile of the user. With lineations and "fit", also the axis profile of the user. | | - | - | x, or - x x | Not possible. An S-plot needs a fold frame coordinate. The S-plot check box is disabled until the user selects an axial surface. The tooltip says why. | Rules: @@ -719,8 +740,11 @@ with orientations: `s2`, `s1` and `s0`. stratigraphic group: - `fold_event`: the name of the fold event, or `None`; - `axial_surface`: on or off, and `fold_orientation` weight; - - `fold_axis`: off, `constant` (plunge, azimuth), `average`, or `data` - (layer dicts), and `fold_axis_w` weight; + - `fold_axis`: `source` (off, `constant` or `lineations`), the plunge and + azimuth for `constant`, the `fold_axis_w` weight, and for `lineations`: + - `lineation_source`: `layer` (layer dict with the plunge and trend + fields) or `intersection` (calculated); + - `use`: `average` or `fit`; - `splot`: off or on; for the limb and for the axis profile: the type, the wavelength (or "automatic"), and the fixed parameters; - `fold_normalisation`, `fold_norm`, `fold_regularisation`. @@ -741,7 +765,8 @@ with orientations: `s2`, `s1` and `s0`. The S-plot needs the fold frame before the folded feature is built. Thus: - "Calculate rotation angles" builds the fold frame (and the fold events that - fold it) if it is not current. Then it calculates the rotation angles of the + fold it) if it is not current. Then it calculates the lineations (for the + "calculated intersection lineations" source) and the rotation angles of the data of the feature. It does not build the folded feature. Run it in a background task with progress. - The S-plot panel shows the data points, the fitted curve, the S-variogram @@ -791,8 +816,15 @@ Do the parts in this order. Each part is a separate pull request. the three controls with their check boxes and "Advanced" weights. - [ ] Apply the rules of the combination table. Disable the S-plot control when there is no axial surface. -- [ ] Fold axis "lineation data": a layer picker with the trend and plunge - fields. +- [ ] Fold axis "lineations": a choice of the two lineation sources. For + "lineation layer", a point layer picker with the plunge and trend + fields. For "calculated intersection lineations", no input (disabled + without an axial surface). +- [ ] Fold axis "average" or "fit" for the lineations. Disable "fit" without + an axial surface. +- [ ] Convert the plunge and trend of the layer to vectors, and give them + with their points (N x 6) to the fold axis calculation + (`main/fold_spec.py`). - [ ] Fold axis without axial surface: make the tangent constraints on a grid (`main/fold_spec.py`). The grid step comes from the bounding box and the number of elements. @@ -827,6 +859,11 @@ Tests: F2 -> F1 -> S0; a cycle is an error; each row of the combination table gives the correct arguments (a control that is off gives `None`); the S-plot is refused without an axial surface; old specs with `folded_feature_name` load. +- Unit tests for the two lineation sources: the plunge and trend of a layer + give the correct vectors; the calculated intersection lineation of a known + folded foliation and a known axial foliation is their cross product; the + "average" of both sources gives the same axis for the same data; "fit" and + the calculated source are refused without an axial surface. - A QGIS test that builds the refolded fold from the test layers and compares the S0 scalar field on a coarse grid with the result of the LoopStructural calls of the example (same data, same arguments). @@ -837,6 +874,22 @@ depend on 7.1. 7.4 depends on 7.3 (the S-plot control). 7.5 comes last. Phase 7 depends on 6.1 (the workflow mode), 6.3 (`SectionStack`) and 6.4 (the number of elements of each feature, which is also used for the fold frames). +### Phase 8: Advanced 3D viewer + +A second 3D viewer for hard 3D problems. It is a standalone application +(Rust, Bevy) with a live link to QGIS. The PyVista viewer stays as the default +viewer. The new viewer shows geological data objects that are similar to +those of Geoscience ANALYST (`geoh5` types): points, curves, surfaces, +sections, block models, drillholes and orientations. + +This phase is large, and most of the work is in a separate repository. The +tasks, the protocol, the risks and the open questions are in the +[viewer plan](viewer-plan.md). + +Acceptance: from step 5, a user opens the advanced viewer. The viewer shows +the model and its input data, and updates after each build. A pick in the +viewer selects the feature in the dock and shows the point on the map. + ## Risks - **Large UI change.** Users of the current version must learn the new layout. @@ -897,3 +950,20 @@ elements of each feature, which is also used for the fold frames). Does a fault need a different rule from a foliation? 8. (6.1) Must a manual unconformity or fold in the "constraints" mode use the column? Recommendation: no. It uses only the features of that mode. +9. (7) How does the plugin give the lineations of a layer to LoopStructural? + `FoldFrame.calculate_fold_axis_rotation` has a `fold_axis` argument (an + N x 6 array of points and lineations), but + `FoldedFeatureBuilder.set_fold_axis` does not use it, and there is no build + argument for it. "Average" of a lineation layer is easy: the plugin + calculates the mean and gives it as `fold_axis`. For "fit", either the + plugin calls `calculate_fold_axis_rotation` itself and sets + `fold.fold_axis_rotation`, or LoopStructural gets a build argument for the + lineations. Recommendation: add the build argument to LoopStructural, and + use the plugin call until that version is released. +10. (7) Is "fold axis without axial surface" (tangent constraints on a grid) + good enough, or must the plugin make a simple fold frame from the fold + axis? Recommendation: tangent constraints first. Compare the two on the + single fold example (`load_noddy_single_fold`). +11. (7) Must a fault that cuts a fold frame also cut the folded features? + LoopStructural calculates the rotation angles in the restored space. + Recommendation: yes, use the same faults. Test it before 7.5. diff --git a/docs/usage/interface.md b/docs/usage/interface.md index 54611f7..a9709a3 100644 --- a/docs/usage/interface.md +++ b/docs/usage/interface.md @@ -115,4 +115,11 @@ Units that don't match the stratigraphic column will have null values, helping y Once the layers have been selected, stratigraphic column defined and the fault topology relationships set, the LoopStructural model can be initialised. Initialise model will create a LoopStructural model with all of the geological features in the model. For each feature in the model the number of interpolation elements (degrees of freedom), the weighting of the regularisation, contact points and orientation weight can be changed. + +#### Number of interpolation elements + +By default the number of elements is "Automatic". The plugin chooses it for each feature from the data of the feature: more constraints, more surfaces and a wider spread of orientations give more elements. The number is rounded to 1 000 and kept between 5 000 and 250 000. The panel of the feature shows the number that was used. + +- To keep a number for one feature, clear the "Automatic" check box in the panel of the feature and enter the number. The feature keeps it, and the project file saves it, until you select "Automatic" again. +- To use one fixed number for all features, clear "Automatic" in the plugin settings and enter the number. Features that you set by hand keep their own number. ![Model Parameters](../static/model-setup.png) diff --git a/loopstructural/gui/dlg_settings.py b/loopstructural/gui/dlg_settings.py index 439768a..d719475 100644 --- a/loopstructural/gui/dlg_settings.py +++ b/loopstructural/gui/dlg_settings.py @@ -81,6 +81,10 @@ def __init__(self, parent): if hasattr(self, "btn_open_debug_directory"): self.btn_open_debug_directory.pressed.connect(self._open_debug_directory) + self.n_elements_auto_check_box.toggled.connect( + lambda checked: self.n_elements_spin_box.setEnabled(not checked) + ) + # load previously saved settings self.load_settings() @@ -92,6 +96,7 @@ def apply(self): settings.debug_mode = self.opt_debug.isChecked() settings.separate_dock_widgets = self.opt_separate_dock_widgets.isChecked() settings.interpolator_nelements = self.n_elements_spin_box.value() + settings.interpolator_nelements_auto = self.n_elements_auto_check_box.isChecked() settings.interpolator_npw = self.npw_spin_box.value() settings.interpolator_cpw = self.cpw_spin_box.value() settings.interpolator_regularisation = self.regularisation_spin_box.value() @@ -121,6 +126,8 @@ def load_settings(self): self.lbl_version_saved_value.setText(settings.version) # self.interpolator_type_combo.setCurrentText(settings.interpolator_type) self.n_elements_spin_box.setValue(settings.interpolator_nelements) + self.n_elements_auto_check_box.setChecked(settings.interpolator_nelements_auto) + self.n_elements_spin_box.setEnabled(not settings.interpolator_nelements_auto) self.regularisation_spin_box.setValue(settings.interpolator_regularisation) self.cpw_spin_box.setValue(settings.interpolator_cpw) self.npw_spin_box.setValue(settings.interpolator_npw) diff --git a/loopstructural/gui/dlg_settings.ui b/loopstructural/gui/dlg_settings.ui index 477f844..6a7c50e 100644 --- a/loopstructural/gui/dlg_settings.ui +++ b/loopstructural/gui/dlg_settings.ui @@ -280,6 +280,16 @@ + + + + Choose the number of elements of each feature from its data. If this is off, every feature uses the fixed number. + + + Automatic + + + diff --git a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py index 675eebc..300f59b 100644 --- a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py +++ b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py @@ -4,6 +4,7 @@ from qgis.gui import QgsCollapsibleGroupBox, QgsMapLayerComboBox from qgis.PyQt.QtCore import Qt from qgis.PyQt.QtWidgets import ( + QCheckBox, QComboBox, QDoubleSpinBox, QFormLayout, @@ -134,8 +135,16 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage self.n_elements_spinbox.setRange(100, 1000000) self.n_elements_spinbox.setValue(self.getNelements(feature)) self.n_elements_spinbox.setPrefix("Number of Elements: ") + self.n_elements_auto_check = QCheckBox("Automatic") + self.n_elements_auto_check.setToolTip( + "Choose the number of elements from the data of the feature. " + "Clear the check box to keep the number that you enter." + ) + self.n_elements_auto_check.setChecked(not self._has_nelements_override()) + self.n_elements_spinbox.setEnabled(not self.n_elements_auto_check.isChecked()) self.n_elements_spinbox.valueChanged.connect(self.updateNelements) + self.n_elements_auto_check.toggled.connect(self._on_nelements_auto_toggled) table_group_box = QgsCollapsibleGroupBox('Data Layers') self.layer_table = LayerSelectionTable( @@ -160,6 +169,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage form_layout = QFormLayout() form_layout.addRow(self.interpolator_type_label, self.interpolator_type_combo) form_layout.addRow("Number of Elements:", self.n_elements_spinbox) + form_layout.addRow("", self.n_elements_auto_check) form_layout.addRow('Regularisation', self.regularisation_spin_box) form_layout.addRow('Contact points weight', self.cpw_spin_box) form_layout.addRow('Orientation point weight', self.npw_spin_box) @@ -172,6 +182,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage summary=lambda: ( f"{self.interpolator_type_combo.currentText()}, " f"{int(self.n_elements_spinbox.value())} elements" + f"{' (automatic)' if self.n_elements_auto_check.isChecked() else ''}" ), ) self.layout.add_section(self._build_preview_widget(), 'preview', 'Preview', collapsed=True) @@ -565,8 +576,38 @@ def _on_bounding_box_updated(self, bounding_box): except Exception: pass - def updateNelements(self, value): - """Update the number of elements in the feature's interpolator.""" + def _has_nelements_override(self): + """Return True if the user set the number of elements of this feature.""" + manager = self.model_manager + name = getattr(self.feature, 'name', None) + return manager is not None and name in getattr(manager, 'nelements_overrides', {}) + + def _on_nelements_auto_toggled(self, automatic): + """Use the automatic number again, or keep the number that is shown.""" + manager = self.model_manager + name = getattr(self.feature, 'name', None) + self.n_elements_spinbox.setEnabled(not automatic) + if manager is None or name is None: + return + if automatic: + manager.set_nelements_override(name, None) + used = manager.nelements_used.get(name) + if used is not None: + self.n_elements_spinbox.blockSignals(True) + self.n_elements_spinbox.setValue(used[0]) + self.n_elements_spinbox.blockSignals(False) + self.updateNelements(used[0], keep=False) + else: + manager.set_nelements_override(name, int(self.n_elements_spinbox.value())) + + def updateNelements(self, value, keep=True): + """Update the number of elements in the feature's interpolator. + + With `keep`, the feature keeps this number in the next builds, until + the user selects "Automatic". + """ + if keep and self.model_manager is not None and self.feature is not None: + self.model_manager.set_nelements_override(self.feature.name, int(value)) if self.feature: if issubclass(type(self.feature), StructuralFrame): for i in range(3): @@ -587,6 +628,11 @@ def updateNelements(self, value): def getNelements(self, feature): """Get the number of elements from the feature's interpolator.""" if feature: + used = getattr(self.model_manager, 'nelements_used', {}).get( + getattr(feature, 'name', None) + ) + if used is not None: + return used[0] if issubclass(type(feature), StructuralFrame): return feature[0].interpolator.n_elements elif feature.interpolator is not None: diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index 7e7fd7d..281bd18 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -1848,6 +1848,7 @@ def save_state(self, filepath): state['manual_foliations'] = self._model_manager.manual_foliations_to_dict() state['detached_features'] = self._model_manager.detached_to_dict() state['parametric_faults'] = self._model_manager.parametric_faults_to_dict() + state['nelements_overrides'] = self._model_manager.nelements_overrides_to_dict() with open(path, 'w') as f: json.dump(state, f, indent=2) @@ -1882,6 +1883,7 @@ def load_state(self, filepath): self._model_manager.manual_foliations_from_dict(state.get('manual_foliations', {})) self._model_manager.detached_from_dict(state.get('detached_features', {})) self._model_manager.parametric_faults_from_dict(state.get('parametric_faults', {})) + self._model_manager.nelements_overrides_from_dict(state.get('nelements_overrides', {})) # the data was just read from the layers self._changed_layer_ids.clear() self.refresh_layer_watchers() diff --git a/loopstructural/main/interpolation_size.py b/loopstructural/main/interpolation_size.py new file mode 100644 index 0000000..c01819b --- /dev/null +++ b/loopstructural/main/interpolation_size.py @@ -0,0 +1,132 @@ +"""Choose the number of interpolation elements from the data of a feature. + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +from dataclasses import dataclass + +import numpy as np + +MIN_NELEMENTS = 5_000 +MAX_NELEMENTS = 250_000 +ROUND_TO = 1_000 +ELEMENTS_PER_EQUATION = 25 +MAX_SURFACES = 10 +SURFACE_STEP = 0.1 + + +@dataclass +class DataSummary: + """The measures of the data of one feature. + + `spread` is 0 for parallel planes and 1 for orientations in all + directions. + """ + + n_value: int = 0 + n_orientation: int = 0 + n_surfaces: int = 0 + spread: float = 0.0 + + +def _normals_from_strike_dip(strike, dip) -> np.ndarray: + """Return the unit normals of planes (right-hand rule strike, dip in degrees).""" + strike = np.radians(np.asarray(strike, dtype=float)) + dip = np.radians(np.asarray(dip, dtype=float)) + dip_direction = strike + np.pi / 2 + return np.column_stack( + [ + -np.sin(dip) * np.sin(dip_direction), + -np.sin(dip) * np.cos(dip_direction), + np.cos(dip), + ] + ) + + +def orientation_spread(normals) -> float: + """Return the spread of orientations: 0 (parallel) to 1 (all directions). + + The sign of a normal does not change the result: the measure uses the + largest eigenvalue of the mean orientation tensor, which is 1 for parallel + planes and 1/3 for uniform directions. + """ + normals = np.asarray(normals, dtype=float).reshape(-1, 3) + normals = normals[np.all(np.isfinite(normals), axis=1)] + lengths = np.linalg.norm(normals, axis=1) + normals = normals[lengths > 0] / lengths[lengths > 0, None] + if len(normals) < 2: + return 0.0 + tensor = normals.T @ normals / len(normals) + largest = np.linalg.eigvalsh(tensor)[-1] + return float(np.clip((1.0 - largest) / (2.0 / 3.0), 0.0, 1.0)) + + +def _has(df, columns) -> np.ndarray: + """Return a mask of the rows where all `columns` exist and are not NaN.""" + mask = np.ones(len(df), dtype=bool) + for column in columns: + if column not in df.columns: + return np.zeros(len(df), dtype=bool) + mask &= df[column].notna().to_numpy() + return mask + + +def summarise_data(df) -> DataSummary: + """Return the `DataSummary` of a feature data frame (or None for no data). + + Value, interface and inequality rows are value constraints. Rows with a + normal (nx, ny, nz or gx, gy, gz), a strike and a dip, or a tangent (tx, + ty, tz) are orientation constraints. + """ + if df is None or len(df) == 0: + return DataSummary() + value = _has(df, ['val']) | _has(df, ['interface']) | _has(df, ['l', 'u']) + normal = _has(df, ['nx', 'ny', 'nz']) + gradient = _has(df, ['gx', 'gy', 'gz']) & ~normal + strike_dip = _has(df, ['strike', 'dip']) & ~normal & ~gradient + tangent = _has(df, ['tx', 'ty', 'tz']) + orientation = normal | gradient | strike_dip | tangent + + vectors = [] + if normal.any(): + vectors.append(df.loc[normal, ['nx', 'ny', 'nz']].to_numpy(dtype=float)) + if gradient.any(): + vectors.append(df.loc[gradient, ['gx', 'gy', 'gz']].to_numpy(dtype=float)) + if strike_dip.any(): + vectors.append( + _normals_from_strike_dip( + df.loc[strike_dip, 'strike'].to_numpy(), df.loc[strike_dip, 'dip'].to_numpy() + ) + ) + spread = orientation_spread(np.vstack(vectors)) if vectors else 0.0 + + n_surfaces = 0 + if _has(df, ['val']).any(): + n_surfaces = int(df.loc[_has(df, ['val']), 'val'].nunique()) + if 'interface' in df.columns and df['interface'].notna().any(): + n_surfaces = max(n_surfaces, int(df['interface'].nunique())) + return DataSummary( + n_value=int(value.sum()), + n_orientation=int(orientation.sum()), + n_surfaces=n_surfaces, + spread=spread, + ) + + +def suggest_nelements(summary: DataSummary) -> int: + """Return the number of elements for a feature with the data `summary`. + + elements = equations x 25 x surface factor x spread factor, where an + orientation is two equations, the surface factor is 1 + 0.1 for each + surface after the first (at most 10 surfaces) and the spread factor is + 1 + spread. The result is rounded to 1 000 and limited to 5 000 .. 250 000. + """ + equations = summary.n_value + 2 * summary.n_orientation + if equations <= 0: + return MIN_NELEMENTS + surfaces = min(max(summary.n_surfaces, 1), MAX_SURFACES) + surface_factor = 1 + SURFACE_STEP * (surfaces - 1) + spread_factor = 1 + min(max(summary.spread, 0.0), 1.0) + elements = equations * ELEMENTS_PER_EQUATION * surface_factor * spread_factor + elements = int(round(elements / ROUND_TO)) * ROUND_TO + return int(min(max(elements, MIN_NELEMENTS), MAX_NELEMENTS)) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index a241295..cc59778 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -27,11 +27,12 @@ from LoopStructural.utils.observer import Observable from LoopStructural import GeologicalModel -from loopstructural.toolbelt.preferences import PlgSettingsStructure +from loopstructural.toolbelt.preferences import PlgOptionsManager from ..main import constraints, parametric_fault from ..main.data_types import FaultEntry, StratigraphyEntry from ..main.helpers import qgisAttributeIsNone +from ..main.interpolation_size import suggest_nelements, summarise_data from ..main.workflow_mode import WORKFLOW_MODE_CONSTRAINTS, WORKFLOW_MODE_MAP, WORKFLOW_MODES @@ -236,6 +237,13 @@ def __init__(self, debug_manager=None): self._cancel_requested = False # How the user builds the model; see `set_workflow_mode`. self.workflow_mode = WORKFLOW_MODE_MAP + # feature name -> number of elements that the user set for the + # feature. A feature that is not here uses the automatic number or the + # fixed number of the settings; see `interpolator_arguments`. + self.nelements_overrides: Dict[str, int] = {} + # feature name -> (number of elements, True if automatic) of the last + # build, for the feature panel. + self.nelements_used: Dict[str, tuple] = {} @property def uses_map_data(self) -> bool: @@ -309,6 +317,8 @@ def reset(self): self.detached = {} self.extra_constraints = {} self.parametric_faults = {} + self.nelements_overrides = {} + self.nelements_used = {} self.dem_function = lambda x, y: 0 self._topology_dirty = False self._data_dirty = False @@ -977,14 +987,53 @@ def _build_domain_fault_boundary(self, fault_name, groupname, flipped=False): self.model.data = data_for_fault self.model.create_and_add_domain_fault( fault_name, - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(fault_name, data_for_fault), ) self._strip_domain_fault_region_from_boundary_above(fault_name) return True + def interpolator_arguments(self, name, data) -> dict: + """Return `nelements`, `npw`, `cpw` and `regularisation` for the + feature `name` that is built from the data frame `data`. + + The number of elements is, in this order: the number that the user set + for the feature, the number from the data (if the saved setting + "automatic" is on), or the saved fixed number. All values come from the + saved settings, not from the class defaults. + """ + settings = PlgOptionsManager.get_plg_settings() + automatic = False + if self.nelements_overrides.get(name): + nelements = int(self.nelements_overrides[name]) + elif settings.interpolator_nelements_auto: + nelements = suggest_nelements(summarise_data(data)) + automatic = True + else: + nelements = int(settings.interpolator_nelements) + self.nelements_used[name] = (nelements, automatic) + return { + 'nelements': nelements, + 'npw': settings.interpolator_npw, + 'cpw': settings.interpolator_cpw, + 'regularisation': settings.interpolator_regularisation, + } + + def nelements_overrides_to_dict(self) -> dict: + """Return `nelements_overrides` in a form that `json.dump` can write.""" + return {name: int(n) for name, n in self.nelements_overrides.items()} + + def nelements_overrides_from_dict(self, overrides: dict): + """Replace `nelements_overrides` with the numbers from `nelements_overrides_to_dict`.""" + self.nelements_overrides = {name: int(n) for name, n in (overrides or {}).items()} + + def set_nelements_override(self, name: str, nelements: Optional[int]): + """Keep `nelements` for the feature `name` in the next builds, or use + the automatic (or fixed) number again if `nelements` is None.""" + if nelements is None: + self.nelements_overrides.pop(name, None) + else: + self.nelements_overrides[name] = int(nelements) + def _strip_domain_fault_region_from_boundary_above(self, fault_name): """Work around a LoopStructural core bug that crops the boundary above a domain fault. @@ -1064,10 +1113,7 @@ def update_foliation_features(self): data=data, force_constrained=True, **self._extra_kwargs(groupname), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(groupname, data), ) fault_name, flipped = self._base_fault_boundary(group) if fault_name is None or not self._build_domain_fault_boundary( @@ -1413,10 +1459,7 @@ def update_fault_features(self): fault_pitch=pitch, data=data, faults=self._cutting_faults_for(fault_name), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + **self.interpolator_arguments(fault_name, data), ) self._build_parametric_faults() self.apply_fault_abutting_relationships() @@ -1458,13 +1501,11 @@ def _build_parametric_faults(self): for name, spec in self.parametric_faults.items(): self._report_progress(f"Building fault '{name}'") try: + fault_data = parametric_fault.frame_data(spec) self.model.create_and_add_fault( name, - data=parametric_fault.frame_data(spec), - nelements=PlgSettingsStructure.interpolator_nelements, - npw=PlgSettingsStructure.interpolator_npw, - cpw=PlgSettingsStructure.interpolator_cpw, - regularisation=PlgSettingsStructure.interpolator_regularisation, + data=fault_data, + **self.interpolator_arguments(name, fault_data), **parametric_fault.fault_arguments(spec), ) except Exception as e: @@ -2105,7 +2146,9 @@ def _create_foliation( ): """Create the foliation feature in the model; see `add_foliation`.""" data, kwargs = self._foliation_data(name, data, sampler, use_z_coordinate) - foliation = self.model.create_and_add_foliation(name, data=data, **kwargs) + foliation = self.model.create_and_add_foliation( + name, data=data, **kwargs, **self.interpolator_arguments(name, data) + ) if not restrict_to_stratigraphic_domain: foliation.regions = [ r for r in foliation.regions if not isinstance(r, UnconformityFeature) diff --git a/loopstructural/toolbelt/preferences.py b/loopstructural/toolbelt/preferences.py index 083ce8b..158f453 100644 --- a/loopstructural/toolbelt/preferences.py +++ b/loopstructural/toolbelt/preferences.py @@ -32,6 +32,7 @@ class PlgSettingsStructure: version: str = __version__ interpolator_type: str = 'FDI' interpolator_nelements: int = 50000 + interpolator_nelements_auto: bool = True interpolator_regularisation: float = 1.0 interpolator_cpw: float = 1.0 interpolator_npw: float = 1.0 diff --git a/tests/unit/test_interpolation_size.py b/tests/unit/test_interpolation_size.py new file mode 100644 index 0000000..4eb9f8b --- /dev/null +++ b/tests/unit/test_interpolation_size.py @@ -0,0 +1,107 @@ +"""Pytest tests for the choice of the number of interpolation elements. + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +import numpy as np +import pandas as pd + +from loopstructural.main.interpolation_size import ( + MAX_NELEMENTS, + MIN_NELEMENTS, + DataSummary, + orientation_spread, + suggest_nelements, + summarise_data, +) + + +class TestSuggestNelements: + def test_no_data_gives_the_minimum(self): + assert suggest_nelements(DataSummary()) == MIN_NELEMENTS + + def test_more_data_gives_more_or_equal_elements(self): + previous = 0 + for n in (1, 10, 100, 1_000, 10_000, 100_000): + value = suggest_nelements(DataSummary(n_value=n, n_orientation=n, n_surfaces=3)) + assert value >= previous + previous = value + + def test_limits(self): + assert suggest_nelements(DataSummary(n_value=1)) == MIN_NELEMENTS + assert suggest_nelements(DataSummary(n_value=10**7)) == MAX_NELEMENTS + + def test_result_is_a_multiple_of_1000(self): + assert suggest_nelements(DataSummary(n_value=777, n_orientation=33)) % 1000 == 0 + + def test_spread_gives_more_elements(self): + parallel = suggest_nelements(DataSummary(n_value=200, n_orientation=200, spread=0)) + spread = suggest_nelements(DataSummary(n_value=200, n_orientation=200, spread=1)) + assert spread > parallel + + def test_surfaces_give_more_elements_and_stop_at_ten(self): + one = suggest_nelements(DataSummary(n_value=200, n_surfaces=1)) + five = suggest_nelements(DataSummary(n_value=200, n_surfaces=5)) + ten = suggest_nelements(DataSummary(n_value=200, n_surfaces=10)) + twenty = suggest_nelements(DataSummary(n_value=200, n_surfaces=20)) + assert one < five < ten + assert ten == twenty + + def test_small_and_large_data_differ(self): + small = summarise_data(_values(20)) + large = summarise_data(_values(5_000)) + assert suggest_nelements(small) < suggest_nelements(large) + + +def _values(n): + return pd.DataFrame({'X': range(n), 'val': [0.0] * n, 'feature_name': 'a'}) + + +class TestOrientationSpread: + def test_parallel_planes_have_no_spread(self): + assert orientation_spread(np.tile([0, 0, 1], (10, 1))) == 0 + + def test_all_directions_have_spread_near_one(self): + rng = np.random.default_rng(0) + normals = rng.normal(size=(5_000, 3)) + assert orientation_spread(normals) > 0.95 + + def test_sign_does_not_matter(self): + rng = np.random.default_rng(1) + normals = rng.normal(size=(50, 3)) + flipped = normals.copy() + flipped[::2] *= -1 + assert np.isclose(orientation_spread(normals), orientation_spread(flipped)) + + def test_one_orientation_has_no_spread(self): + assert orientation_spread([[1, 0, 0]]) == 0 + + +class TestSummariseData: + def test_none_and_empty(self): + assert summarise_data(None) == DataSummary() + assert summarise_data(pd.DataFrame()) == DataSummary() + + def test_counts_and_surfaces(self): + df = pd.DataFrame( + { + 'val': [0.0, 0.0, 10.0, np.nan], + 'nx': [np.nan, np.nan, np.nan, 0.0], + 'ny': [np.nan, np.nan, np.nan, 0.0], + 'nz': [np.nan, np.nan, np.nan, 1.0], + } + ) + summary = summarise_data(df) + assert summary.n_value == 3 + assert summary.n_orientation == 1 + assert summary.n_surfaces == 2 + + def test_strike_dip_spread(self): + parallel = pd.DataFrame({'strike': [10.0] * 5, 'dip': [30.0] * 5}) + mixed = pd.DataFrame({'strike': [0.0, 90.0, 180.0, 270.0], 'dip': [60.0] * 4}) + assert summarise_data(parallel).spread == 0 + assert summarise_data(mixed).spread > 0.3 + + def test_opposite_dip_directions_spread(self): + a = summarise_data(pd.DataFrame({'strike': [0.0, 180.0], 'dip': [30.0, 30.0]})) + assert a.spread > 0 From 4c318fe261611a810825bff8d22a26f2ed378a86 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 30 Oct 2026 19:28:58 +1030 Subject: [PATCH 2/3] fix: element count from the extent, and the Automatic check box - Calculate the element count from the model extent as well as from the number of data rows: the volume of the model divided by the cube of the median spacing of the data points. Use the larger number. - When the user selects Automatic again, show the number of the settings, not the last number of the user. - The check box follows the global Automatic setting. Without it, the box uses the fixed number of the settings. - Correct the sign of the normal in _normals_from_strike_dip. - Add tests for interpolator_arguments and for the save and load of the overrides. --- .../feature_details_panel/_base.py | 44 +++++++--- loopstructural/main/interpolation_size.py | 54 ++++++++++-- loopstructural/main/model_manager.py | 38 +++++++-- tests/qgis/test_interpolator_arguments.py | 84 +++++++++++++++++++ tests/unit/test_interpolation_size.py | 39 +++++++++ 5 files changed, 231 insertions(+), 28 deletions(-) create mode 100644 tests/qgis/test_interpolator_arguments.py diff --git a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py index 83342ec..0310417 100644 --- a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py +++ b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py @@ -22,6 +22,7 @@ from ....background_task import finish_background_task, start_background_task from ....messages import push_info, push_warning +from .....toolbelt.preferences import PlgOptionsManager from ...steps.section_stack import SectionStack from ..bounding_box_widget import BoundingBoxWidget from ..layer_selection_table import LayerSelectionTable @@ -134,11 +135,20 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage self.n_elements_spinbox.setRange(100, 1000000) self.n_elements_spinbox.setValue(self.getNelements(feature)) self.n_elements_spinbox.setPrefix("Number of Elements: ") - self.n_elements_auto_check = QCheckBox("Automatic") - self.n_elements_auto_check.setToolTip( - "Choose the number of elements from the data of the feature. " - "Clear the check box to keep the number that you enter." - ) + # The check box uses the number of the settings: automatic from the + # data, or the fixed number, as the global "Automatic" setting gives + if PlgOptionsManager.get_plg_settings().interpolator_nelements_auto: + self.n_elements_auto_check = QCheckBox("Automatic") + self.n_elements_auto_check.setToolTip( + "Choose the number of elements from the data of the feature. " + "Clear the check box to keep the number that you enter." + ) + else: + self.n_elements_auto_check = QCheckBox("Use the number of the settings") + self.n_elements_auto_check.setToolTip( + "Use the fixed number of elements of the plugin settings. " + "Clear the check box to keep the number that you enter." + ) self.n_elements_auto_check.setChecked(not self._has_nelements_override()) self.n_elements_spinbox.setEnabled(not self.n_elements_auto_check.isChecked()) @@ -181,7 +191,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage summary=lambda: ( f"{self.interpolator_type_combo.currentText()}, " f"{int(self.n_elements_spinbox.value())} elements" - f"{' (automatic)' if self.n_elements_auto_check.isChecked() else ''}" + f"{' (automatic)' if self._nelements_are_automatic() else ''}" ), ) self.layout.add_section(self._build_preview_widget(), 'preview', 'Preview', collapsed=True) @@ -581,8 +591,15 @@ def _has_nelements_override(self): name = getattr(self.feature, 'name', None) return manager is not None and name in getattr(manager, 'nelements_overrides', {}) + def _nelements_are_automatic(self): + """Return True if the number of elements comes from the data.""" + return ( + self.n_elements_auto_check.isChecked() + and PlgOptionsManager.get_plg_settings().interpolator_nelements_auto + ) + def _on_nelements_auto_toggled(self, automatic): - """Use the automatic number again, or keep the number that is shown.""" + """Use the number of the settings again, or keep the number that is shown.""" manager = self.model_manager name = getattr(self.feature, 'name', None) self.n_elements_spinbox.setEnabled(not automatic) @@ -590,12 +607,13 @@ def _on_nelements_auto_toggled(self, automatic): return if automatic: manager.set_nelements_override(name, None) - used = manager.nelements_used.get(name) - if used is not None: - self.n_elements_spinbox.blockSignals(True) - self.n_elements_spinbox.setValue(used[0]) - self.n_elements_spinbox.blockSignals(False) - self.updateNelements(used[0], keep=False) + # not `nelements_used`: after a build with the number of the user, + # it has that number + nelements, _automatic = manager.default_nelements(name) + self.n_elements_spinbox.blockSignals(True) + self.n_elements_spinbox.setValue(nelements) + self.n_elements_spinbox.blockSignals(False) + self.updateNelements(nelements, keep=False) else: manager.set_nelements_override(name, int(self.n_elements_spinbox.value())) diff --git a/loopstructural/main/interpolation_size.py b/loopstructural/main/interpolation_size.py index c01819b..f875aa2 100644 --- a/loopstructural/main/interpolation_size.py +++ b/loopstructural/main/interpolation_size.py @@ -20,24 +20,30 @@ class DataSummary: """The measures of the data of one feature. `spread` is 0 for parallel planes and 1 for orientations in all - directions. + directions. `spacing` is the median distance from a data point to the + nearest other data point (0 if there are fewer than two points). """ n_value: int = 0 n_orientation: int = 0 n_surfaces: int = 0 spread: float = 0.0 + spacing: float = 0.0 def _normals_from_strike_dip(strike, dip) -> np.ndarray: - """Return the unit normals of planes (right-hand rule strike, dip in degrees).""" + """Return the upward unit normals of planes (right-hand rule strike, dip in degrees). + + The horizontal part of a normal points in the dip direction, as in + LoopStructural's `strikedip2vector`. + """ strike = np.radians(np.asarray(strike, dtype=float)) dip = np.radians(np.asarray(dip, dtype=float)) dip_direction = strike + np.pi / 2 return np.column_stack( [ - -np.sin(dip) * np.sin(dip_direction), - -np.sin(dip) * np.cos(dip_direction), + np.sin(dip) * np.sin(dip_direction), + np.sin(dip) * np.cos(dip_direction), np.cos(dip), ] ) @@ -61,6 +67,23 @@ def orientation_spread(normals) -> float: return float(np.clip((1.0 - largest) / (2.0 / 3.0), 0.0, 1.0)) +def point_spacing(points) -> float: + """Return the median distance from each point to the nearest other point. + + Points at the same location count as one point. Returns 0 for fewer than + two points. + """ + from scipy.spatial import cKDTree + + points = np.asarray(points, dtype=float).reshape(-1, 3) + points = points[np.all(np.isfinite(points), axis=1)] + points = np.unique(points, axis=0) + if len(points) < 2: + return 0.0 + distances, _ = cKDTree(points).query(points, k=2) + return float(np.median(distances[:, 1])) + + def _has(df, columns) -> np.ndarray: """Return a mask of the rows where all `columns` exist and are not NaN.""" mask = np.ones(len(df), dtype=bool) @@ -105,21 +128,32 @@ def summarise_data(df) -> DataSummary: n_surfaces = int(df.loc[_has(df, ['val']), 'val'].nunique()) if 'interface' in df.columns and df['interface'].notna().any(): n_surfaces = max(n_surfaces, int(df['interface'].nunique())) + spacing = 0.0 + if all(column in df.columns for column in ('X', 'Y', 'Z')): + spacing = point_spacing(df[['X', 'Y', 'Z']].to_numpy(dtype=float)) return DataSummary( n_value=int(value.sum()), n_orientation=int(orientation.sum()), n_surfaces=n_surfaces, spread=spread, + spacing=spacing, ) -def suggest_nelements(summary: DataSummary) -> int: +def suggest_nelements(summary: DataSummary, volume: float = 0.0) -> int: """Return the number of elements for a feature with the data `summary`. - elements = equations x 25 x surface factor x spread factor, where an - orientation is two equations, the surface factor is 1 + 0.1 for each - surface after the first (at most 10 surfaces) and the spread factor is - 1 + spread. The result is rounded to 1 000 and limited to 5 000 .. 250 000. + From the number of data rows: elements = equations x 25 x surface factor + x spread factor, where an orientation is two equations, the surface + factor is 1 + 0.1 for each surface after the first (at most 10 surfaces) + and the spread factor is 1 + spread. + + From the model extent: if `volume` (the volume of the model box) and the + data spacing are known, elements = volume / spacing^3, so that an element + is not larger than the distance between the data points. + + The result is the larger of the two, rounded to 1 000 and limited to + 5 000 .. 250 000. """ equations = summary.n_value + 2 * summary.n_orientation if equations <= 0: @@ -128,5 +162,7 @@ def suggest_nelements(summary: DataSummary) -> int: surface_factor = 1 + SURFACE_STEP * (surfaces - 1) spread_factor = 1 + min(max(summary.spread, 0.0), 1.0) elements = equations * ELEMENTS_PER_EQUATION * surface_factor * spread_factor + if volume > 0 and summary.spacing > 0: + elements = max(elements, min(volume / summary.spacing**3, MAX_NELEMENTS)) elements = int(round(elements / ROUND_TO)) * ROUND_TO return int(min(max(elements, MIN_NELEMENTS), MAX_NELEMENTS)) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 0a76366..d7bc1f7 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -32,7 +32,7 @@ from ..main import constraints, parametric_fault from ..main.data_types import FaultEntry, StratigraphyEntry from ..main.helpers import qgisAttributeIsNone -from ..main.interpolation_size import suggest_nelements, summarise_data +from ..main.interpolation_size import DataSummary, suggest_nelements, summarise_data from ..main.workflow_mode import WORKFLOW_MODE_CONSTRAINTS, WORKFLOW_MODE_MAP, WORKFLOW_MODES @@ -246,6 +246,8 @@ def __init__(self, debug_manager=None): # feature name -> (number of elements, True if automatic) of the last # build, for the feature panel. self.nelements_used: Dict[str, tuple] = {} + # feature name -> DataSummary of the data of the last build + self._data_summaries: Dict[str, DataSummary] = {} @property def uses_map_data(self) -> bool: @@ -322,6 +324,7 @@ def reset(self): self.parametric_faults = {} self.nelements_overrides = {} self.nelements_used = {} + self._data_summaries = {} self.dem_function = lambda x, y: 0 self._topology_dirty = False self._data_dirty = False @@ -1005,14 +1008,12 @@ def interpolator_arguments(self, name, data) -> dict: saved settings, not from the class defaults. """ settings = PlgOptionsManager.get_plg_settings() - automatic = False + self._data_summaries[name] = summarise_data(data) if self.nelements_overrides.get(name): nelements = int(self.nelements_overrides[name]) - elif settings.interpolator_nelements_auto: - nelements = suggest_nelements(summarise_data(data)) - automatic = True + automatic = False else: - nelements = int(settings.interpolator_nelements) + nelements, automatic = self.default_nelements(name) self.nelements_used[name] = (nelements, automatic) return { 'nelements': nelements, @@ -1021,6 +1022,31 @@ def interpolator_arguments(self, name, data) -> dict: 'regularisation': settings.interpolator_regularisation, } + def default_nelements(self, name) -> tuple: + """Return (number of elements, True if automatic) for the feature + `name` without the number that the user set for it. + + With the saved setting "automatic", the number comes from the data of + the last build of the feature. Otherwise it is the saved fixed number. + """ + settings = PlgOptionsManager.get_plg_settings() + if not settings.interpolator_nelements_auto: + return int(settings.interpolator_nelements), False + summary = self._data_summaries.get(name, DataSummary()) + return suggest_nelements(summary, self._model_volume()), True + + def _model_volume(self) -> float: + """Return the volume of the model box, or 0 if the box is not set.""" + bounding_box = getattr(self.model, 'bounding_box', None) + try: + lengths = np.asarray(bounding_box.maximum, dtype=float) - np.asarray( + bounding_box.origin, dtype=float + ) + except (AttributeError, TypeError, ValueError): + return 0.0 + volume = float(np.prod(lengths)) + return volume if np.isfinite(volume) and volume > 0 else 0.0 + def nelements_overrides_to_dict(self) -> dict: """Return `nelements_overrides` in a form that `json.dump` can write.""" return {name: int(n) for name, n in self.nelements_overrides.items()} diff --git a/tests/qgis/test_interpolator_arguments.py b/tests/qgis/test_interpolator_arguments.py new file mode 100644 index 0000000..45d5b56 --- /dev/null +++ b/tests/qgis/test_interpolator_arguments.py @@ -0,0 +1,84 @@ +"""Pytest tests for the number of elements that `interpolator_arguments` gives. + +The order is: the number that the user set for the feature, then the number +from the data (with the setting "automatic"), then the fixed number. +""" + +import json + +import pandas as pd +import pytest +from LoopStructural.datatypes import BoundingBox + +from loopstructural.main import model_manager as model_manager_module +from loopstructural.main.interpolation_size import suggest_nelements, summarise_data +from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.toolbelt.preferences import PlgSettingsStructure + +DATA = pd.DataFrame( + { + 'X': [10.0, 50.0, 90.0], + 'Y': [10.0, 50.0, 90.0], + 'Z': [0.0, 0.0, 0.0], + 'val': [0.0, 0.0, 1.0], + } +) + + +class _DebugManager: + def log(self, *args, **kwargs): + pass + + +def _use_settings(monkeypatch, **values): + settings = PlgSettingsStructure(**values) + monkeypatch.setattr( + model_manager_module.PlgOptionsManager, 'get_plg_settings', staticmethod(lambda: settings) + ) + + +@pytest.fixture +def manager(): + manager = GeologicalModelManager(debug_manager=_DebugManager()) + manager.update_bounding_box(BoundingBox(origin=[0, 0, -50], maximum=[100, 100, 50])) + return manager + + +def test_the_number_of_the_user_comes_first(manager, monkeypatch): + _use_settings(monkeypatch, interpolator_nelements_auto=True, interpolator_nelements=7000) + manager.set_nelements_override('s1', 12345) + assert manager.interpolator_arguments('s1', DATA)['nelements'] == 12345 + assert manager.nelements_used['s1'] == (12345, False) + + +def test_automatic_uses_the_data(manager, monkeypatch): + _use_settings(monkeypatch, interpolator_nelements_auto=True, interpolator_nelements=7000) + expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0 * 100.0) + assert manager.interpolator_arguments('s1', DATA)['nelements'] == expected + assert manager.nelements_used['s1'] == (expected, True) + + +def test_the_fixed_number_without_automatic(manager, monkeypatch): + _use_settings(monkeypatch, interpolator_nelements_auto=False, interpolator_nelements=7000) + assert manager.interpolator_arguments('s1', DATA)['nelements'] == 7000 + assert manager.nelements_used['s1'] == (7000, False) + + +def test_the_default_is_not_the_number_of_the_user(manager, monkeypatch): + _use_settings(monkeypatch, interpolator_nelements_auto=True, interpolator_nelements=7000) + manager.set_nelements_override('s1', 12345) + manager.interpolator_arguments('s1', DATA) + expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0 * 100.0) + assert manager.default_nelements('s1') == (expected, True) + + +def test_the_overrides_are_saved_and_loaded(manager): + manager.set_nelements_override('s1', 12345) + manager.set_nelements_override('f1', 6000) + saved = json.loads(json.dumps(manager.nelements_overrides_to_dict())) + + other = GeologicalModelManager(debug_manager=_DebugManager()) + other.nelements_overrides_from_dict(saved) + assert other.nelements_overrides == {'s1': 12345, 'f1': 6000} + other.set_nelements_override('s1', None) + assert other.nelements_overrides == {'f1': 6000} diff --git a/tests/unit/test_interpolation_size.py b/tests/unit/test_interpolation_size.py index 4eb9f8b..08ff0e2 100644 --- a/tests/unit/test_interpolation_size.py +++ b/tests/unit/test_interpolation_size.py @@ -10,7 +10,9 @@ MAX_NELEMENTS, MIN_NELEMENTS, DataSummary, + _normals_from_strike_dip, orientation_spread, + point_spacing, suggest_nelements, summarise_data, ) @@ -105,3 +107,40 @@ def test_strike_dip_spread(self): def test_opposite_dip_directions_spread(self): a = summarise_data(pd.DataFrame({'strike': [0.0, 180.0], 'dip': [30.0, 30.0]})) assert a.spread > 0 + + +class TestNormals: + def test_the_normal_points_in_the_dip_direction(self): + # strike 0 (north) dips to the east, so the normal points east and up + normal = _normals_from_strike_dip([0.0], [30.0])[0] + assert normal[0] > 0 + assert abs(normal[1]) < 1e-12 + assert normal[2] > 0 + assert np.isclose(np.linalg.norm(normal), 1.0) + + +class TestExtent: + def test_point_spacing(self): + points = [(x, 0.0, 0.0) for x in (0.0, 10.0, 20.0, 30.0)] + assert point_spacing(points) == 10.0 + + def test_repeated_points_count_once(self): + assert point_spacing([(0.0, 0.0, 0.0)] * 3 + [(5.0, 0.0, 0.0)]) == 5.0 + assert point_spacing([(0.0, 0.0, 0.0)] * 3) == 0.0 + + def test_summary_has_the_spacing(self): + df = pd.DataFrame({'X': [0.0, 50.0, 100.0], 'Y': 0.0, 'Z': 0.0, 'val': 0.0}) + assert summarise_data(df).spacing == 50.0 + + def test_a_large_model_with_dense_data_gives_more_elements(self): + summary = DataSummary(n_value=20, spacing=100.0) + from_rows = suggest_nelements(summary) + small = suggest_nelements(summary, volume=1000.0**3) + large = suggest_nelements(summary, volume=10_000.0 * 10_000.0 * 2_000.0) + assert small == from_rows + assert large > from_rows + assert large <= MAX_NELEMENTS + + def test_no_spacing_uses_only_the_rows(self): + summary = DataSummary(n_value=20) + assert suggest_nelements(summary, volume=1e12) == suggest_nelements(summary) From 7ec25a1aed0fdf5a3361069d6eb991ad9a88d733 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sat, 31 Oct 2026 11:32:29 +1030 Subject: [PATCH 3/3] fix: limit the extent factor of the automatic element count volume / spacing^3 gave the maximum (250 000) for most map data: contacts are sampled densely along lines, so the spacing is small. The solve was then slow. The extent factor is now the map area of the model divided by the map area of the data, limited to 1 .. 2. --- loopstructural/main/interpolation_size.py | 57 +++++++++++------------ loopstructural/main/model_manager.py | 16 +++---- tests/qgis/test_interpolator_arguments.py | 4 +- tests/unit/test_interpolation_size.py | 49 +++++++++++-------- 4 files changed, 65 insertions(+), 61 deletions(-) diff --git a/loopstructural/main/interpolation_size.py b/loopstructural/main/interpolation_size.py index f875aa2..6e1e550 100644 --- a/loopstructural/main/interpolation_size.py +++ b/loopstructural/main/interpolation_size.py @@ -13,6 +13,7 @@ ELEMENTS_PER_EQUATION = 25 MAX_SURFACES = 10 SURFACE_STEP = 0.1 +MAX_EXTENT_FACTOR = 2.0 @dataclass @@ -20,15 +21,15 @@ class DataSummary: """The measures of the data of one feature. `spread` is 0 for parallel planes and 1 for orientations in all - directions. `spacing` is the median distance from a data point to the - nearest other data point (0 if there are fewer than two points). + directions. `area` is the map (XY) area of the box around the data points + (0 if it is not known). """ n_value: int = 0 n_orientation: int = 0 n_surfaces: int = 0 spread: float = 0.0 - spacing: float = 0.0 + area: float = 0.0 def _normals_from_strike_dip(strike, dip) -> np.ndarray: @@ -67,21 +68,14 @@ def orientation_spread(normals) -> float: return float(np.clip((1.0 - largest) / (2.0 / 3.0), 0.0, 1.0)) -def point_spacing(points) -> float: - """Return the median distance from each point to the nearest other point. - - Points at the same location count as one point. Returns 0 for fewer than - two points. - """ - from scipy.spatial import cKDTree - - points = np.asarray(points, dtype=float).reshape(-1, 3) +def xy_area(points) -> float: + """Return the map (XY) area of the box around the points, or 0.""" + points = np.asarray(points, dtype=float).reshape(-1, 2) points = points[np.all(np.isfinite(points), axis=1)] - points = np.unique(points, axis=0) if len(points) < 2: return 0.0 - distances, _ = cKDTree(points).query(points, k=2) - return float(np.median(distances[:, 1])) + width, height = points.max(axis=0) - points.min(axis=0) + return float(width * height) def _has(df, columns) -> np.ndarray: @@ -128,32 +122,33 @@ def summarise_data(df) -> DataSummary: n_surfaces = int(df.loc[_has(df, ['val']), 'val'].nunique()) if 'interface' in df.columns and df['interface'].notna().any(): n_surfaces = max(n_surfaces, int(df['interface'].nunique())) - spacing = 0.0 - if all(column in df.columns for column in ('X', 'Y', 'Z')): - spacing = point_spacing(df[['X', 'Y', 'Z']].to_numpy(dtype=float)) + area = 0.0 + if 'X' in df.columns and 'Y' in df.columns: + area = xy_area(df[['X', 'Y']].to_numpy(dtype=float)) return DataSummary( n_value=int(value.sum()), n_orientation=int(orientation.sum()), n_surfaces=n_surfaces, spread=spread, - spacing=spacing, + area=area, ) -def suggest_nelements(summary: DataSummary, volume: float = 0.0) -> int: +def suggest_nelements(summary: DataSummary, model_area: float = 0.0) -> int: """Return the number of elements for a feature with the data `summary`. - From the number of data rows: elements = equations x 25 x surface factor - x spread factor, where an orientation is two equations, the surface - factor is 1 + 0.1 for each surface after the first (at most 10 surfaces) - and the spread factor is 1 + spread. + elements = equations x 25 x surface factor x spread factor x extent + factor, where an orientation is two equations, the surface factor is + 1 + 0.1 for each surface after the first (at most 10 surfaces) and the + spread factor is 1 + spread. - From the model extent: if `volume` (the volume of the model box) and the - data spacing are known, elements = volume / spacing^3, so that an element - is not larger than the distance between the data points. + The extent factor is the map area of the model (`model_area`) divided by + the map area of the data, limited to 1 .. 2. When the data cover only a + part of the model, the elements are spread over a larger area, so more + elements keep the resolution near the data. The limit keeps the solve + time low. The factor is 1 if an area is not known. - The result is the larger of the two, rounded to 1 000 and limited to - 5 000 .. 250 000. + The result is rounded to 1 000 and limited to 5 000 .. 250 000. """ equations = summary.n_value + 2 * summary.n_orientation if equations <= 0: @@ -162,7 +157,7 @@ def suggest_nelements(summary: DataSummary, volume: float = 0.0) -> int: surface_factor = 1 + SURFACE_STEP * (surfaces - 1) spread_factor = 1 + min(max(summary.spread, 0.0), 1.0) elements = equations * ELEMENTS_PER_EQUATION * surface_factor * spread_factor - if volume > 0 and summary.spacing > 0: - elements = max(elements, min(volume / summary.spacing**3, MAX_NELEMENTS)) + if model_area > 0 and summary.area > 0: + elements *= min(max(model_area / summary.area, 1.0), MAX_EXTENT_FACTOR) elements = int(round(elements / ROUND_TO)) * ROUND_TO return int(min(max(elements, MIN_NELEMENTS), MAX_NELEMENTS)) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index d7bc1f7..543b647 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -1033,19 +1033,19 @@ def default_nelements(self, name) -> tuple: if not settings.interpolator_nelements_auto: return int(settings.interpolator_nelements), False summary = self._data_summaries.get(name, DataSummary()) - return suggest_nelements(summary, self._model_volume()), True + return suggest_nelements(summary, self._model_area()), True - def _model_volume(self) -> float: - """Return the volume of the model box, or 0 if the box is not set.""" + def _model_area(self) -> float: + """Return the map (XY) area of the model box, or 0 if the box is not set.""" bounding_box = getattr(self.model, 'bounding_box', None) try: - lengths = np.asarray(bounding_box.maximum, dtype=float) - np.asarray( + lengths = np.asarray(bounding_box.maximum, dtype=float)[:2] - np.asarray( bounding_box.origin, dtype=float - ) - except (AttributeError, TypeError, ValueError): + )[:2] + except (AttributeError, TypeError, ValueError, IndexError): return 0.0 - volume = float(np.prod(lengths)) - return volume if np.isfinite(volume) and volume > 0 else 0.0 + area = float(np.prod(lengths)) + return area if np.isfinite(area) and area > 0 else 0.0 def nelements_overrides_to_dict(self) -> dict: """Return `nelements_overrides` in a form that `json.dump` can write.""" diff --git a/tests/qgis/test_interpolator_arguments.py b/tests/qgis/test_interpolator_arguments.py index 45d5b56..6e17ba3 100644 --- a/tests/qgis/test_interpolator_arguments.py +++ b/tests/qgis/test_interpolator_arguments.py @@ -53,7 +53,7 @@ def test_the_number_of_the_user_comes_first(manager, monkeypatch): def test_automatic_uses_the_data(manager, monkeypatch): _use_settings(monkeypatch, interpolator_nelements_auto=True, interpolator_nelements=7000) - expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0 * 100.0) + expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0) assert manager.interpolator_arguments('s1', DATA)['nelements'] == expected assert manager.nelements_used['s1'] == (expected, True) @@ -68,7 +68,7 @@ def test_the_default_is_not_the_number_of_the_user(manager, monkeypatch): _use_settings(monkeypatch, interpolator_nelements_auto=True, interpolator_nelements=7000) manager.set_nelements_override('s1', 12345) manager.interpolator_arguments('s1', DATA) - expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0 * 100.0) + expected = suggest_nelements(summarise_data(DATA), 100.0 * 100.0) assert manager.default_nelements('s1') == (expected, True) diff --git a/tests/unit/test_interpolation_size.py b/tests/unit/test_interpolation_size.py index 08ff0e2..579322d 100644 --- a/tests/unit/test_interpolation_size.py +++ b/tests/unit/test_interpolation_size.py @@ -12,7 +12,7 @@ DataSummary, _normals_from_strike_dip, orientation_spread, - point_spacing, + xy_area, suggest_nelements, summarise_data, ) @@ -120,27 +120,36 @@ def test_the_normal_points_in_the_dip_direction(self): class TestExtent: - def test_point_spacing(self): - points = [(x, 0.0, 0.0) for x in (0.0, 10.0, 20.0, 30.0)] - assert point_spacing(points) == 10.0 + def test_xy_area(self): + assert xy_area([(0.0, 0.0), (10.0, 5.0), (3.0, 2.0)]) == 50.0 + assert xy_area([(1.0, 1.0)]) == 0.0 - def test_repeated_points_count_once(self): - assert point_spacing([(0.0, 0.0, 0.0)] * 3 + [(5.0, 0.0, 0.0)]) == 5.0 - assert point_spacing([(0.0, 0.0, 0.0)] * 3) == 0.0 + def test_summary_has_the_area(self): + df = pd.DataFrame({'X': [0.0, 50.0, 100.0], 'Y': [0.0, 20.0, 40.0], 'Z': 0.0, 'val': 0.0}) + assert summarise_data(df).area == 4000.0 - def test_summary_has_the_spacing(self): - df = pd.DataFrame({'X': [0.0, 50.0, 100.0], 'Y': 0.0, 'Z': 0.0, 'val': 0.0}) - assert summarise_data(df).spacing == 50.0 - - def test_a_large_model_with_dense_data_gives_more_elements(self): - summary = DataSummary(n_value=20, spacing=100.0) + def test_data_in_part_of_the_model_gives_more_elements(self): + summary = DataSummary(n_value=2000, area=100.0 * 100.0) + from_rows = suggest_nelements(summary) + same = suggest_nelements(summary, model_area=100.0 * 100.0) + quarter = suggest_nelements(summary, model_area=200.0 * 200.0) + assert same == from_rows + # the factor is at most 2 + assert quarter == 2 * from_rows + + def test_dense_map_data_do_not_give_the_maximum(self): + # five contact lines sampled every 50 m on a 10 km map + x = np.arange(0.0, 10_000.0, 50.0) + df = pd.concat( + pd.DataFrame({'X': x, 'Y': y, 'Z': 0.0, 'val': k * 100.0}) + for k, y in enumerate(np.linspace(1_000.0, 9_000.0, 5)) + ) + summary = summarise_data(df) from_rows = suggest_nelements(summary) - small = suggest_nelements(summary, volume=1000.0**3) - large = suggest_nelements(summary, volume=10_000.0 * 10_000.0 * 2_000.0) - assert small == from_rows - assert large > from_rows - assert large <= MAX_NELEMENTS + with_extent = suggest_nelements(summary, model_area=10_000.0 * 10_000.0) + assert from_rows <= with_extent <= 2 * from_rows + assert with_extent < MAX_NELEMENTS - def test_no_spacing_uses_only_the_rows(self): + def test_no_area_uses_only_the_rows(self): summary = DataSummary(n_value=20) - assert suggest_nelements(summary, volume=1e12) == suggest_nelements(summary) + assert suggest_nelements(summary, model_area=1e12) == suggest_nelements(summary)