EX: new kernel based concentration analysis example.

examples/kernel_interpolation.py

@@ -0,0 +1,116 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import skimage
+
+import darsia
+
+# ! ---- DATA MANAGEMENT ---- !
+
+# The example images originate from a water tracer experiment with a multi-colored
+# indicator. Green color corresponds to concentrations close to 100%, while blue
+# color corresponds to concentrations ~99-0%. The images have been cropped and resized
+# mainly for data storage reasons. It is recommended to use full resolution images.
+baseline_path = (
+    f"{os.path.dirname(__file__)}/images/kernel_interpolation_example_base.npz"
+)
+image_path = f"{os.path.dirname(__file__)}/images/kernel_interpolation_example_test.npz"
+baseline = darsia.imread(baseline_path)
+image = darsia.imread(image_path)
+
+# ! ---- MAIN CONCENTRATION ANALYSIS AND CALIBRATION ---- !
+"""
+Before using kernel interpolation to transform the given image to a concentration image,
+the support data on which the interpolation depends has to be defined.
+In the thesis these were predefined samples of the smoothed and preprocessed image and
+corresponding concentrations.
+Now it is possible to select these samples in a GUI using darsia.BoxSelectionAssistant().
+In these samples of the image, the most common colour was identified using a histogram
+analysis, which is now available through darsia.extract_characteristic_data().
+"""
+
+# Predefine concentration analysis for now without any model (to be defined later).
+analysis = darsia.ConcentrationAnalysis(
+    base=baseline.img_as(float),
+    restoration=darsia.TVD(
+        weight=0.025, eps=1e-4, max_num_iter=100, method="isotropic Bregman"
+    ),
+    **{"diff option": "plain"},
+)
+
+# The goal is to define ome ROIs for which physical information is known.
+# One possibility is to use a GUI for interactive use. This option can
+# be activated on demand. For testing purposes this example by default
+# uses a pre-defined sample selection.
+if True:
+    samples = [
+        (slice(15, 40), slice(20, 45)),
+        (slice(15, 40), slice(220, 245)),
+        (slice(15, 40), slice(420, 445)),
+        (slice(15, 40), slice(720, 745)),
+    ]
+else:
+    # Same but under the use of a graphical user interface.
+    # Ask user to provide characteristic regions with expected concentration values
+    assistant = darsia.BoxSelectionAssistant(image)
+    samples = assistant()
+
+# For simplicity the concentrations are pre-defined. These could also be defined
+# by the user. Here, we use concentrations in %, such that 100 means fully saturated.
+n = len(samples)
+concentrations = 100 * np.append(np.linspace(1, 0.99, n - 1), 0)
+
+# Now add kernel interpolation as model trained by the extracted information.
+"""
+comments:
+-   Here, the support points as defined above are used to build a Kernel Interpolation.
+    For this, the 'concentration' without any model is used, correpsonding to the difference
+    to the baseline.
+-   This Kernel interpolation is then used as the model in the darsia.ConcentrationAnalysis
+-   use plain difference to keep all information (no cut off or norm)
+    this is the advantage of the kernel interpolation - it can handle negative colours
+-   finding the right kernel parameters is part of the modeling
+-   also clip the interpolated values at 100%
+"""
+smooth_RGB = analysis(image.img_as(float)).img
+colours_RGB = darsia.extract_characteristic_data(
+    signal=smooth_RGB, samples=samples, verbosity=True, surpress_plot=True
+)
+analysis.model = darsia.CombinedModel(
+    [
+        darsia.KernelInterpolation(
+            darsia.GaussianKernel(gamma=9.73), colours_RGB, concentrations
+        ),
+        darsia.ClipModel(**{"max value": 100.0}),
+    ]
+)
+
+# Finally, apply the (full) concentration analysis to analyze the test image
+concentration_image = analysis(image.img_as(float)).img
+
+# ! ----- VISUALISATION ---- !
+
+plt.figure("concentration average along one dimension")
+plt.plot(np.average(concentration_image, axis=0))
+plt.xlabel("horizontal pixel")
+plt.ylabel("concentration")
+
+fig = plt.figure()
+fig.suptitle("original image and resulting concentration")
+ax = plt.subplot(212)
+im = ax.imshow(concentration_image)
+ax.set_ylabel("vertical pixel")
+ax.set_xlabel("horizontal pixel")
+cb = fig.colorbar(im, orientation="horizontal", label="concentration [%]")
+ax = plt.subplot(211)
+ax.imshow(skimage.img_as_ubyte(image.img))
+ax.set_ylabel("vertical pixel")
+ax.set_xlabel("horizontal pixel")
+
+# To enable the example as test, the plots are closed after short time.
+# Pause longer if it is desired to keep the images on the screen.
+print("Warning: The plot is closed after short time to enable testing.")
+plt.show(block=False)
+plt.pause(5)
+plt.close()

src/darsia/__init__.py

@@ -27,6 +27,8 @@
 from darsia.utils.dtype import *
 from darsia.utils.grid import *
 from darsia.utils.fv import *
+from darsia.utils.kernels import *
+from darsia.utils.extractcharacteristicdata import *
 from darsia.corrections.basecorrection import *
 from darsia.corrections.shape.curvature import *
 from darsia.corrections.shape.affine import *
@@ -46,6 +48,7 @@
 from darsia.signals.models.staticthresholdmodel import *
 from darsia.signals.models.dynamicthresholdmodel import *
 from darsia.signals.models.binarydataselector import *
+from darsia.signals.models.kernelinterpolation import *
 from darsia.signals.reduction.signalreduction import *
 from darsia.signals.reduction.monochromatic import *
 from darsia.signals.reduction.dimensionreduction import *
@@ -68,6 +71,7 @@
 from darsia.manager.co2analysis import *
 from darsia.assistants.base_assistant import *
 from darsia.assistants.point_selection_assistant import *
+from darsia.assistants.box_selection_assistant import *
 from darsia.assistants.rotation_correction_assistant import *
 from darsia.assistants.subregion_assistant import *
 from darsia.assistants.crop_assistant import *

src/darsia/assistants/box_selection_assistant.py

@@ -0,0 +1,62 @@
+"""Module for defining subregions interactively."""
+
+from typing import Optional
+
+import numpy as np
+
+import darsia
+
+
+class BoxSelectionAssistant(darsia.PointSelectionAssistant):
+    def __init__(self, img: darsia.Image, **kwargs) -> None:
+        super().__init__(img)
+
+        self.name = "Box selection assistant"
+        """Name of assistant."""
+
+        self.width = kwargs.get("width", 100)
+        """Width and height of selected box(es)."""
+
+        self.boxes = None
+        """Output of assistant."""
+
+        self.shape = img.shape
+        """Image shape."""
+
+        self.dim = img.space_dim
+        """Dimension of image."""
+
+        if self.dim != 2:
+            raise NotImplementedError
+
+    def _convert_pts(self) -> None:
+        """Convert list of points to list of boxes in terms of slices."""
+
+        self.boxes = []
+        half_width = self.width / 2
+        # Point selection assistant produces points with cartesian ordering
+        pts_matrix_indexing = [np.flip(pt) for pt in self.pts]
+        for pt in pts_matrix_indexing:
+            self.boxes.append(
+                tuple(
+                    [
+                        slice(
+                            max(int(pt[d] - half_width), 0),
+                            min(int(pt[d] + half_width), self.shape[d]),
+                        )
+                        for d in range(self.dim)
+                    ]
+                )
+            )
+
+    def __call__(self) -> Optional[np.ndarray]:
+        """Call the assistant."""
+
+        super().__call__()
+        self._convert_pts()
+        return self.boxes
+
+    def _print_info(self) -> None:
+        """Print info about boxes so far assigned by the assistant."""
+        self._convert_pts()
+        print(self.boxes)

src/darsia/image/image.py

@@ -421,6 +421,36 @@ def astype(self, data_type) -> Any:
 
         return copy_image
 
+    def img_as(self, data_type) -> Any:
+        """Change data type via skimage.
+
+        Args:
+            data_type: target data type
+
+        Returns:
+            Image: image with transformed data type
+
+        """
+        copy_image = self.copy()
+        if data_type in [bool]:
+            copy_image.img = skimage.img_as_bool(copy_image.img)
+        elif data_type in [float]:
+            copy_image.img = skimage.img_as_float(copy_image.img)
+        elif data_type in [np.float32]:
+            copy_image.img = skimage.img_as_float32(copy_image.img)
+        elif data_type in [np.float64]:
+            copy_image.img = skimage.img_as_float64(copy_image.img)
+        elif data_type in [int]:
+            copy_image.img = skimage.img_as_int(copy_image.img)
+        elif data_type in [np.uint8]:
+            copy_image.img = skimage.img_as_ubyte(copy_image.img)
+        elif data_type in [np.uint16]:
+            copy_image.img = skimage.img_as_uint(copy_image.img)
+        else:
+            raise NotImplementedError
+
+        return copy_image
+
     # ! ---- Extraction routines
 
     def metadata(self) -> dict:

src/darsia/multi_image_analysis/concentrationanalysis.py

@@ -151,6 +151,8 @@ def find_cleaning_filter(
 
             # Use internal baseline images, if available.
             baseline_images = self._base_collection[1:]
+            if len(baseline_images) == 0:
+                baseline_images = None
 
         self.threshold_cleaning_filter = None
         """Cleaning filter."""
@@ -171,7 +173,7 @@ def find_cleaning_filter(
                 diff = self._subtract_background(probe_img)
 
                 # Extract scalar version
-                monochromatic_diff = self._extract_scalar_information(diff)
+                monochromatic_diff = self._reduce_signal(diff)
 
                 # Consider elementwise max
                 self.threshold_cleaning_filter = np.maximum(
@@ -230,7 +232,7 @@ def __call__(self, img: darsia.Image) -> darsia.Image:
         self._inspect_diff(diff)
 
         # Extract monochromatic version and take difference wrt the baseline image
-        signal = self._extract_scalar_information(diff)
+        signal = self._reduce_signal(diff)
 
         # Provide possibility for tuning and inspection of intermediate results
         self._inspect_scalar_signal(signal)
@@ -255,7 +257,11 @@ def __call__(self, img: darsia.Image) -> darsia.Image:
             plt.show()
 
         metadata = img.metadata()
-        return darsia.ScalarImage(concentration, **metadata)
+        is_scalar = len(concentration.shape) == len(img.shape) - 1
+        if is_scalar:
+            return darsia.ScalarImage(concentration, **metadata)
+        else:
+            return type(img)(concentration, **metadata)
 
     # ! ---- Inspection routines
     def _inspect_diff(self, img: np.ndarray) -> None:
@@ -318,6 +324,8 @@ def _subtract_background(self, img: darsia.Image) -> darsia.Image:
                 diff = np.clip(-img.img, 0, None)
             elif self._diff_option == "absolute":
                 diff = np.absolute(img.img)
+            elif self._diff_option == "plain":
+                diff = img.img
             else:
                 raise ValueError(f"Diff option {self._diff_option} not supported")
         else:
@@ -330,12 +338,14 @@ def _subtract_background(self, img: darsia.Image) -> darsia.Image:
                 diff = skimage.util.compare_images(
                     img.img, self.base.img, method="diff"
                 )
+            elif self._diff_option == "plain":
+                diff = img.img - self.base.img
             else:
                 raise ValueError(f"Diff option {self._diff_option} not supported")
 
         return diff
 
-    def _extract_scalar_information(self, img: np.ndarray) -> np.ndarray:
+    def _reduce_signal(self, img: np.ndarray) -> np.ndarray:
         """
         Make a scalar image from potentially multi-colored image.
 
@@ -382,10 +392,7 @@ def _balance_signal(self, img: np.ndarray) -> np.ndarray:
             np.ndarray: balanced image
 
         """
-        if self.balancing is not None:
-            return self.balancing(img)
-        else:
-            return img
+        return img if self.balancing is None else self.balancing(img)
 
     def _prepare_signal(self, signal: np.ndarray) -> np.ndarray:
         """

src/darsia/restoration/tvd.py

@@ -28,7 +28,7 @@ def __init__(self, key: str = "", **kwargs) -> None:
 
         """
         # Determine method
-        self.method = kwargs.pop(key + "method", "chambolle")
+        self.method = kwargs.pop(key + "method", "chambolle").lower()
 
         # Method-specific parameters
         if self.method == "heterogeneous bregman":