Implement LoopStructural interpolation algorithms

- Complete processAlgorithm for LoopStructuralInterpolationAlgorithm
- Create new LoopStructuralInterpolationEvaluationAlgorithm
- Register new algorithm in provider
- Add support for building interpolators from value and gradient constraints
- Add support for evaluating interpolators on raster grids, 3D grids, and point sets

Co-authored-by: lachlangrose <7371904+lachlangrose@users.noreply.github.com>

Author: Copilot

loopstructural/processing/algorithms/interpolation/interpolation_algorithm.py

@@ -10,22 +10,26 @@
 """
 
 from typing import Any, Optional
+import tempfile
+import numpy as np
+import dill as pickle
 
-from qgis import processing
 from qgis.core import (
-    QgsFeatureSink,
     QgsProcessing,
     QgsProcessingAlgorithm,
     QgsProcessingContext,
     QgsProcessingException,
     QgsProcessingFeedback,
     QgsProcessingParameterExtent,
     QgsProcessingParameterNumber,
-    QgsProcessingParameterFeatureSink,
     QgsProcessingParameterFeatureSource,
-    QgsProcessingParameterField
+    QgsProcessingParameterField,
+    QgsProcessingParameterFileDestination,
 )
 
+from LoopStructural import InterpolatorBuilder
+from LoopStructural.datatypes import BoundingBox
+
 
 class LoopStructuralInterpolationAlgorithm(QgsProcessingAlgorithm):
     """Processing algorithm to create basal contacts."""
@@ -146,33 +150,217 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
         )
 
         self.addParameter(
-            QgsProcessingParameterFeatureSink(
+            QgsProcessingParameterFileDestination(
                 self.OUTPUT,
-                "Interpolated Surface",
+                "Output interpolator file (pickle)",
+                fileFilter='Pickle files (*.pkl)'
             )
         )
-        pass
 
     def processAlgorithm(
         self,
         parameters: dict[str, Any],
         context: QgsProcessingContext,
         feedback: QgsProcessingFeedback,
     ) -> dict[str, Any]:
+        # Get input parameters
         value_layer = self.parameterAsSource(parameters, self.VALUE, context)
         value_field = self.parameterAsFields(parameters, self.VALUE_FIELD, context)
 
         gradient_layer = self.parameterAsSource(parameters, self.GRADIENT, context)
         strike_field = self.parameterAsFields(parameters, self.STRIKE_FIELD, context)
         dip_field = self.parameterAsFields(parameters, self.DIP_FIELD, context)
-
-        inequality_layer = self.parameterAsSource(parameters, self.INEQUALITY, context)
-        upper_field = self.parameterAsFields(parameters, self.INEQUALITY_UPPER_FIELD, context)
-        lower_field = self.parameterAsFields(parameters, self.INEQUALITY_LOWER_FIELD, context)
 
         extent = self.parameterAsExtent(parameters, self.EXTENT, context)
+        pixel_size = self.parameterAsDouble(parameters, self.PIXEL_SIZE, context)
+        
+        # Validate that we have at least some data
+        if value_layer is None and gradient_layer is None:
+            raise QgsProcessingException(
+                "At least one of value or gradient layers must be provided."
+            )
+        
+        # Extract field names from lists
+        value_field_name = value_field[0] if value_field else None
+        strike_field_name = strike_field[0] if strike_field else None
+        dip_field_name = dip_field[0] if dip_field else None
+        
+        feedback.pushInfo("Building interpolator from constraints...")
+        
+        # Create bounding box from extent
+        if extent is None:
+            raise QgsProcessingException("Output extent must be provided.")
+        
+        # Calculate z extent based on pixel size (arbitrary choice: 10x pixel size for depth)
+        z_min = 0
+        z_max = pixel_size * 10
+        
+        bbox = BoundingBox(
+            origin=[extent.xMinimum(), extent.yMinimum(), z_min],
+            maximum=[extent.xMaximum(), extent.yMaximum(), z_max]
+        )
+        
+        # Create interpolator builder
+        interpolator_builder = InterpolatorBuilder(
+            interpolatortype='PLI',  # Piecewise Linear Interpolator
+            bounding_box=bbox
+        )
+        
+        # Add value constraints
+        if value_layer is not None and value_field_name:
+            feedback.pushInfo(f"Adding value constraints from field '{value_field_name}'...")
+            value_constraints = self._extract_value_constraints(
+                value_layer, value_field_name, feedback
+            )
+            if len(value_constraints) > 0:
+                interpolator_builder.add_value_constraints(value_constraints)
+                feedback.pushInfo(f"Added {len(value_constraints)} value constraints.")
+            else:
+                feedback.pushWarning("No valid value constraints found.")
+        
+        # Add gradient constraints from strike/dip
+        if gradient_layer is not None and strike_field_name and dip_field_name:
+            feedback.pushInfo(
+                f"Adding gradient constraints from strike/dip fields '{strike_field_name}'/'{dip_field_name}'..."
+            )
+            gradient_constraints = self._extract_gradient_constraints(
+                gradient_layer, strike_field_name, dip_field_name, feedback
+            )
+            if len(gradient_constraints) > 0:
+                interpolator_builder.add_gradient_constraints(gradient_constraints)
+                feedback.pushInfo(f"Added {len(gradient_constraints)} gradient constraints.")
+            else:
+                feedback.pushWarning("No valid gradient constraints found.")
+        
+        # Build the interpolator
+        feedback.pushInfo("Building interpolator...")
+        try:
+            interpolator = interpolator_builder.build()
+        except Exception as e:
+            raise QgsProcessingException(f"Failed to build interpolator: {e}")
+        
+        # Save interpolator to pickle file
+        out_path = self.parameterAsString(parameters, self.OUTPUT, context)
+        if not out_path:
+            tmp = tempfile.NamedTemporaryFile(suffix='.pkl', delete=False)
+            out_path = tmp.name
+            tmp.close()
+        
+        feedback.pushInfo(f"Saving interpolator to {out_path}...")
+        try:
+            with open(out_path, 'wb') as fh:
+                pickle.dump(interpolator, fh)
+        except Exception as e:
+            raise QgsProcessingException(f"Failed to save interpolator to '{out_path}': {e}")
+        
+        feedback.pushInfo("Interpolator saved successfully.")
+        return {self.OUTPUT: out_path}
+    
+    def _extract_value_constraints(
+        self, 
+        source,
+        value_field: str,
+        feedback: QgsProcessingFeedback
+    ) -> np.ndarray:
+        """Extract value constraints as numpy array (x, y, z, value)."""
+        constraints = []
+        for feat in source.getFeatures():
+            geom = feat.geometry()
+            if geom is None or geom.isEmpty():
+                continue
+            
+            # Get value
+            try:
+                value = feat.attribute(value_field)
+                if value is None:
+                    continue
+                value = float(value)
+            except (ValueError, TypeError):
+                continue
+            
+            # Extract points from geometry
+            points = []
+            if geom.isMultipart():
+                if geom.type() == 0:  # Point
+                    points = geom.asMultiPoint()
+                elif geom.type() == 1:  # Line
+                    for line in geom.asMultiPolyline():
+                        points.extend(line)
+            else:
+                if geom.type() == 0:  # Point
+                    points = [geom.asPoint()]
+                elif geom.type() == 1:  # Line
+                    points = geom.asPolyline()
+            
+            # Add constraint for each point
+            for p in points:
+                try:
+                    x = p.x()
+                    y = p.y()
+                    z = p.z() if hasattr(p, 'z') else 0.0
+                    constraints.append([x, y, z, value])
+                except Exception:
+                    continue
+        
+        return np.array(constraints) if constraints else np.array([]).reshape(0, 4)
+    
+    def _extract_gradient_constraints(
+        self,
+        source,
+        strike_field: str,
+        dip_field: str,
+        feedback: QgsProcessingFeedback
+    ) -> np.ndarray:
+        """Extract gradient constraints from strike/dip as numpy array (x, y, z, gx, gy, gz)."""
+        constraints = []
+        for feat in source.getFeatures():
+            geom = feat.geometry()
+            if geom is None or geom.isEmpty():
+                continue
+            
+            # Get strike and dip
+            try:
+                strike = feat.attribute(strike_field)
+                dip = feat.attribute(dip_field)
+                if strike is None or dip is None:
+                    continue
+                strike = float(strike)
+                dip = float(dip)
+            except (ValueError, TypeError):
+                continue
+            
+            # Extract points from geometry (should be point geometry)
+            points = []
+            if geom.isMultipart():
+                if geom.type() == 0:  # Point
+                    points = geom.asMultiPoint()
+            else:
+                if geom.type() == 0:  # Point
+                    points = [geom.asPoint()]
+            
+            # Convert strike/dip to gradient vector
+            # Strike is azimuth (0-360), dip is angle from horizontal (0-90)
+            # This follows geological convention
+            strike_rad = np.radians(strike)
+            dip_rad = np.radians(dip)
+            
+            # Gradient vector (normal to the plane)
+            # This is the gradient direction pointing updip
+            gx = np.sin(dip_rad) * np.sin(strike_rad)
+            gy = -np.sin(dip_rad) * np.cos(strike_rad)
+            gz = np.cos(dip_rad)
+            
+            # Add constraint for each point
+            for p in points:
+                try:
+                    x = p.x()
+                    y = p.y()
+                    z = p.z() if hasattr(p, 'z') else 0.0
+                    constraints.append([x, y, z, gx, gy, gz])
+                except Exception:
+                    continue
 
-        print(extent)
+        return np.array(constraints) if constraints else np.array([]).reshape(0, 6)
 
     def createInstance(self) -> QgsProcessingAlgorithm:
         """Create a new instance of the algorithm."""