feat: choose the number of elements from the data (phase 6.4)

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.

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)

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)

loopstructural/gui/dlg_settings.ui

@@ -280,6 +280,16 @@
         </property>
        </widget>
       </item>
+      <item row="4" column="1">
+       <widget class="QCheckBox" name="n_elements_auto_check_box">
+        <property name="toolTip">
+         <string>Choose the number of elements of each feature from its data. If this is off, every feature uses the fixed number.</string>
+        </property>
+        <property name="text">
+         <string>Automatic</string>
+        </property>
+       </widget>
+      </item>
      </layout>
     </widget>
    </item>

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:

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()