From 58757ad5b9a2a8b96b97caff80b8589893ad8331 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Mon, 18 Nov 2024 16:38:22 +0000
Subject: [PATCH 01/34] WIP: Scope out refactor for grains.py

---
 topostats/grains.py       | 342 +++++++++++++++++++++++++++-----------
 topostats/processing.py   |  18 +-
 topostats/unet_masking.py |  84 ++++++++++
 3 files changed, 339 insertions(+), 105 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 159c8f6e8e..ba382e21cb 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -6,6 +6,7 @@
 import logging
 import sys
 from collections import defaultdict
+from dataclasses import dataclass, field
 
 import keras
 import numpy as np
@@ -21,6 +22,8 @@
 from topostats.unet_masking import (
     make_bounding_box_square,
     pad_bounding_box,
+    pad_crop,
+    make_crop_square,
     predict_unet,
 )
 from topostats.utils import _get_mask, get_thresholds
@@ -37,6 +40,102 @@
 # pylint: disable=too-many-public-methods
 
 
+@dataclass
+class GrainCrop:
+    """
+    Dataclass for storing the crops of grains.
+
+    Attributes
+    ----------
+    image : npt.NDArray[np.float32]
+        2-D Numpy array of the cropped image.
+    mask : npt.NDArray[np.bool_]
+        2-D Numpy array of the cropped mask.
+    padding : int
+        Padding added to the bounding box of the grain during cropping.
+    bbox: tuple[int, int, int, int]
+        Bounding box of the crop including padding.
+    """
+
+    image: npt.NDArray[np.float32]
+    mask: npt.NDArray[np.bool_]
+    padding: int
+    bbox: tuple[int, int, int, int]
+
+
+def validate_full_mask_tensor_shape(array: npt.NDArray[np.bool_]) -> npt.NDArray[np.bool_]:
+    """
+    Validate the shape of the full mask tensor.
+
+    Parameters
+    ----------
+    array : npt.NDArray
+        Numpy array to validate.
+
+    Returns
+    -------
+    npt.NDArray
+        Numpy array if valid.
+    """
+    if len(array.shape) != 3 or array.shape[2] != 3 or array.shape[1] != array.shape[0]:
+        raise ValueError(f"Full mask tensor must be NxNx3 but has shape {array.shape}")
+    return array
+
+
+@dataclass
+class GrainCropsDirection:
+    """
+    Dataclass for storing the crops of grains in a particular imaging direction.
+
+    Attributes
+    ----------
+    full_mask_tensor : npt.NDArray[np.bool_]
+        Boolean NxNx3 array of the full mask tensor.
+    crops : GrainCrops
+        Grain crops.
+    """
+
+    full_mask_tensor: npt.NDArray[np.bool_]
+    crops: dict[int, GrainCrop]
+
+    def __post_init__(self):
+        """
+        Validate the full mask tensor shape.
+
+        Raises
+        ------
+        ValueError
+            If the full mask tensor shape is invalid.
+        """
+        self._full_mask_tensor = validate_full_mask_tensor_shape(self.full_mask_tensor)
+
+    @property
+    def full_mask_tensor(self) -> npt.NDArray[np.bool_]:
+        """Return the full mask tensor."""
+        return self._full_mask_tensor
+
+    @full_mask_tensor.setter
+    def full_mask_tensor(self, value: npt.NDArray[np.bool_]):
+        self._full_mask_tensor = validate_full_mask_tensor_shape(value)
+
+
+@dataclass
+class ImageGrainCrops:
+    """
+    Dataclass for storing the crops of grains in an image.
+
+    Attributes
+    ----------
+    above : GrainCropDirection | None
+        Grains in the above direction.
+    below : GrainCropDirection | None
+        Grains in the below direction.
+    """
+
+    above: GrainCropsDirection | None
+    below: GrainCropsDirection | None
+
+
 class Grains:
     """
     Find grains in an image.
@@ -98,6 +197,7 @@ def __init__(
         remove_edge_intersecting_grains: bool = True,
         classes_to_merge: list[tuple[int, int]] | None = None,
         vetting: dict | None = None,
+        grain_crop_padding: int = 1,
     ):
         """
         Initialise the class.
@@ -177,6 +277,7 @@ def __init__(
         self.region_properties = defaultdict()
         self.bounding_boxes = defaultdict()
         self.grainstats = None
+        self.grain_crop_padding = grain_crop_padding
         self.unet_config = unet_config
         self.vetting = vetting
         self.classes_to_merge = classes_to_merge
@@ -550,28 +651,31 @@ def find_grains(self):
                 np.int32
             )  # Will produce an NxNx2 array
 
+            full_mask_tensor = self.directions[direction]["labelled_regions_02"]
+
+            graincrops = self.extract_grains_from_full_image_mask(
+                image=self.image,
+                labelled_full_mask_tensor=self.directions[direction]["labelled_regions_02"],
+                padding=self.grain_crop_padding,
+            )
+
             # Check whether to run the UNet model
             if self.unet_config["model_path"] is not None:
 
                 # Run unet segmentation on only the class 1 layer of the labelled_regions_02. Need to make this configurable
                 # later on along with all the other hardcoded class 1s.
-                unet_mask, unet_labelled_regions = Grains.improve_grain_segmentation_unet(
+                graincrops = Grains.improve_grain_segmentation_unet(
                     filename=self.filename,
                     direction=direction,
                     unet_config=self.unet_config,
                     image=self.image,
-                    labelled_grain_regions=self.directions[direction]["labelled_regions_02"][:, :, 1],
+                    graincrops=graincrops,
                 )
 
-                # Update the image masks to be the unet masks instead
-                self.directions[direction]["removed_small_objects"] = unet_mask
-                self.directions[direction]["labelled_regions_02"] = unet_labelled_regions
-
-                class_counts = [
-                    unet_labelled_regions[class_idx].max() for class_idx in range(unet_labelled_regions.shape[2])
-                ]
-                LOGGER.info(
-                    f"[{self.filename}] : Overridden {thresholding_grain_count} grains with {class_counts} UNet predictions ({direction})"
+                # Construct full mask from the crops
+                full_mask_tensor = Grains.construct_full_mask_from_crops(
+                    graincrops=graincrops,
+                    image_shape=self.image.shape,
                 )
 
             # Vet the grains
@@ -604,6 +708,28 @@ def find_grains(self):
             self.directions[direction]["removed_small_objects"] = labelled_final_grains.astype(bool)
             self.directions[direction]["labelled_regions_02"] = labelled_final_grains.astype(np.int32)
 
+            self.directions[direction]["grain_crops_direction"] = GrainCropsDirection(
+                full_mask_tensor=full_mask_tensor,
+                crops=graincrops,
+            )
+
+        if "above" in self.directions and "below" in self.directions:
+            return ImageGrainCrops(
+                above=self.directions["above"]["grain_crops_direction"],
+                below=self.directions["below"]["grain_crops_direction"],
+            )
+        if "above" in self.directions:
+            return ImageGrainCrops(
+                above=self.directions["above"]["grain_crops_direction"],
+                below=None,
+            )
+        if "below" in self.directions:
+            return ImageGrainCrops(
+                above=None,
+                below=self.directions["below"]["grain_crops_direction"],
+            )
+        raise ValueError("No grains found in either direction")
+
     # pylint: disable=too-many-locals
     @staticmethod
     def improve_grain_segmentation_unet(
@@ -611,8 +737,8 @@ def improve_grain_segmentation_unet(
         direction: str,
         unet_config: dict[str, str | int | float | tuple[int | None, int, int, int] | None],
         image: npt.NDArray,
-        labelled_grain_regions: npt.NDArray,
-    ) -> tuple[npt.NDArray, npt.NDArray]:
+        graincrops: dict[int, GrainCrop],
+    ) -> dict[int, GrainCrop]:
         """
         Use a UNet model to re-segment existing grains to improve their accuracy.
 
@@ -667,59 +793,16 @@ def improve_grain_segmentation_unet(
         # unet_model = keras.models.load_model(unet_config["model_path"], custom_objects={"mean_iou": mean_iou})
         LOGGER.debug(f"Output shape of UNet model: {unet_model.output_shape}")
 
-        # Initialise an empty mask to iteratively add to for each grain, with the correct number of class channels based on
-        # the loaded model's output shape
-        # Note that the minimum number of classes is 2, as even for binary outputs, we will force categorical type
-        # data, so we have a class for background.
-        unet_mask = np.zeros((image.shape[0], image.shape[1], np.max([2, unet_model.output_shape[-1]]))).astype(
-            np.bool_
-        )
-        # Set the background class to be all 1s by default since not all of the image will be covered by the
-        # u-net predictions.
-        unet_mask[:, :, 0] = 1
-        # Labelled regions will be the same by default, but will be overwritten if there are any grains present.
-        unet_labelled_regions = np.zeros_like(unet_mask).astype(np.int32)
-        # For each detected molecule, create an image of just that molecule and run the UNet
-        # on that image to segment it
-        grain_region_properties = regionprops(labelled_grain_regions)
-        for grain_number, region in enumerate(grain_region_properties):
-            LOGGER.debug(f"Unet predicting mask for grain {grain_number} of {len(grain_region_properties)}")
-
-            # Get the bounding box for the region
-            bounding_box: tuple[int, int, int, int] = tuple(region.bbox)  # min_row, min_col, max_row, max_col
-
-            # Pad the bounding box
-            bounding_box = pad_bounding_box(
-                crop_min_row=bounding_box[0],
-                crop_min_col=bounding_box[1],
-                crop_max_row=bounding_box[2],
-                crop_max_col=bounding_box[3],
-                image_shape=(image.shape[0], image.shape[1]),
-                padding=unet_config["grain_crop_padding"],
-            )
-
-            # Make the bounding box square within the confines of the image
-            bounding_box = make_bounding_box_square(
-                crop_min_row=bounding_box[0],
-                crop_min_col=bounding_box[1],
-                crop_max_row=bounding_box[2],
-                crop_max_col=bounding_box[3],
-                image_shape=(image.shape[0], image.shape[1]),
-            )
-
-            # Grab the cropped image. Using slice since the bounding box from skimage is
-            # half-open, so the max_row and max_col are not included in the region.
-            region_image = image[
-                bounding_box[0] : bounding_box[2],
-                bounding_box[1] : bounding_box[3],
-            ]
+        new_graincrops = {}
+        for grain_number, graincrop in graincrops.items():
+            LOGGER.debug(f"Unet predicting mask for grain {grain_number} of {len(graincrops)}")
 
             # Run the UNet on the region. This is allowed to be a single channel
             # as we can add a background channel afterwards if needed.
             # Remember that this region is cropped from the original image, so it's not
             # the same size as the original image.
             predicted_mask = predict_unet(
-                image=region_image,
+                image=graincrop.image,
                 model=unet_model,
                 confidence=0.1,
                 model_input_shape=unet_model.input_shape,
@@ -730,45 +813,16 @@ def improve_grain_segmentation_unet(
             assert len(predicted_mask.shape) == 3
             LOGGER.debug(f"Predicted mask shape: {predicted_mask.shape}")
 
-            # Add each class of the predicted mask to the overall full image mask
-            for class_index in range(unet_mask.shape[2]):
+            new_graincrops[grain_number] = GrainCrop(
+                image=graincrop.image,
+                mask=predicted_mask,
+                padding=graincrop.padding,
+                bbox=graincrop.bbox,
+            )
 
                 # Grab the unet mask for the class
                 unet_predicted_mask_labelled = morphology.label(predicted_mask[:, :, class_index])
-
-                # Directly set the background to be equal instead of logical or since they are by default
-                # 1, and should be turned off if any other class is on
-                if class_index == 0:
-                    unet_mask[
-                        bounding_box[0] : bounding_box[2],
-                        bounding_box[1] : bounding_box[3],
-                        class_index,
-                    ] = unet_predicted_mask_labelled
-                else:
-                    unet_mask[
-                        bounding_box[0] : bounding_box[2],
-                        bounding_box[1] : bounding_box[3],
-                        class_index,
-                    ] = np.logical_or(
-                        unet_mask[
-                            bounding_box[0] : bounding_box[2],
-                            bounding_box[1] : bounding_box[3],
-                            class_index,
-                        ],
-                        unet_predicted_mask_labelled,
-                    )
-
-            assert len(unet_mask.shape) == 3, f"Unet mask shape: {unet_mask.shape}"
-            assert unet_mask.shape[-1] >= 2, f"Unet mask shape: {unet_mask.shape}"
-
-            # For each class in the unet mask tensor, label the mask and add to unet_labelled_regions
-            # The labelled background class will be identical to the binary one from the unet mask.
-            unet_labelled_regions[:, :, 0] = unet_mask[:, :, 0]
-            # Iterate over each class and label the regions
-            for class_index in range(unet_mask.shape[2]):
-                unet_labelled_regions[:, :, class_index] = Grains.label_regions(unet_mask[:, :, class_index])
-
-        return unet_mask, unet_labelled_regions
+        return graincrops
 
     @staticmethod
     def keep_largest_labelled_region(
@@ -1573,3 +1627,99 @@ def merge_classes(
             grain_mask_tensor = np.dstack([grain_mask_tensor, combined_mask])
 
         return grain_mask_tensor.astype(bool)
+
+    @staticmethod
+    def construct_full_mask_from_crops(graincrops: dict[int, GrainCrop], image_shape: tuple[int, int]) -> npt.NDArray:
+        """
+        Construct a full mask tensor from the grain crops.
+
+        Parameters
+        ----------
+        graincrops : dict[int, GrainCrop]
+            Dictionary of grain crops.
+        image_shape : tuple[int, int]
+            Shape of the original image.
+
+        Returns
+        -------
+        npt.NDArray
+            NxNx3 Numpy array of the full mask tensor.
+        """
+        full_mask_tensor = np.zeros((image_shape[0], image_shape[1], 3), dtype=np.bool_)
+        for grain_number, graincrop in graincrops.items():
+            bounding_box = graincrop.bbox
+            crop = graincrop.mask
+            padding = graincrop.padding
+
+            # Add the crop to the full mask tensor without overriding anything else, for all classes
+            for class_index in range(crop.shape[2]):
+                full_mask_tensor[
+                    bounding_box[0] + padding : bounding_box[2] - padding,
+                    bounding_box[1] + padding : bounding_box[3] - padding,
+                    class_index,
+                ] += crop[:, :, class_index]
+
+        return full_mask_tensor
+
+    @staticmethod
+    def extract_grains_from_full_image_mask(
+        image: npt.NDArray[np.float32],
+        labelled_full_mask_tensor: npt.NDArray[np.bool_],
+        padding: int,
+    ) -> dict[int, GrainCrop]:
+        """
+        Extract grains from the full image mask tensor.
+
+        Parameters
+        ----------
+        image : npt.NDArray[np.float32]
+            2-D Numpy array of the image.
+        labelled_full_mask_tensor : npt.NDArray[np.bool_]
+            Labelled 3-D Numpy array of the full mask tensor.
+        padding : int
+            Padding added to the bounding box of the grain before cropping.
+
+        Returns
+        -------
+        dict[int, GrainCrop]
+            Dictionary of grain crops.
+        """
+        grain_region_properties = regionprops(labelled_full_mask_tensor)
+        graincrops = {}
+        for grain_number, region in enumerate(grain_region_properties):
+            # Get the bounding box for the region
+            bounding_box: tuple[int, int, int, int] = tuple(region.bbox)  # min_row, min_col, max_row, max_col
+            mask_crop = labelled_full_mask_tensor[
+                bounding_box[0] : bounding_box[2],
+                bounding_box[1] : bounding_box[3],
+            ]
+
+            # Pad the mask
+            padded_mask, bounding_box = pad_crop(
+                crop=mask_crop,
+                bbox=bounding_box,
+                image_shape=(image.shape[0], image.shape[1]),
+                padding=padding,
+            )
+
+            square_padded_mask, bounding_box = make_crop_square(
+                crop=padded_mask,
+                bbox=bounding_box,
+                image_shape=(image.shape[0], image.shape[1]),
+            )
+
+            # Grab the cropped image. Using slice since the bounding box from skimage is
+            # half-open, so the max_row and max_col are not included in the region.
+            region_image = image[
+                bounding_box[0] : bounding_box[2],
+                bounding_box[1] : bounding_box[3],
+            ]
+
+            graincrops[grain_number] = GrainCrop(
+                image=region_image,
+                mask=square_padded_mask,
+                padding=padding,
+                bbox=bounding_box,
+            )
+
+        return graincrops
diff --git a/topostats/processing.py b/topostats/processing.py
index 8d6388a25f..cd9153dbfc 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -833,7 +833,9 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
+                    splining_molstats = create_empty_dataframe(
+                        column_set="mol_statistics", index_col="molecule_number"
+                    )
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1014,7 +1016,7 @@ def process_scan(
     )
 
     # Find Grains :
-    grain_masks = run_grains(
+    grain_dict = run_grains(
         image=topostats_object["image_flattened"],
         pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
         filename=topostats_object["filename"],
@@ -1024,16 +1026,12 @@ def process_scan(
         grains_config=grains_config,
     )
 
-    # Update grain masks if new grain masks are returned. Else keep old grain masks. Topostats object's "grain_masks"
-    # defaults to an empty dictionary so this is safe.
-    topostats_object["grain_masks"] = grain_masks if grain_masks is not None else topostats_object["grain_masks"]
-
     if "above" in topostats_object["grain_masks"].keys() or "below" in topostats_object["grain_masks"].keys():
         # Grainstats :
         grainstats_df, height_profiles = run_grainstats(
             image=topostats_object["image_flattened"],
             pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
-            grain_masks=topostats_object["grain_masks"],
+            grain_dict=grain_dict,
             filename=topostats_object["filename"],
             basename=topostats_object["img_path"],
             grainstats_config=grainstats_config,
@@ -1045,7 +1043,7 @@ def process_scan(
         # Disordered Tracing
         disordered_traces_data, grainstats_df, disordered_tracing_stats = run_disordered_tracing(
             image=topostats_object["image_flattened"],
-            grain_masks=topostats_object["grain_masks"],
+            grain_dict=grain_dict,
             pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
             filename=topostats_object["filename"],
             basename=topostats_object["img_path"],
@@ -1104,7 +1102,9 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
+        disordered_tracing_stats = create_empty_dataframe(
+            column_set="disordered_tracing_statistics", index_col="index"
+        )
         height_profiles = {}
 
     # Get image statistics
diff --git a/topostats/unet_masking.py b/topostats/unet_masking.py
index a9734fe83b..01fd808cb3 100644
--- a/topostats/unet_masking.py
+++ b/topostats/unet_masking.py
@@ -325,3 +325,87 @@ def pad_bounding_box(
     new_crop_max_col: int = min(image_shape[1], crop_max_col + padding)
 
     return new_crop_min_row, new_crop_min_col, new_crop_max_row, new_crop_max_col
+
+
+def pad_crop(
+    crop: npt.NDArray,
+    bbox: tuple[int, int, int, int],
+    image_shape: tuple[int, int],
+    padding: int,
+) -> np.NDArray:
+    """
+    Pad a crop.
+
+    Parameters
+    ----------
+    crop : npt.NDArray
+        The crop to pad.
+    bbox : tuple[int, int, int, int]
+        The bounding box of the crop.
+    image_shape : tuple[int, int]
+        The shape of the image.
+    padding : int
+        The padding to apply to the crop.
+
+    Returns
+    -------
+    np.NDArray
+        The padded crop.
+    """
+    new_bounding_box = pad_bounding_box(
+        crop_min_row=bbox[0],
+        crop_min_col=bbox[1],
+        crop_max_row=bbox[2],
+        crop_max_col=bbox[3],
+        image_shape=image_shape,
+        padding=padding,
+    )
+
+    # Pad the crop with the new bounding box
+    padded_crop = np.zeros((new_bounding_box[2] - new_bounding_box[0], new_bounding_box[3] - new_bounding_box[1]))
+    padded_crop[
+        bbox[0] - new_bounding_box[0] : bbox[2] - new_bounding_box[0],
+        bbox[1] - new_bounding_box[1] : bbox[3] - new_bounding_box[1],
+    ] = crop
+
+    return padded_crop, new_bounding_box
+
+
+def make_crop_square(
+    crop: npt.NDArray,
+    bbox: tuple[int, int, int, int],
+    image_shape: tuple[int, int],
+) -> np.NDArray:
+    """
+    Make a crop square.
+
+    Parameters
+    ----------
+    crop : npt.NDArray
+        The crop to make square.
+    bbox : tuple[int, int, int, int]
+        The bounding box of the crop.
+    image_shape : tuple[int, int]
+        The shape of the image.
+
+    Returns
+    -------
+    np.NDArray
+        The square crop.
+    """
+    new_bounding_box = make_bounding_box_square(
+        crop_min_row=bbox[0],
+        crop_min_col=bbox[1],
+        crop_max_row=bbox[2],
+        crop_max_col=bbox[3],
+        image_shape=image_shape,
+    )
+
+    # Make the crop square
+    square_crop = np.zeros((new_bounding_box[2] - new_bounding_box[0], new_bounding_box[3] - new_bounding_box[1]))
+    square_crop[
+        bbox[0] - new_bounding_box[0] : bbox[2] - new_bounding_box[0],
+        bbox[1] - new_bounding_box[1] : bbox[3] - new_bounding_box[1],
+    ] = crop
+
+    return square_crop, new_bounding_box

From f49b1e6b0bf47dbe76408d54a62d0e554341a25d Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 19 Nov 2024 11:13:27 +0000
Subject: [PATCH 02/34] WIP: Begin grains > grainstats pipeline overhaul.
 Outline data structures.

---
 topostats/processing.py | 81 ++++++++++++++++++-----------------------
 1 file changed, 36 insertions(+), 45 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index cd9153dbfc..c234e6232a 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -13,7 +13,7 @@
 
 from topostats import __version__
 from topostats.filters import Filters
-from topostats.grains import Grains
+from topostats.grains import Grains, ImageGrainCrops, GrainCropsDirection
 from topostats.grainstats import GrainStats
 from topostats.io import get_out_path, save_topostats_file
 from topostats.logs.logs import LOGGER_NAME
@@ -135,7 +135,7 @@ def run_grains(  # noqa: C901
     core_out_path: Path,
     plotting_config: dict,
     grains_config: dict,
-) -> dict | None:
+) -> ImageGrainCrops | None:
     """
     Identify grains (molecules) and optionally plots the results.
 
@@ -250,7 +250,7 @@ def run_grains(  # noqa: C901
 def run_grainstats(
     image: npt.NDArray,
     pixel_to_nm_scaling: float,
-    grain_masks: dict,
+    image_grain_crops: ImageGrainCrops,
     filename: str,
     basename: Path,
     grainstats_config: dict,
@@ -300,49 +300,40 @@ def run_grainstats(
             }
             grainstats_dict = {}
             height_profiles_dict = {}
+
             # There are two layers to process those above the given threshold and those below
-            for direction, _ in grain_masks.items():
-                # Get the DNA class mask from the tensor
-                LOGGER.debug(f"[{filename}] : Full Mask dimensions: {grain_masks[direction].shape}")
-                assert len(grain_masks[direction].shape) == 3, "Grain masks should be 3D tensors"
-                dna_class_mask = grain_masks[direction][:, :, 1]
-                LOGGER.debug(f"[{filename}] : DNA Mask dimensions: {dna_class_mask.shape}")
-                # Check if there are grains
-                if np.max(dna_class_mask) == 0:
-                    LOGGER.warning(
-                        f"[{filename}] : No grains exist for the {direction} direction. Skipping grainstats for {direction}."
-                    )
-                    grainstats_dict[direction] = create_empty_dataframe(
-                        column_set="grainstats", index_col="grain_number"
-                    )
-                else:
-                    grainstats_calculator = GrainStats(
-                        data=image,
-                        labelled_data=dna_class_mask,
-                        pixel_to_nanometre_scaling=pixel_to_nm_scaling,
-                        direction=direction,
-                        base_output_dir=grain_out_path,
-                        image_name=filename,
-                        plot_opts=grain_plot_dict,
-                        **grainstats_config,
-                    )
-                    grainstats_dict[direction], grains_plot_data, height_profiles_dict[direction] = (
-                        grainstats_calculator.calculate_stats()
-                    )
-                    grainstats_dict[direction]["threshold"] = direction
-                    # Plot grains if required
-                    if plotting_config["image_set"] == "all":
-                        LOGGER.info(f"[{filename}] : Plotting grain images for direction: {direction}.")
-                        for plot_data in grains_plot_data:
-                            LOGGER.debug(
-                                f"[{filename}] : Plotting grain image {plot_data['filename']} for direction: {direction}."
-                            )
-                            Images(
-                                data=plot_data["data"],
-                                output_dir=plot_data["output_dir"],
-                                filename=plot_data["filename"],
-                                **plotting_config["plot_dict"][plot_data["name"]],
-                            ).plot_and_save()
+            grain_crops_direction: GrainCropsDirection
+            for direction, grain_crops_direction in image_grain_crops.__dict__.items():
+                if grain_crops_direction is None:
+                    continue
+                grain_crops = grain_crops_direction.crops
+
+                grainstats_calculator = GrainStats(
+                    grain_crops=grain_crops,
+                    pixel_to_nanometre_scaling=pixel_to_nm_scaling,
+                    direction=direction,
+                    base_output_dir=grain_out_path,
+                    image_name=filename,
+                    plot_opts=grain_plot_dict,
+                    **grainstats_config,
+                )
+                grainstats_dict[direction], grains_plot_data, height_profiles_dict[direction] = (
+                    grainstats_calculator.calculate_stats()
+                )
+                grainstats_dict[direction]["threshold"] = direction
+                # Plot grains if required
+                if plotting_config["image_set"] == "all":
+                    LOGGER.info(f"[{filename}] : Plotting grain images for direction: {direction}.")
+                    for plot_data in grains_plot_data:
+                        LOGGER.debug(
+                            f"[{filename}] : Plotting grain image {plot_data['filename']} for direction: {direction}."
+                        )
+                        Images(
+                            data=plot_data["data"],
+                            output_dir=plot_data["output_dir"],
+                            filename=plot_data["filename"],
+                            **plotting_config["plot_dict"][plot_data["name"]],
+                        ).plot_and_save()
             # Create results dataframe from above and below results
             # Appease pylint and ensure that grainstats_df is always created
             grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")

From 433402757b8ffb764feb23913c5d73934f6bae87 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 19 Nov 2024 12:53:33 +0000
Subject: [PATCH 03/34] WIP: Scope out changes to GrainStats.calculate_stats to
 allow for multi class & subgrains

---
 topostats/grainstats.py | 40 +++++++++++++++++++++++++++++++---------
 1 file changed, 31 insertions(+), 9 deletions(-)

diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index c62ef35089..a008235f6a 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -18,6 +18,7 @@
 from topostats.logs.logs import LOGGER_NAME
 from topostats.measure import feret, height_profiles
 from topostats.utils import create_empty_dataframe
+from topostats.grains import GrainCrop
 
 # pylint: disable=too-many-lines
 # pylint: disable=too-many-instance-attributes
@@ -92,8 +93,7 @@ class GrainStats:
 
     def __init__(
         self,
-        data: npt.NDArray,
-        labelled_data: npt.NDArray,
+        grain_crops: dict[int, GrainCrop],
         pixel_to_nanometre_scaling: float,
         direction: str,
         base_output_dir: str | Path,
@@ -134,8 +134,7 @@ def __init__(
             Multiplier to convert the current length scale to metres. Default: 1e-9 for the
             usual AFM length scale of nanometres.
         """
-        self.data = data
-        self.labelled_data = labelled_data
+        self.grain_crops = grain_crops
         self.pixel_to_nanometre_scaling = pixel_to_nanometre_scaling
         self.direction = direction
         self.base_output_dir = Path(base_output_dir)
@@ -203,14 +202,35 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
         """
         grains_plot_data = []
         all_height_profiles = {}
-        if self.labelled_data is None:
+        if len(self.grain_crops) == 0:
             LOGGER.warning(
-                f"[{self.image_name}] : No labelled regions for this image, grain statistics can not be calculated."
+                f"[{self.image_name}] : No grain crops for this image, grain statistics can not be calculated."
             )
             return pd.DataFrame(columns=GRAIN_STATS_COLUMNS), grains_plot_data, all_height_profiles
 
-        # Calculate region properties
-        region_properties = skimage_measure.regionprops(self.labelled_data)
+
+        # Iterate over each grain
+        for grain_crop in self.grain_crops.values():
+            image = grain_crop.image
+            mask = grain_crop.mask
+
+            # Iterate over all the classes except background
+
+                # Split the class into connected components
+
+                # Iterate over all the sub_grains in the class
+
+                    # Skip subgrain if too small to calculate stats for
+
+                    # Create directory for each subgrain's plots
+
+                    # Obtain cropped grain mask and image
+
+                    # Get cropped image and mask
+
+                    # Calculate all the stats
+
+                    # Construct the dictionary of stats for the grain and subgrain
 
         # Iterate over all the grains in the image
         stats_array = []
@@ -598,7 +618,9 @@ def sort_points(self, points: list) -> list:
         # relative to that and _then_ sort it.
         # pivot_angles = self.get_angle(points, self.start_point)
         # Recursively sort the arrays until each point is sorted
-        return self.sort_points(smaller) + sorted(equal, key=self.calculate_squared_distance) + self.sort_points(larger)
+        return (
+            self.sort_points(smaller) + sorted(equal, key=self.calculate_squared_distance) + self.sort_points(larger)
+        )
         # Return sorted array where equal angle points are sorted by distance
 
     def get_start_point(self, edges: npt.NDArray) -> None:

From d106b0fa6f864bae1ea5226efa1b54e77c564b91 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 19 Nov 2024 16:22:25 +0000
Subject: [PATCH 04/34] WIP: Initial proposal for grainstats using grain
 dictionary refactor

---
 topostats/grainstats.py | 250 +++++++++++++++++-----------------------
 1 file changed, 107 insertions(+), 143 deletions(-)

diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index a008235f6a..0ea83bcbd5 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -208,162 +208,126 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
             )
             return pd.DataFrame(columns=GRAIN_STATS_COLUMNS), grains_plot_data, all_height_profiles
 
+        grainstats_dict: dict[int, dict[int, dict[int, dict[int, dict]]]] = {}
 
         # Iterate over each grain
-        for grain_crop in self.grain_crops.values():
+        for grain_index, grain_crop in self.grain_crops.items():
+
+            grainstats_dict[grain_index] = {}
+            all_height_profiles[grain_index] = {}
+
             image = grain_crop.image
             mask = grain_crop.mask
+            grain_bbox = grain_crop.bbox
+            grain_padding = grain_crop.padding
+            grain_anchor = (grain_bbox[0] + grain_padding, grain_bbox[1] + grain_padding)
+
+            # Create directory for grain's plots
+            output_grain = self.base_output_dir / self.direction / f"grain_{grain_index}"
 
             # Iterate over all the classes except background
+            for class_index in range(1, mask.shape[2]):
 
-                # Split the class into connected components
+                grainstats_dict[grain_index][class_index] = {}
+                all_height_profiles[grain_index][class_index] = {}
 
-                # Iterate over all the sub_grains in the class
+                class_mask = mask[:, :, class_index]
 
-                    # Skip subgrain if too small to calculate stats for
+                # Split the class into connected components
+                class_mask_regionprops = skimage_measure.regionprops(skimage_measure.label(class_mask))
 
-                    # Create directory for each subgrain's plots
+                # Iterate over all the sub_grains in the class
+                for subgrain_index, subgrain_region in enumerate(class_mask_regionprops):
 
-                    # Obtain cropped grain mask and image
+                    subgrain_mask = subgrain_region.image
+                    subgrain_mask_image = np.ma.masked_array(
+                        image, mask=np.invert(subgrain_mask), fill_value=np.nan
+                    ).filled()
 
-                    # Get cropped image and mask
+                    # Skip subgrain if too small to calculate stats for
+                    if min(subgrain_mask.shape) < 5:
+                        LOGGER.debug(
+                            f"[{self.image_name}] : Skipping subgrain due to being too small "
+                            f"(size: {subgrain_mask.shape}) to calculate stats for."
+                        )
 
                     # Calculate all the stats
 
-                    # Construct the dictionary of stats for the grain and subgrain
-
-        # Iterate over all the grains in the image
-        stats_array = []
-        # List to hold all the plot data for all the grains. Each entry is a dictionary of plotting data.
-        # There are multiple entries for each grain.
-        for index, region in enumerate(region_properties):
-            LOGGER.debug(f"[{self.image_name}] : Processing grain: {index}")
-
-            # Skip grain if too small to calculate stats for
-            LOGGER.debug(f"[{self.image_name}] : Grain size: {region.image.size}")
-            if min(region.image.shape) < 5:
-                LOGGER.debug(
-                    f"[{self.image_name}] : Skipping grain due to being too small (size: {region.image.shape}) to calculate stats for."
-                )
-                continue
-
-            # Create directory for each grain's plots
-            output_grain = self.base_output_dir / self.direction
-            # Obtain cropped grain mask and image
-            minr, minc, maxr, maxc = region.bbox
-            grain_mask = np.array(region.image)
-            grain_image = self.data[minr:maxr, minc:maxc]
-            grain_mask_image = np.ma.masked_array(grain_image, mask=np.invert(grain_mask), fill_value=np.nan).filled()
-
-            if self.cropped_size == -1:
-                for name, image in {
-                    "grain_image": grain_image,
-                    "grain_mask": grain_mask,
-                    "grain_mask_image": grain_mask_image,
-                }.items():
-                    grains_plot_data.append(
-                        {
-                            "data": image,
-                            "output_dir": output_grain,
-                            "filename": f"{self.image_name}_{name}_{index}",
-                            "name": name,
-                        }
+                    points = self.calculate_points(subgrain_mask)
+                    edges = self.calculate_edges(subgrain_mask, edge_detection_method=self.edge_detection_method)
+                    radius_stats = self.calculate_radius_stats(edges, points)
+                    # hull, hull_indices, hull_simplexes = self.convex_hull(edges, output_grain)
+                    _, _, hull_simplexes = self.convex_hull(edges, output_grain)
+                    centroid = self._calculate_centroid(points)
+                    # Centroids for the grains (grain anchor added because centroid returns values local to the
+                    # cropped grain images)
+                    centre_x = centroid[0] + grain_anchor[0]
+                    centre_y = centroid[1] + grain_anchor[1]
+                    (
+                        smallest_bounding_width,
+                        smallest_bounding_length,
+                        aspect_ratio,
+                    ) = self.calculate_aspect_ratio(
+                        edges=edges,
+                        hull_simplices=hull_simplexes,
+                        path=output_grain,
                     )
 
-            else:
-                # Get cropped image and mask
-                grain_centre = int((minr + maxr) / 2), int((minc + maxc) / 2)
-                length = int(self.cropped_size / (2 * self.pixel_to_nanometre_scaling))
-                solo_mask = self.labelled_data.copy()
-                solo_mask[solo_mask != index + 1] = 0
-                solo_mask[solo_mask == index + 1] = 1
-                cropped_grain_image = self.get_cropped_region(self.data, length, np.asarray(grain_centre))
-                cropped_grain_mask = self.get_cropped_region(solo_mask, length, np.asarray(grain_centre)).astype(bool)
-                cropped_grain_mask_image = np.ma.masked_array(
-                    grain_image, mask=np.invert(grain_mask), fill_value=np.nan
-                ).filled()
-                for name, image in {
-                    "grain_image": cropped_grain_image,
-                    "grain_mask": cropped_grain_mask,
-                    "grain_mask_image": cropped_grain_mask_image,
-                }.items():
-                    grains_plot_data.append(
-                        {
-                            "data": image,
-                            "output_dir": output_grain,
-                            "filename": f"{self.image_name}_{name}_{index}",
-                            "name": name,
-                        }
-                    )
-
-            points = self.calculate_points(grain_mask)
-            edges = self.calculate_edges(grain_mask, edge_detection_method=self.edge_detection_method)
-            radius_stats = self.calculate_radius_stats(edges, points)
-            # hull, hull_indices, hull_simplexes = self.convex_hull(edges, output_grain)
-            _, _, hull_simplexes = self.convex_hull(edges, output_grain)
-            centroid = self._calculate_centroid(points)
-            # Centroids for the grains (minc and minr added because centroid returns values local to the cropped grain images)
-            centre_x = centroid[0] + minc
-            centre_y = centroid[1] + minr
-            (
-                smallest_bounding_width,
-                smallest_bounding_length,
-                aspect_ratio,
-            ) = self.calculate_aspect_ratio(
-                edges=edges,
-                hull_simplices=hull_simplexes,
-                path=output_grain,
-            )
-
-            # Calculate scaling factors
-            length_scaling_factor = self.pixel_to_nanometre_scaling * self.metre_scaling_factor
-            area_scaling_factor = length_scaling_factor**2
-
-            # Calculate minimum and maximum feret diameters and scale the distances
-            feret_statistics = feret.min_max_feret(points)
-            feret_statistics["min_feret"] = feret_statistics["min_feret"] * length_scaling_factor
-            feret_statistics["max_feret"] = feret_statistics["max_feret"] * length_scaling_factor
-
-            if self.extract_height_profile:
-                all_height_profiles[index] = height_profiles.interpolate_height_profile(
-                    img=grain_image, mask=grain_mask
-                )
-                LOGGER.debug(f"[{self.image_name}] : Height profiles extracted.")
-
-            # Save the stats to dictionary. Note that many of the stats are multiplied by a scaling factor to convert
-            # from pixel units to nanometres.
-            # Removed formatting, better to keep accurate until the end, including in CSV, then shorten display
-            stats = {
-                "centre_x": centre_x * length_scaling_factor,
-                "centre_y": centre_y * length_scaling_factor,
-                "radius_min": radius_stats["min"] * length_scaling_factor,
-                "radius_max": radius_stats["max"] * length_scaling_factor,
-                "radius_mean": radius_stats["mean"] * length_scaling_factor,
-                "radius_median": radius_stats["median"] * length_scaling_factor,
-                "height_min": np.nanmin(grain_mask_image) * self.metre_scaling_factor,
-                "height_max": np.nanmax(grain_mask_image) * self.metre_scaling_factor,
-                "height_median": np.nanmedian(grain_mask_image) * self.metre_scaling_factor,
-                "height_mean": np.nanmean(grain_mask_image) * self.metre_scaling_factor,
-                # [volume] = [pixel] * [pixel] * [height] = px * px * nm.
-                # To turn into m^3, multiply by pixel_to_nanometre_scaling^2 and metre_scaling_factor^3.
-                "volume": np.nansum(grain_mask_image)
-                * self.pixel_to_nanometre_scaling**2
-                * (self.metre_scaling_factor**3),
-                "area": region.area * area_scaling_factor,
-                "area_cartesian_bbox": region.area_bbox * area_scaling_factor,
-                "smallest_bounding_width": smallest_bounding_width * length_scaling_factor,
-                "smallest_bounding_length": smallest_bounding_length * length_scaling_factor,
-                "smallest_bounding_area": smallest_bounding_length * smallest_bounding_width * area_scaling_factor,
-                "aspect_ratio": aspect_ratio,
-                "threshold": self.direction,
-                "max_feret": feret_statistics["max_feret"],
-                "min_feret": feret_statistics["min_feret"],
-            }
-            stats_array.append(stats)
-        if len(stats_array) > 0:
-            grainstats_df = pd.DataFrame(data=stats_array)
-        else:
+                    # Calculate scaling factors
+                    length_scaling_factor = self.pixel_to_nanometre_scaling * self.metre_scaling_factor
+                    area_scaling_factor = length_scaling_factor**2
+
+                    # Calculate minimum and maximum feret diameters and scale the distances
+                    feret_statistics = feret.min_max_feret(points)
+                    feret_statistics["min_feret"] = feret_statistics["min_feret"] * length_scaling_factor
+                    feret_statistics["max_feret"] = feret_statistics["max_feret"] * length_scaling_factor
+
+                    if self.extract_height_profile:
+                        all_height_profiles[grain_index][class_index][subgrain_index] = (
+                            height_profiles.interpolate_height_profile(img=image, mask=subgrain_mask)
+                        )
+                        LOGGER.debug(f"[{self.image_name}] : Height profiles extracted.")
+
+                    # Save the stats to dictionary. Note that many of the stats are multiplied by a scaling factor to convert
+                    # from pixel units to nanometres.
+                    # Removed formatting, better to keep accurate until the end, including in CSV, then shorten display
+                    stats = {
+                        "centre_x": centre_x * length_scaling_factor,
+                        "centre_y": centre_y * length_scaling_factor,
+                        "radius_min": radius_stats["min"] * length_scaling_factor,
+                        "radius_max": radius_stats["max"] * length_scaling_factor,
+                        "radius_mean": radius_stats["mean"] * length_scaling_factor,
+                        "radius_median": radius_stats["median"] * length_scaling_factor,
+                        "height_min": np.nanmin(subgrain_mask_image) * self.metre_scaling_factor,
+                        "height_max": np.nanmax(subgrain_mask_image) * self.metre_scaling_factor,
+                        "height_median": np.nanmedian(subgrain_mask_image) * self.metre_scaling_factor,
+                        "height_mean": np.nanmean(subgrain_mask_image) * self.metre_scaling_factor,
+                        # [volume] = [pixel] * [pixel] * [height] = px * px * nm.
+                        # To turn into m^3, multiply by pixel_to_nanometre_scaling^2 and metre_scaling_factor^3.
+                        "volume": np.nansum(subgrain_mask_image)
+                        * self.pixel_to_nanometre_scaling**2
+                        * (self.metre_scaling_factor**3),
+                        "area": subgrain_region.area * area_scaling_factor,
+                        "area_cartesian_bbox": subgrain_region.area_bbox * area_scaling_factor,
+                        "smallest_bounding_width": smallest_bounding_width * length_scaling_factor,
+                        "smallest_bounding_length": smallest_bounding_length * length_scaling_factor,
+                        "smallest_bounding_area": smallest_bounding_length
+                        * smallest_bounding_width
+                        * area_scaling_factor,
+                        "aspect_ratio": aspect_ratio,
+                        "threshold": self.direction,
+                        "max_feret": feret_statistics["max_feret"],
+                        "min_feret": feret_statistics["min_feret"],
+                    }
+
+                    grainstats_dict[grain_index][class_index][subgrain_index] = stats
+
+        grainstats_df = pd.DataFrame(data=grainstats_dict)
+
+        # Check if the dataframe is empty
+        if grainstats_df.empty:
             grainstats_df = create_empty_dataframe()
+
         grainstats_df.index.name = "grain_number"
         grainstats_df["image"] = self.image_name
 
@@ -434,7 +398,7 @@ def calculate_edges(grain_mask: npt.NDArray, edge_detection_method: str) -> list
         # return np.stack(nonzero_coordinates, axis=1)
         return [list(vector) for vector in np.transpose(nonzero_coordinates)]
 
-    def calculate_radius_stats(self, edges: list, points: list) -> tuple[float]:
+    def calculate_radius_stats(self, edges: list, points: list) -> dict[str, float]:
         """
         Calculate the radius of grains.
 
@@ -449,8 +413,8 @@ def calculate_radius_stats(self, edges: list, points: list) -> tuple[float]:
 
         Returns
         -------
-        tuple[float]
-            A tuple of the minimum, maximum, mean and median radius of the grain.
+        dict[str, float]
+            A dictionary containing the minimum, maximum, mean and median radius of the grain.
         """
         # Calculate the centroid of the grain
         centroid = self._calculate_centroid(points)
@@ -461,8 +425,8 @@ def calculate_radius_stats(self, edges: list, points: list) -> tuple[float]:
         return {
             "min": np.min(radii),
             "max": np.max(radii),
-            "mean": np.mean(radii),
-            "median": np.median(radii),
+            "mean": float(np.mean(radii)),
+            "median": float(np.median(radii)),
         }
 
     @staticmethod

From 1acee6792b0380c53cdbcd65fa8dbb6b4a0cd23b Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 09:41:51 +0000
Subject: [PATCH 05/34] Add function: graincrops_merge_classes

---
 topostats/grains.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/topostats/grains.py b/topostats/grains.py
index ba382e21cb..e8cbb71957 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -1723,3 +1723,34 @@ def extract_grains_from_full_image_mask(
             )
 
         return graincrops
+
+    @staticmethod
+    def graincrops_merge_classes(
+        graincrops: dict[int, GrainCrop],
+        classes_to_merge: list[tuple[int]] | None,
+    ) -> dict[int, GrainCrop]:
+        """
+        Merge classes in the grain crops.
+
+        Parameters
+        ----------
+        graincrops : dict[int, GrainCrop]
+            Dictionary of grain crops.
+        classes_to_merge : list | None
+            List of tuples for classes to merge, can be any number of classes.
+
+        Returns
+        -------
+        dict[int, GrainCrop]
+            Dictionary of grain crops with classes merged.
+        """
+        if classes_to_merge is None:
+            return graincrops
+
+        for grain_number, graincrop in graincrops.items():
+            graincrop.mask = Grains.merge_classes(
+                grain_mask_tensor=graincrop.mask,
+                classes_to_merge=classes_to_merge,
+            )
+
+        return graincrops

From a7b90754b72e177e51a5860b27b47a79468fbd65 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 09:42:20 +0000
Subject: [PATCH 06/34] Add function: graincrops_update_background_class

---
 topostats/grains.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/topostats/grains.py b/topostats/grains.py
index e8cbb71957..f4fb779b29 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -1754,3 +1754,25 @@ def graincrops_merge_classes(
             )
 
         return graincrops
+
+    @staticmethod
+    def graincrops_update_background_class(
+        graincrops: dict[int, GrainCrop],
+    ) -> dict[int, GrainCrop]:
+        """
+        Update the background class in the grain crops.
+
+        Parameters
+        ----------
+        graincrops : dict[int, GrainCrop]
+            Dictionary of grain crops.
+
+        Returns
+        -------
+        dict[int, GrainCrop]
+            Dictionary of grain crops with updated background class.
+        """
+        for grain_number, graincrop in graincrops.items():
+            graincrop.mask = Grains.update_background_class(graincrop.mask)
+
+        return graincrops

From 913e3b55db4237ef5837963dd4da09ae0152dec0 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 10:18:05 +0000
Subject: [PATCH 07/34] Update: extract_grains_from_full_image now works in
 theory, untested

---
 topostats/grains.py | 77 +++++++++++++++++++++++++++++++--------------
 1 file changed, 53 insertions(+), 24 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index f4fb779b29..09ee61b63b 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -1664,12 +1664,14 @@ def construct_full_mask_from_crops(graincrops: dict[int, GrainCrop], image_shape
     @staticmethod
     def extract_grains_from_full_image_mask(
         image: npt.NDArray[np.float32],
-        labelled_full_mask_tensor: npt.NDArray[np.bool_],
+        full_mask_tensor: npt.NDArray[np.bool_],
         padding: int,
     ) -> dict[int, GrainCrop]:
         """
         Extract grains from the full image mask tensor.
 
+        Grains are detected using connected components across all classes in the full mask tensor.
+
         Parameters
         ----------
         image : npt.NDArray[np.float32]
@@ -1684,42 +1686,69 @@ def extract_grains_from_full_image_mask(
         dict[int, GrainCrop]
             Dictionary of grain crops.
         """
-        grain_region_properties = regionprops(labelled_full_mask_tensor)
+        # Flatten the mask tensor
+        flat_mask = Grains.flatten_multi_class_tensor(full_mask_tensor)
+        labelled_flat_full_mask = label(flat_mask)
+        flat_regionprops_full_mask = regionprops(labelled_flat_full_mask)
         graincrops = {}
-        for grain_number, region in enumerate(grain_region_properties):
+        for grain_number, flat_region in enumerate(flat_regionprops_full_mask):
+            # Get a flattened binary mask for the whole grain and no other grains
+            flattened_grain_binary_mask = labelled_flat_full_mask == flat_region.label
+
+            # For each class, set all pixels to zero that are not in the current region
+            for class_index in range(1, full_mask_tensor.shape[2]):
+                # Set all pixels to zero that are not in the current region's pixels by multiplying by a binary mask
+                # for the whole flattened grain mask
+                grain_tensor_full_mask = np.zeros_like(full_mask_tensor)
+                grain_tensor_full_mask[:, :, class_index] = (
+                    flattened_grain_binary_mask * full_mask_tensor[:, :, class_index]
+                )
+
+            # Crop the tensor
             # Get the bounding box for the region
-            bounding_box: tuple[int, int, int, int] = tuple(region.bbox)  # min_row, min_col, max_row, max_col
-            mask_crop = labelled_full_mask_tensor[
-                bounding_box[0] : bounding_box[2],
-                bounding_box[1] : bounding_box[3],
-            ]
+            flat_bounding_box: tuple[int, int, int, int] = tuple(
+                flat_region.bbox
+            )  # min_row, min_col, max_row, max_col
 
             # Pad the mask
-            padded_mask, bounding_box = pad_crop(
-                crop=mask_crop,
-                bbox=bounding_box,
-                image_shape=(image.shape[0], image.shape[1]),
+            padded_flat_bounding_box = pad_bounding_box(
+                crop_min_row=flat_bounding_box[0],
+                crop_min_col=flat_bounding_box[1],
+                crop_max_row=flat_bounding_box[2],
+                crop_max_col=flat_bounding_box[3],
+                image_shape=(full_mask_tensor.shape[0], full_mask_tensor.shape[1]),
                 padding=padding,
             )
 
-            square_padded_mask, bounding_box = make_crop_square(
-                crop=padded_mask,
-                bbox=bounding_box,
-                image_shape=(image.shape[0], image.shape[1]),
+            # Make the mask square
+            square_flat_bounding_box = make_bounding_box_square(
+                crop_min_row=padded_flat_bounding_box[0],
+                crop_min_col=padded_flat_bounding_box[1],
+                crop_max_row=padded_flat_bounding_box[2],
+                crop_max_col=padded_flat_bounding_box[3],
+                image_shape=(full_mask_tensor.shape[0], full_mask_tensor.shape[1]),
             )
 
-            # Grab the cropped image. Using slice since the bounding box from skimage is
-            # half-open, so the max_row and max_col are not included in the region.
-            region_image = image[
-                bounding_box[0] : bounding_box[2],
-                bounding_box[1] : bounding_box[3],
+            # Grab image and mask for the cropped region
+            grain_cropped_image = image[
+                square_flat_bounding_box[0] : square_flat_bounding_box[2],
+                square_flat_bounding_box[1] : square_flat_bounding_box[3],
             ]
 
+            grain_cropped_tensor = grain_tensor_full_mask[
+                square_flat_bounding_box[0] : square_flat_bounding_box[2],
+                square_flat_bounding_box[1] : square_flat_bounding_box[3],
+                :,
+            ]
+
+            # Update background class to reflect the removal of any non-connected grains
+            grain_cropped_tensor = Grains.update_background_class(grain_mask_tensor=grain_cropped_tensor)
+
             graincrops[grain_number] = GrainCrop(
-                image=region_image,
-                mask=square_padded_mask,
+                image=grain_cropped_image,
+                mask=grain_cropped_tensor,
                 padding=padding,
-                bbox=bounding_box,
+                bbox=flat_bounding_box,
             )
 
         return graincrops

From 843a606cc098ed847e21015ca92a4b2e38a54228 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 10:37:32 +0000
Subject: [PATCH 08/34] Fix: extract_grains_from_full_image_mask: allocating
 region to empty mask required bool. Tested in debugger.

---
 topostats/grains.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 09ee61b63b..b6e35c1b11 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -1696,13 +1696,13 @@ def extract_grains_from_full_image_mask(
             flattened_grain_binary_mask = labelled_flat_full_mask == flat_region.label
 
             # For each class, set all pixels to zero that are not in the current region
+            grain_tensor_full_mask = np.zeros_like(full_mask_tensor).astype(bool)
             for class_index in range(1, full_mask_tensor.shape[2]):
                 # Set all pixels to zero that are not in the current region's pixels by multiplying by a binary mask
                 # for the whole flattened grain mask
-                grain_tensor_full_mask = np.zeros_like(full_mask_tensor)
                 grain_tensor_full_mask[:, :, class_index] = (
                     flattened_grain_binary_mask * full_mask_tensor[:, :, class_index]
-                )
+                ).astype(bool)
 
             # Crop the tensor
             # Get the bounding box for the region

From 535bea8607ba2c344f42f297c6038aeeda4edf02 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 10:40:09 +0000
Subject: [PATCH 09/34] WIP: Switch vetting, merging and update background to
 work on dicts of GrainCrops

---
 topostats/grains.py | 40 +++++++++++++++-------------------------
 1 file changed, 15 insertions(+), 25 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index b6e35c1b11..0e84f1c3fb 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -95,8 +95,8 @@ class GrainCropsDirection:
         Grain crops.
     """
 
-    full_mask_tensor: npt.NDArray[np.bool_]
     crops: dict[int, GrainCrop]
+    full_mask_tensor: npt.NDArray[np.bool_]
 
     def __post_init__(self):
         """
@@ -616,7 +616,6 @@ def find_grains(self):
             )
             self.bounding_boxes[direction] = self.get_bounding_boxes(direction=direction)
             LOGGER.debug(f"[{self.filename}] : Extracted bounding boxes ({direction})")
-            thresholding_grain_count = self.directions[direction]["labelled_regions_02"].max()
 
             # Force labelled_regions_02 to be of shape NxNx2, where the two classes are a binary background mask and the second is a binary grain mask.
             # This is because we want to support multiple classes, and so we standardise so that the first layer is background mask, then feature mask 1, then feature mask 2 etc.
@@ -655,7 +654,7 @@ def find_grains(self):
 
             graincrops = self.extract_grains_from_full_image_mask(
                 image=self.image,
-                labelled_full_mask_tensor=self.directions[direction]["labelled_regions_02"],
+                full_mask_tensor=self.directions[direction]["labelled_regions_02"],
                 padding=self.grain_crop_padding,
             )
 
@@ -680,37 +679,28 @@ def find_grains(self):
 
             # Vet the grains
             if self.vetting is not None:
-                vetted_grains = Grains.vet_grains(
-                    grain_mask_tensor=self.directions[direction]["labelled_regions_02"].astype(bool),
+                graincrops_vetted = Grains.vet_grains(
+                    graincrops=graincrops,
                     pixel_to_nm_scaling=self.pixel_to_nm_scaling,
                     **self.vetting,
                 )
             else:
-                vetted_grains = self.directions[direction]["labelled_regions_02"].astype(bool)
+                graincrops_vetted = graincrops
 
             # Merge classes if necessary
-            merged_classes = Grains.merge_classes(
-                vetted_grains,
-                self.classes_to_merge,
+            graincrops_merged_classes = Grains.graincrops_merge_classes(
+                graincrops=graincrops_vetted,
+                classes_to_merge=self.classes_to_merge,
             )
 
             # Update the background class
-            final_grains = Grains.update_background_class(grain_mask_tensor=merged_classes)
-
-            # Label each class in the tensor
-            labelled_final_grains = np.zeros_like(final_grains).astype(int)
-            # The background class will be the same as the binary mask
-            labelled_final_grains[:, :, 0] = final_grains[:, :, 0]
-            # Iterate over each class and label the regions
-            for class_index in range(final_grains.shape[2]):
-                labelled_final_grains[:, :, class_index] = Grains.label_regions(final_grains[:, :, class_index])
-
-            self.directions[direction]["removed_small_objects"] = labelled_final_grains.astype(bool)
-            self.directions[direction]["labelled_regions_02"] = labelled_final_grains.astype(np.int32)
+            graincrops_updated_background = Grains.graincrops_update_background_class(
+                graincrops=graincrops_merged_classes
+            )
 
             self.directions[direction]["grain_crops_direction"] = GrainCropsDirection(
                 full_mask_tensor=full_mask_tensor,
-                crops=graincrops,
+                crops=graincrops_updated_background,
             )
 
         if "above" in self.directions and "below" in self.directions:
@@ -820,8 +810,6 @@ def improve_grain_segmentation_unet(
                 bbox=graincrop.bbox,
             )
 
-                # Grab the unet mask for the class
-                unet_predicted_mask_labelled = morphology.label(predicted_mask[:, :, class_index])
         return graincrops
 
     @staticmethod
@@ -969,7 +957,9 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
+                vetting_criteria[1:]
+                for vetting_criteria in class_size_thresholds
+                if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:

From 36c07e25ef8f82ff5c80da183671cebebfd93476 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 10:40:45 +0000
Subject: [PATCH 10/34] WIP: Update vet_grains to take / return dicts of
 GrainCrops

---
 topostats/grains.py | 44 +++++++++++++++++---------------------------
 1 file changed, 17 insertions(+), 27 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 0e84f1c3fb..008f3c1fb2 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -1472,7 +1472,7 @@ def convert_classes_when_too_big_or_small(
 
     @staticmethod
     def vet_grains(
-        grain_mask_tensor: npt.NDArray,
+        graincrops: dict[int, GrainCrop],
         pixel_to_nm_scaling: float,
         class_conversion_size_thresholds: list[tuple[tuple[int, int, int], tuple[int, int]]] | None,
         class_size_thresholds: list[tuple[int, int, int]] | None,
@@ -1481,14 +1481,14 @@ def vet_grains(
         class_touching_threshold: int,
         keep_largest_labelled_regions_classes: list[int] | None,
         class_connection_point_thresholds: list[tuple[tuple[int, int], tuple[int, int]]] | None,
-    ) -> npt.NDArray:
+    ) -> dict[int, GrainCrop]:
         """
         Vet grains in a grain mask tensor based on a variety of criteria.
 
         Parameters
         ----------
-        grain_mask_tensor : npt.NDArray
-            3-D Numpy array of the grain mask tensor.
+        graincrops : dict[int, GrainCrop]
+            Dictionary of grain crops.
         pixel_to_nm_scaling : float
             Scaling of pixels to nanometres.
         class_conversion_size_thresholds : list
@@ -1509,16 +1509,15 @@ def vet_grains(
 
         Returns
         -------
-        npt.NDArray
-            3-D Numpy array of the vetted grain mask tensor.
+        dict[int, GrainCrop]
+            Dictionary of grain crops that passed the vetting.
         """
-        # Get individual grain crops
-        grain_tensor_crops, bounding_boxes, padding = Grains.get_individual_grain_crops(grain_mask_tensor)
-
-        passed_grain_crops_and_bounding_boxes = []
+        passed_graincrops: dict[int, GrainCrop] = {}
 
         # Iterate over the grain crops
-        for _, (single_grain_mask_tensor, bounding_box) in enumerate(zip(grain_tensor_crops, bounding_boxes)):
+        for grain_number, graincrop in graincrops.items():
+
+            single_grain_mask_tensor = graincrop.mask
 
             # Convert small / big areas to other classes
             single_grain_mask_tensor = Grains.convert_classes_when_too_big_or_small(
@@ -1565,24 +1564,15 @@ def vet_grains(
             ):
                 continue
 
-            # If passed all vetting steps, add to the list of passed grain crops
-            passed_grain_crops_and_bounding_boxes.append(
-                {
-                    "grain_tensor": largest_only_single_grain_mask_tensor,
-                    "bounding_box": bounding_box,
-                    "padding": padding,
-                }
+            # If passed all vetting steps, add to the dictionary of passed grain crops
+            passed_graincrops[grain_number] = GrainCrop(
+                image=graincrop.image,
+                mask=largest_only_single_grain_mask_tensor,
+                padding=graincrop.padding,
+                bbox=graincrop.bbox,
             )
 
-        # Construct a new grain mask tensor from the passed grains
-        return Grains.assemble_grain_mask_tensor_from_crops(
-            grain_mask_tensor_shape=(
-                grain_mask_tensor.shape[0],
-                grain_mask_tensor.shape[1],
-                grain_mask_tensor.shape[2],
-            ),
-            grain_crops_and_bounding_boxes=passed_grain_crops_and_bounding_boxes,
-        )
+        return passed_graincrops
 
     @staticmethod
     def merge_classes(

From d1e4c1d5838898040b8d694acb42c8d8b198825e Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Thu, 21 Nov 2024 10:41:34 +0000
Subject: [PATCH 11/34] WIP: Update grainstats handling of dataframe to use
 list of dicts for each row

---
 topostats/grainstats.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index 0ea83bcbd5..56b4b49dce 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -208,12 +208,10 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
             )
             return pd.DataFrame(columns=GRAIN_STATS_COLUMNS), grains_plot_data, all_height_profiles
 
-        grainstats_dict: dict[int, dict[int, dict[int, dict[int, dict]]]] = {}
+        grainstats_rows: list[dict] = []
 
         # Iterate over each grain
         for grain_index, grain_crop in self.grain_crops.items():
-
-            grainstats_dict[grain_index] = {}
             all_height_profiles[grain_index] = {}
 
             image = grain_crop.image
@@ -227,8 +225,6 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
 
             # Iterate over all the classes except background
             for class_index in range(1, mask.shape[2]):
-
-                grainstats_dict[grain_index][class_index] = {}
                 all_height_profiles[grain_index][class_index] = {}
 
                 class_mask = mask[:, :, class_index]
@@ -292,6 +288,9 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
                     # from pixel units to nanometres.
                     # Removed formatting, better to keep accurate until the end, including in CSV, then shorten display
                     stats = {
+                        "grain_number": grain_index,
+                        "class_number": class_index,
+                        "subgrain_number": subgrain_index,
                         "centre_x": centre_x * length_scaling_factor,
                         "centre_y": centre_y * length_scaling_factor,
                         "radius_min": radius_stats["min"] * length_scaling_factor,
@@ -320,9 +319,7 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
                         "min_feret": feret_statistics["min_feret"],
                     }
 
-                    grainstats_dict[grain_index][class_index][subgrain_index] = stats
-
-        grainstats_df = pd.DataFrame(data=grainstats_dict)
+                    grainstats_rows.append(stats)
 
         # Check if the dataframe is empty
         if grainstats_df.empty:

From ace46ea912fb4011e550d87d1c579903b0570f37 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 10:00:07 +0000
Subject: [PATCH 12/34] Fix: validate_full_mask_tensor_shape: Require
 len(shape) == 3, shape[2] >= 2, shape[1]==shape[2]

---
 topostats/grains.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 008f3c1fb2..2459da98d2 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -77,8 +77,8 @@ def validate_full_mask_tensor_shape(array: npt.NDArray[np.bool_]) -> npt.NDArray
     npt.NDArray
         Numpy array if valid.
     """
-    if len(array.shape) != 3 or array.shape[2] != 3 or array.shape[1] != array.shape[0]:
-        raise ValueError(f"Full mask tensor must be NxNx3 but has shape {array.shape}")
+    if len(array.shape) != 3 or array.shape[2] < 2 or array.shape[1] != array.shape[0]:
+        raise ValueError(f"Full mask tensor must be NxNxC with C >= 2 but has shape {array.shape}")
     return array
 
 

From e2824f9643d680dd8463ee3c1394d518fd24279c Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 10:02:16 +0000
Subject: [PATCH 13/34] Edit: find_grains now stores grains in
 self.image_grain_crops: ImageGrainCrops. Locally debugged working

---
 topostats/grains.py | 96 ++++++++++++++++++---------------------------
 1 file changed, 38 insertions(+), 58 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 2459da98d2..8a47763abe 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -287,6 +287,11 @@ def __init__(
         self.minimum_grain_size_px = 10
         self.minimum_bbox_size_px = 5
 
+        self.image_grain_crops = ImageGrainCrops(
+            above=None,
+            below=None,
+        )
+
     def tidy_border(self, image: npt.NDArray, **kwargs) -> npt.NDArray:
         """
         Remove grains touching the border.
@@ -547,7 +552,7 @@ def get_bounding_boxes(self, direction: str) -> dict:
         """
         return {region.area: region.area_bbox for region in self.region_properties[direction]}
 
-    def find_grains(self):
+    def find_grains(self) -> None:
         """Find grains."""
         LOGGER.debug(f"[{self.filename}] : Thresholding method (grains) : {self.threshold_method}")
         self.thresholds = get_thresholds(
@@ -558,6 +563,9 @@ def find_grains(self):
             absolute=self.threshold_absolute,
         )
 
+        # Create an ImageGrainCrops object to store the grain crops
+        image_grain_crops = ImageGrainCrops(above=None, below=None)
+
         for direction in self.direction:
             LOGGER.debug(f"[{self.filename}] : Finding {direction} grains, threshold: ({self.thresholds[direction]})")
             self.directions[direction] = {}
@@ -617,50 +625,32 @@ def find_grains(self):
             self.bounding_boxes[direction] = self.get_bounding_boxes(direction=direction)
             LOGGER.debug(f"[{self.filename}] : Extracted bounding boxes ({direction})")
 
-            # Force labelled_regions_02 to be of shape NxNx2, where the two classes are a binary background mask and the second is a binary grain mask.
-            # This is because we want to support multiple classes, and so we standardise so that the first layer is background mask, then feature mask 1, then feature mask 2 etc.
+            # Create a tensor out of the grain mask of shape NxNx2, where the two classes are a binary background
+            # mask and the second is a binary grain mask. This is because we want to support multiple classes, and
+            # so we standardise so that the first layer is background mask, then feature mask 1, then feature mask
+            # 2 etc.
 
             # Get a binary mask where 1s are background and 0s are grains
             labelled_regions_background_mask = np.where(self.directions[direction]["labelled_regions_02"] == 0, 1, 0)
-            # keep only the largest region
-            labelled_regions_background_mask = label(labelled_regions_background_mask)
-            areas = [region.area for region in regionprops(labelled_regions_background_mask)]
-            labelled_regions_background_mask = np.where(
-                labelled_regions_background_mask == np.argmax(areas) + 1, labelled_regions_background_mask, 0
-            )
 
-            self.directions[direction]["labelled_regions_02"] = np.stack(
+            # Create a tensor out of the background and foreground masks
+            full_mask_tensor = np.stack(
                 [
                     labelled_regions_background_mask,
                     self.directions[direction]["labelled_regions_02"],
                 ],
                 axis=-1,
-            ).astype(
-                np.int32
-            )  # Will produce an NxNx2 array
-
-            # Do the same for removed_small_objects, using the same labelled_regions_backgroudn_mask as the background since they will be the same
-            self.directions[direction]["removed_small_objects"] = np.stack(
-                [
-                    labelled_regions_background_mask,
-                    self.directions[direction]["removed_small_objects"],
-                ],
-                axis=-1,
-            ).astype(
-                np.int32
-            )  # Will produce an NxNx2 array
+            ).astype(np.int32)
 
-            full_mask_tensor = self.directions[direction]["labelled_regions_02"]
-
-            graincrops = self.extract_grains_from_full_image_mask(
+            # Extract tensor mask crops of each grain.
+            graincrops = self.extract_grains_from_full_image_tensor(
                 image=self.image,
-                full_mask_tensor=self.directions[direction]["labelled_regions_02"],
+                full_mask_tensor=full_mask_tensor,
                 padding=self.grain_crop_padding,
             )
 
-            # Check whether to run the UNet model
+            # Optionally run a user-supplied u-net model on the grains to improve the segmentation
             if self.unet_config["model_path"] is not None:
-
                 # Run unet segmentation on only the class 1 layer of the labelled_regions_02. Need to make this configurable
                 # later on along with all the other hardcoded class 1s.
                 graincrops = Grains.improve_grain_segmentation_unet(
@@ -670,8 +660,7 @@ def find_grains(self):
                     image=self.image,
                     graincrops=graincrops,
                 )
-
-                # Construct full mask from the crops
+                # Construct full masks from the crops
                 full_mask_tensor = Grains.construct_full_mask_from_crops(
                     graincrops=graincrops,
                     image_shape=self.image.shape,
@@ -687,38 +676,32 @@ def find_grains(self):
             else:
                 graincrops_vetted = graincrops
 
-            # Merge classes if necessary
+            # Merge classes as specified by the user
             graincrops_merged_classes = Grains.graincrops_merge_classes(
                 graincrops=graincrops_vetted,
                 classes_to_merge=self.classes_to_merge,
             )
 
-            # Update the background class
+            # Update the background class to ensure the background is accurate
             graincrops_updated_background = Grains.graincrops_update_background_class(
                 graincrops=graincrops_merged_classes
             )
 
-            self.directions[direction]["grain_crops_direction"] = GrainCropsDirection(
-                full_mask_tensor=full_mask_tensor,
-                crops=graincrops_updated_background,
-            )
+            # Store the grain crops
+            if direction == "above":
+                image_grain_crops.above = GrainCropsDirection(
+                    crops=graincrops_updated_background,
+                    full_mask_tensor=full_mask_tensor,
+                )
+            elif direction == "below":
+                image_grain_crops.below = GrainCropsDirection(
+                    crops=graincrops_updated_background,
+                    full_mask_tensor=full_mask_tensor,
+                )
+            else:
+                raise ValueError(f"Invalid direction: {direction}. Allowed values are 'above' and 'below'")
 
-        if "above" in self.directions and "below" in self.directions:
-            return ImageGrainCrops(
-                above=self.directions["above"]["grain_crops_direction"],
-                below=self.directions["below"]["grain_crops_direction"],
-            )
-        if "above" in self.directions:
-            return ImageGrainCrops(
-                above=self.directions["above"]["grain_crops_direction"],
-                below=None,
-            )
-        if "below" in self.directions:
-            return ImageGrainCrops(
-                above=None,
-                below=self.directions["below"]["grain_crops_direction"],
-            )
-        raise ValueError("No grains found in either direction")
+        self.image_grain_crops = image_grain_crops
 
     # pylint: disable=too-many-locals
     @staticmethod
@@ -726,7 +709,6 @@ def improve_grain_segmentation_unet(
         filename: str,
         direction: str,
         unet_config: dict[str, str | int | float | tuple[int | None, int, int, int] | None],
-        image: npt.NDArray,
         graincrops: dict[int, GrainCrop],
     ) -> dict[int, GrainCrop]:
         """
@@ -748,8 +730,6 @@ def improve_grain_segmentation_unet(
                 Upper bound for normalising the image.
             lower_norm_bound: float
                 Lower bound for normalising the image.
-        image : npt.NDArray
-            2-D Numpy array of image.
         labelled_grain_regions : npt.NDArray
             2-D Numpy array of labelled grain regions.
 
@@ -1642,7 +1622,7 @@ def construct_full_mask_from_crops(graincrops: dict[int, GrainCrop], image_shape
         return full_mask_tensor
 
     @staticmethod
-    def extract_grains_from_full_image_mask(
+    def extract_grains_from_full_image_tensor(
         image: npt.NDArray[np.float32],
         full_mask_tensor: npt.NDArray[np.bool_],
         padding: int,

From f21eccd15493eb401c17766dca229312176b0bd8 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 10:05:40 +0000
Subject: [PATCH 14/34] WIP: Handle ImageGrainCrops between run_grains and
 run_grainstats

---
 topostats/processing.py | 35 ++++++++++++++---------------------
 1 file changed, 14 insertions(+), 21 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index c234e6232a..2ff8736e66 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -135,7 +135,7 @@ def run_grains(  # noqa: C901
     core_out_path: Path,
     plotting_config: dict,
     grains_config: dict,
-) -> ImageGrainCrops | None:
+) -> ImageGrainCrops:
     """
     Identify grains (molecules) and optionally plots the results.
 
@@ -184,6 +184,7 @@ def run_grains(  # noqa: C901
             LOGGER.error(
                 f"[{filename}] : An error occurred during grain finding, skipping following steps.", exc_info=e
             )
+            raise e
         else:
             for direction, region_props in grains.region_properties.items():
                 if len(region_props) == 0:
@@ -199,9 +200,6 @@ def run_grains(  # noqa: C901
                         grain_out_path_direction.mkdir(parents=True, exist_ok=True)
                         LOGGER.debug(f"[{filename}] : Target grain directory created : {grain_out_path_direction}")
                     for plot_name, array in image_arrays.items():
-                        if len(array.shape) == 3:
-                            # Use the DNA class mask from the tensor. Hardcoded to 1 as this implementation is not yet generalised.
-                            array = array[:, :, 1]
                         LOGGER.debug(f"[{filename}] : Plotting {plot_name} image")
                         plotting_config["plot_dict"][plot_name]["output_dir"] = grain_out_path_direction
                         Images(
@@ -215,10 +213,9 @@ def run_grains(  # noqa: C901
                         region_properties=grains.region_properties[direction],
                     ).plot_and_save()
                     plotting_config["plot_dict"]["coloured_boxes"]["output_dir"] = grain_out_path_direction
-                    # hard code to class index 1, as this implementation is not yet generalised.
                     Images(
-                        data=np.zeros_like(grains.directions[direction]["labelled_regions_02"][:, :, 1]),
-                        masked_array=grains.directions[direction]["labelled_regions_02"][:, :, 1],
+                        data=np.zeros_like(grains.directions[direction]["labelled_regions_02"]),
+                        masked_array=grains.directions[direction]["labelled_regions_02"],
                         **plotting_config["plot_dict"]["coloured_boxes"],
                         region_properties=grains.region_properties[direction],
                     ).plot_and_save()
@@ -229,7 +226,7 @@ def run_grains(  # noqa: C901
                     Images(
                         image,
                         filename=f"{filename}_{direction}_masked",
-                        masked_array=grains.directions[direction]["removed_small_objects"][:, :, 1].astype(bool),
+                        masked_array=grains.directions[direction]["removed_small_objects"].astype(bool),
                         **plotting_config["plot_dict"][plot_name],
                         region_properties=grains.region_properties[direction],
                     ).plot_and_save()
@@ -237,26 +234,22 @@ def run_grains(  # noqa: C901
             else:
                 # Otherwise, return None and warn that plotting is disabled for grain finding images
                 LOGGER.info(f"[{filename}] : Plotting disabled for Grain Finding Images")
-            grain_masks = {}
-            for direction in grains.directions:
-                grain_masks[direction] = grains.directions[direction]["labelled_regions_02"]
             LOGGER.info(f"[{filename}] : Grain Finding stage completed successfully.")
-            return grain_masks
+            return grains.image_grain_crops
     # Otherwise, return None and warn grainstats is disabled
     LOGGER.info(f"[{filename}] Detection of grains disabled, GrainStats will not be run.")
-    return None
+    return ImageGrainCrops(above=None, below=None)
 
 
 def run_grainstats(
-    image: npt.NDArray,
-    pixel_to_nm_scaling: float,
     image_grain_crops: ImageGrainCrops,
+    pixel_to_nm_scaling: float,
     filename: str,
     basename: Path,
     grainstats_config: dict,
     plotting_config: dict,
     grain_out_path: Path,
-):
+) -> tuple[pd.DataFrame, dict]:
     """
     Calculate grain statistics for an image and optionally plots the results.
 
@@ -351,10 +344,11 @@ def run_grainstats(
             LOGGER.info(f"[{filename}] : Calculated grainstats for {len(grainstats_df)} grains.")
             LOGGER.info(f"[{filename}] : Grainstats stage completed successfully.")
             return grainstats_df, height_profiles_dict
-        except Exception:
+        except Exception as e:
             LOGGER.info(
                 f"[{filename}] : Errors occurred whilst calculating grain statistics. Returning empty dataframe."
             )
+            raise e
             return create_empty_dataframe(column_set="grainstats", index_col="grain_number"), height_profiles_dict
     else:
         LOGGER.info(
@@ -1007,7 +1001,7 @@ def process_scan(
     )
 
     # Find Grains :
-    grain_dict = run_grains(
+    image_grain_crops = run_grains(
         image=topostats_object["image_flattened"],
         pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
         filename=topostats_object["filename"],
@@ -1017,12 +1011,11 @@ def process_scan(
         grains_config=grains_config,
     )
 
-    if "above" in topostats_object["grain_masks"].keys() or "below" in topostats_object["grain_masks"].keys():
+    if image_grain_crops.above is not None or image_grain_crops.below is not None:
         # Grainstats :
         grainstats_df, height_profiles = run_grainstats(
-            image=topostats_object["image_flattened"],
+            image_grain_crops=image_grain_crops,
             pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
-            grain_dict=grain_dict,
             filename=topostats_object["filename"],
             basename=topostats_object["img_path"],
             grainstats_config=grainstats_config,

From 63d0003673c1a4afef9a88bd5fdd70227b66fab0 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 12:04:11 +0000
Subject: [PATCH 15/34] WIP: Graintstats handles ImageGrainCrops

---
 topostats/grainstats.py | 37 +++++++++++++++++++++----------------
 1 file changed, 21 insertions(+), 16 deletions(-)

diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index 56b4b49dce..5f08f1031f 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -200,13 +200,12 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
             Consists of a pd.DataFrame containing all the grain stats that have been calculated for the labelled image
             and a list of dictionaries containing grain data to be plotted.
         """
-        grains_plot_data = []
-        all_height_profiles = {}
+        all_height_profiles: dict[int, npt.NDArray] = {}
         if len(self.grain_crops) == 0:
             LOGGER.warning(
                 f"[{self.image_name}] : No grain crops for this image, grain statistics can not be calculated."
             )
-            return pd.DataFrame(columns=GRAIN_STATS_COLUMNS), grains_plot_data, all_height_profiles
+            return pd.DataFrame(columns=GRAIN_STATS_COLUMNS), all_height_profiles
 
         grainstats_rows: list[dict] = []
 
@@ -228,29 +227,32 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
                 all_height_profiles[grain_index][class_index] = {}
 
                 class_mask = mask[:, :, class_index]
-
+                labelled_class_mask = skimage_measure.label(class_mask)
                 # Split the class into connected components
-                class_mask_regionprops = skimage_measure.regionprops(skimage_measure.label(class_mask))
+                class_mask_regionprops = skimage_measure.regionprops(labelled_class_mask)
 
                 # Iterate over all the sub_grains in the class
                 for subgrain_index, subgrain_region in enumerate(class_mask_regionprops):
-
-                    subgrain_mask = subgrain_region.image
+                    # Remove all but the current subgrain from the mask
+                    subgrain_only_mask = class_mask * (labelled_class_mask == subgrain_region.label)
+                    # Create a masked image of the subgrain
                     subgrain_mask_image = np.ma.masked_array(
-                        image, mask=np.invert(subgrain_mask), fill_value=np.nan
+                        image, mask=np.invert(subgrain_only_mask), fill_value=np.nan
                     ).filled()
 
+                    # Shape of the subgrain region with no padding and not necessarily square, more accurate measure of
+                    # the bounding box size
+                    subgrain_tight_shape = subgrain_region.image.shape
                     # Skip subgrain if too small to calculate stats for
-                    if min(subgrain_mask.shape) < 5:
+                    if min(subgrain_tight_shape) < 5:
                         LOGGER.debug(
                             f"[{self.image_name}] : Skipping subgrain due to being too small "
-                            f"(size: {subgrain_mask.shape}) to calculate stats for."
+                            f"(size: {subgrain_tight_shape}) to calculate stats for."
                         )
 
                     # Calculate all the stats
-
-                    points = self.calculate_points(subgrain_mask)
-                    edges = self.calculate_edges(subgrain_mask, edge_detection_method=self.edge_detection_method)
+                    points = self.calculate_points(subgrain_only_mask)
+                    edges = self.calculate_edges(subgrain_only_mask, edge_detection_method=self.edge_detection_method)
                     radius_stats = self.calculate_radius_stats(edges, points)
                     # hull, hull_indices, hull_simplexes = self.convex_hull(edges, output_grain)
                     _, _, hull_simplexes = self.convex_hull(edges, output_grain)
@@ -280,7 +282,7 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
 
                     if self.extract_height_profile:
                         all_height_profiles[grain_index][class_index][subgrain_index] = (
-                            height_profiles.interpolate_height_profile(img=image, mask=subgrain_mask)
+                            height_profiles.interpolate_height_profile(img=image, mask=subgrain_only_mask)
                         )
                         LOGGER.debug(f"[{self.image_name}] : Height profiles extracted.")
 
@@ -322,13 +324,16 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
                     grainstats_rows.append(stats)
 
         # Check if the dataframe is empty
-        if grainstats_df.empty:
+        if len(grainstats_rows) == 0:
             grainstats_df = create_empty_dataframe()
+        else:
+            # Create a dataframe from the list of dictionaries
+            grainstats_df = pd.DataFrame(grainstats_rows)
 
         grainstats_df.index.name = "grain_number"
         grainstats_df["image"] = self.image_name
 
-        return grainstats_df, grains_plot_data, all_height_profiles
+        return grainstats_df, all_height_profiles
 
     @staticmethod
     def calculate_points(grain_mask: npt.NDArray) -> list:

From 02ecf43bdc2d62e261e7f36fb90465cc3cca6af1 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 12:22:43 +0000
Subject: [PATCH 16/34] Fix: grainstats: process scan no longer needing grain
 plots returned

---
 topostats/processing.py | 15 +--------------
 1 file changed, 1 insertion(+), 14 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index 2ff8736e66..d9f2ccb5b2 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -310,23 +310,10 @@ def run_grainstats(
                     plot_opts=grain_plot_dict,
                     **grainstats_config,
                 )
-                grainstats_dict[direction], grains_plot_data, height_profiles_dict[direction] = (
+                grainstats_dict[direction], height_profiles_dict[direction] = (
                     grainstats_calculator.calculate_stats()
                 )
                 grainstats_dict[direction]["threshold"] = direction
-                # Plot grains if required
-                if plotting_config["image_set"] == "all":
-                    LOGGER.info(f"[{filename}] : Plotting grain images for direction: {direction}.")
-                    for plot_data in grains_plot_data:
-                        LOGGER.debug(
-                            f"[{filename}] : Plotting grain image {plot_data['filename']} for direction: {direction}."
-                        )
-                        Images(
-                            data=plot_data["data"],
-                            output_dir=plot_data["output_dir"],
-                            filename=plot_data["filename"],
-                            **plotting_config["plot_dict"][plot_data["name"]],
-                        ).plot_and_save()
             # Create results dataframe from above and below results
             # Appease pylint and ensure that grainstats_df is always created
             grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")

From 651e89caa96b8a4284c44df2c2c3f8ba782cdd3f Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 22 Nov 2024 13:07:05 +0000
Subject: [PATCH 17/34] WIP: Begin grains > disordered_tracing pipeline
 overhaul

---
 topostats/processing.py                 | 26 +++-----
 topostats/tracing/disordered_tracing.py | 86 ++++++++++++++-----------
 2 files changed, 57 insertions(+), 55 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index d9f2ccb5b2..33debe70f8 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -299,10 +299,8 @@ def run_grainstats(
             for direction, grain_crops_direction in image_grain_crops.__dict__.items():
                 if grain_crops_direction is None:
                     continue
-                grain_crops = grain_crops_direction.crops
-
                 grainstats_calculator = GrainStats(
-                    grain_crops=grain_crops,
+                    grain_crops=grain_crops_direction.crops,
                     pixel_to_nanometre_scaling=pixel_to_nm_scaling,
                     direction=direction,
                     base_output_dir=grain_out_path,
@@ -310,9 +308,7 @@ def run_grainstats(
                     plot_opts=grain_plot_dict,
                     **grainstats_config,
                 )
-                grainstats_dict[direction], height_profiles_dict[direction] = (
-                    grainstats_calculator.calculate_stats()
-                )
+                grainstats_dict[direction], height_profiles_dict[direction] = grainstats_calculator.calculate_stats()
                 grainstats_dict[direction]["threshold"] = direction
             # Create results dataframe from above and below results
             # Appease pylint and ensure that grainstats_df is always created
@@ -346,7 +342,7 @@ def run_grainstats(
 
 def run_disordered_tracing(
     image: npt.NDArray,
-    grain_masks: dict,
+    image_grain_crops: ImageGrainCrops,
     pixel_to_nm_scaling: float,
     filename: str,
     basename: str,
@@ -399,25 +395,23 @@ def run_disordered_tracing(
         disordered_trace_grainstats = pd.DataFrame()
         disordered_tracing_stats_image = pd.DataFrame()
         try:
-            # run image using directional grain masks
-            for direction, _ in grain_masks.items():
-                # Check if there are grains
-                assert len(grain_masks[direction].shape) == 3, "Grain masks should be 3D tensors"
-                dna_class_mask = grain_masks[direction][:, :, 1]
-                if np.max(dna_class_mask) == 0:
+
+            grain_crop_direction: GrainCropsDirection
+            for direction, grain_crop_direction in image_grain_crops.__dict__.items():
+
+                if grain_crop_direction is None:
                     LOGGER.warning(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     raise ValueError(f"No grains exist for the {direction} direction")
-                # if grains are found
+
                 (
                     disordered_traces_cropped_data,
                     _disordered_trace_grainstats,
                     disordered_tracing_images,
                     disordered_tracing_stats,
                 ) = trace_image_disordered(
-                    image=image,
-                    grains_mask=dna_class_mask,
+                    grain_crops=grain_crop_direction.crops,
                     filename=filename,
                     pixel_to_nm_scaling=pixel_to_nm_scaling,
                     **disordered_tracing_config,
diff --git a/topostats/tracing/disordered_tracing.py b/topostats/tracing/disordered_tracing.py
index 9acce787cd..aceb29201a 100644
--- a/topostats/tracing/disordered_tracing.py
+++ b/topostats/tracing/disordered_tracing.py
@@ -18,6 +18,7 @@
 from topostats.tracing.pruning import prune_skeleton
 from topostats.tracing.skeletonize import getSkeleton
 from topostats.utils import convolve_skeleton
+from topostats.grains import GrainCrop
 
 LOGGER = logging.getLogger(LOGGER_NAME)
 
@@ -256,8 +257,8 @@ def smooth_mask(
 
 
 def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-locals
-    image: npt.NDArray,
-    grains_mask: npt.NDArray,
+    grain_crops: dict[int, GrainCrop],
+    class_index: int,
     filename: str,
     pixel_to_nm_scaling: float,
     min_skeleton_size: int,
@@ -297,58 +298,63 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
     tuple[dict, pd.DataFrame, dict, pd.DataFrame]
         Binary and integer labeled cropped and full-image masks from skeletonising and pruning the grains in the image.
     """
-    # Check both arrays are the same shape - should this be a test instead, why should this ever occur?
-    if image.shape != grains_mask.shape:
-        raise ValueError(f"Image shape ({image.shape}) and Mask shape ({grains_mask.shape}) should match.")
+    # # Check both arrays are the same shape - should this be a test instead, why should this ever occur?
+    # if image.shape != grains_mask.shape:
+    #     raise ValueError(f"Image shape ({image.shape}) and Mask shape ({grains_mask.shape}) should match.")
 
-    cropped_images, cropped_masks, bboxs = prep_arrays(image, grains_mask, pad_width)
-    n_grains = len(cropped_images)
-    img_base = np.zeros_like(image)
+    # cropped_images, cropped_masks, bboxs = prep_arrays(image, grains_mask, pad_width)
+    # n_grains = len(cropped_images)
+    # img_base = np.zeros_like(image)
     disordered_trace_crop_data = {}
     grainstats_additions = {}
-    disordered_tracing_stats = pd.DataFrame()
+    # disordered_tracing_stats = pd.DataFrame()
 
     # want to get each cropped image, use some anchor coords to match them onto the image,
     #   and compile all the grain images onto a single image
-    all_images = {
-        "smoothed_grain": img_base.copy(),
-        "skeleton": img_base.copy(),
-        "pruned_skeleton": img_base.copy(),
-        "branch_indexes": img_base.copy(),
-        "branch_types": img_base.copy(),
-    }
-
-    LOGGER.info(f"[{filename}] : Calculating Disordered Tracing statistics for {n_grains} grains...")
-
-    for cropped_image_index, cropped_image in cropped_images.items():
+    # all_images = {
+    #     "smoothed_grain": img_base.copy(),
+    #     "skeleton": img_base.copy(),
+    #     "pruned_skeleton": img_base.copy(),
+    #     "branch_indexes": img_base.copy(),
+    #     "branch_types": img_base.copy(),
+    # }
+
+    # LOGGER.info(f"[{filename}] : Calculating Disordered Tracing statistics for {n_grains} grains...")
+
+    # for cropped_image_index, cropped_image in cropped_images.items():
+    number_of_grains = len(grain_crops)
+    for grain_number, grain_crop in grain_crops.items():
         try:
-            cropped_mask = cropped_masks[cropped_image_index]
+            grain_crop_tensor = grain_crop.mask
+            grain_crop_class_mask = grain_crop_tensor[:, :, class_index]
+            grain_crop_image = grain_crop.image
+
             disordered_trace_images = disordered_trace_grain(
-                cropped_image=cropped_image,
-                cropped_mask=cropped_mask,
+                cropped_image=grain_crop_image,
+                cropped_mask=grain_crop_class_mask,
                 pixel_to_nm_scaling=pixel_to_nm_scaling,
                 mask_smoothing_params=mask_smoothing_params,
                 skeletonisation_params=skeletonisation_params,
                 pruning_params=pruning_params,
                 filename=filename,
                 min_skeleton_size=min_skeleton_size,
-                n_grain=cropped_image_index,
+                n_grain=grain_number,
             )
-            LOGGER.debug(f"[{filename}] : Disordered Traced grain {cropped_image_index + 1} of {n_grains}")
+            LOGGER.debug(f"[{filename}] : Disordered Traced grain {grain_number + 1} of {number_of_grains}")
 
             if disordered_trace_images is not None:
                 # obtain segment stats
                 try:
                     skan_skeleton = skan.Skeleton(
-                        np.where(disordered_trace_images["pruned_skeleton"] == 1, cropped_image, 0),
+                        np.where(disordered_trace_images["pruned_skeleton"] == 1, grain_crop_image, 0),
                         spacing=pixel_to_nm_scaling,
                     )
                     skan_df = skan.summarize(skan_skeleton)
-                    skan_df = compile_skan_stats(skan_df, skan_skeleton, cropped_image, filename, cropped_image_index)
+                    skan_df = compile_skan_stats(skan_df, skan_skeleton, grain_crop_image, filename, grain_number)
                     total_branch_length = skan_df["branch_distance"].sum() * 1e-9
                 except ValueError:
                     LOGGER.warning(
-                        f"[{filename}] : Skeleton for grain {cropped_image_index} has been pruned out of existence."
+                        f"[{filename}] : Skeleton for grain {grain_number} has been pruned out of existence."
                     )
                     total_branch_length = 0
                     skan_df = pd.DataFrame()
@@ -357,34 +363,36 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
 
                 # obtain stats
                 conv_pruned_skeleton = convolve_skeleton(disordered_trace_images["pruned_skeleton"])
-                grainstats_additions[cropped_image_index] = {
+                grainstats_additions[grain_number] = {
                     "image": filename,
-                    "grain_number": cropped_image_index,
+                    "grain_number": grain_number,
                     "grain_endpoints": np.int64((conv_pruned_skeleton == 2).sum()),
                     "grain_junctions": np.int64((conv_pruned_skeleton == 3).sum()),
                     "total_branch_lengths": total_branch_length,
                 }
 
-                # remap the cropped images back onto the original
-                for image_name, full_image in all_images.items():
-                    crop = disordered_trace_images[image_name]
-                    bbox = bboxs[cropped_image_index]
-                    full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
-                disordered_trace_crop_data[f"grain_{cropped_image_index}"] = disordered_trace_images
-                disordered_trace_crop_data[f"grain_{cropped_image_index}"]["bbox"] = bboxs[cropped_image_index]
+                # # remap the cropped images back onto the original
+                # for image_name, full_image in all_images.items():
+                #     crop = disordered_trace_images[image_name]
+                #     bbox = bboxs[grain_number]
+                #     full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
+                disordered_trace_crop_data[f"grain_{grain_number}"] = disordered_trace_images
+                disordered_trace_crop_data[f"grain_{grain_number}"]["bbox"] = grain_crop.bbox
+                disordered_trace_crop_data[f"grain_{grain_number}"]["padding"] = grain_crop.padding
 
         # when skel too small, pruned to 0's, skan -> ValueError -> skipped
         except Exception as e:  # pylint: disable=broad-exception-caught
             LOGGER.error(  # pylint: disable=logging-not-lazy
                 f"[{filename}] : Disordered tracing of grain "
-                f"{cropped_image_index} failed. Consider raising an issue on GitHub. Error: ",
+                f"{grain_number} failed. Consider raising an issue on GitHub. Error: ",
                 exc_info=e,
             )
 
         # convert stats dict to dataframe
         grainstats_additions_df = pd.DataFrame.from_dict(grainstats_additions, orient="index")
 
-    return disordered_trace_crop_data, grainstats_additions_df, all_images, disordered_tracing_stats
+    # return disordered_trace_crop_data, grainstats_additions_df, all_images, disordered_tracing_stats
+    return disordered_trace_crop_data, grainstats_additions_df, disordered_tracing_stats
 
 
 def compile_skan_stats(

From a03524ebd30933ec92f474b42355d8ed0b71ead0 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 22 Nov 2024 17:22:03 +0000
Subject: [PATCH 18/34] [pre-commit.ci] Fixing issues with pre-commit

---
 topostats/grains.py                     | 12 +++---------
 topostats/grainstats.py                 |  6 ++----
 topostats/processing.py                 | 10 +++-------
 topostats/tracing/disordered_tracing.py |  4 ++--
 4 files changed, 10 insertions(+), 22 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 8a47763abe..6942373116 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -6,7 +6,7 @@
 import logging
 import sys
 from collections import defaultdict
-from dataclasses import dataclass, field
+from dataclasses import dataclass
 
 import keras
 import numpy as np
@@ -22,8 +22,6 @@
 from topostats.unet_masking import (
     make_bounding_box_square,
     pad_bounding_box,
-    pad_crop,
-    make_crop_square,
     predict_unet,
 )
 from topostats.utils import _get_mask, get_thresholds
@@ -937,9 +935,7 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:]
-                for vetting_criteria in class_size_thresholds
-                if vetting_criteria[0] == class_index
+                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1666,9 +1662,7 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(
-                flat_region.bbox
-            )  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(
diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index 5f08f1031f..f21d2870c9 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -15,10 +15,10 @@
 import skimage.measure as skimage_measure
 import skimage.morphology as skimage_morphology
 
+from topostats.grains import GrainCrop
 from topostats.logs.logs import LOGGER_NAME
 from topostats.measure import feret, height_profiles
 from topostats.utils import create_empty_dataframe
-from topostats.grains import GrainCrop
 
 # pylint: disable=too-many-lines
 # pylint: disable=too-many-instance-attributes
@@ -584,9 +584,7 @@ def sort_points(self, points: list) -> list:
         # relative to that and _then_ sort it.
         # pivot_angles = self.get_angle(points, self.start_point)
         # Recursively sort the arrays until each point is sorted
-        return (
-            self.sort_points(smaller) + sorted(equal, key=self.calculate_squared_distance) + self.sort_points(larger)
-        )
+        return self.sort_points(smaller) + sorted(equal, key=self.calculate_squared_distance) + self.sort_points(larger)
         # Return sorted array where equal angle points are sorted by distance
 
     def get_start_point(self, edges: npt.NDArray) -> None:
diff --git a/topostats/processing.py b/topostats/processing.py
index eee24835c6..edca008975 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -13,7 +13,7 @@
 
 from topostats import __version__
 from topostats.filters import Filters
-from topostats.grains import Grains, ImageGrainCrops, GrainCropsDirection
+from topostats.grains import GrainCropsDirection, Grains, ImageGrainCrops
 from topostats.grainstats import GrainStats
 from topostats.io import get_out_path, save_topostats_file
 from topostats.logs.logs import LOGGER_NAME
@@ -802,9 +802,7 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(
-                        column_set="mol_statistics", index_col="molecule_number"
-                    )
+                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1178,9 +1176,7 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(
-            column_set="disordered_tracing_statistics", index_col="index"
-        )
+        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
         height_profiles = {}
 
     # Get image statistics
diff --git a/topostats/tracing/disordered_tracing.py b/topostats/tracing/disordered_tracing.py
index 95428b6043..ecfcadf1db 100644
--- a/topostats/tracing/disordered_tracing.py
+++ b/topostats/tracing/disordered_tracing.py
@@ -14,11 +14,11 @@
 from skimage import filters
 from skimage.morphology import label
 
+from topostats.grains import GrainCrop
 from topostats.logs.logs import LOGGER_NAME
 from topostats.tracing.pruning import prune_skeleton
 from topostats.tracing.skeletonize import getSkeleton
 from topostats.utils import convolve_skeleton
-from topostats.grains import GrainCrop
 
 LOGGER = logging.getLogger(LOGGER_NAME)
 
@@ -372,7 +372,7 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
                 }
 
             # remap the cropped images back onto the original
-            #for image_name, full_image in all_images.items():
+            # for image_name, full_image in all_images.items():
             #    crop = disordered_trace_images[image_name]
             #    bbox = bboxs[cropped_image_index]
             #    full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]

From 8d215fdea9d3455ff22e609d9e1d4048316553b3 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 26 Nov 2024 10:40:26 +0000
Subject: [PATCH 19/34] WIP: grains > disorderd_tracing pipeline | fix typing
 and remove whole image plotting

---
 topostats/processing.py                 | 61 ++++++++++++++-----------
 topostats/tracing/disordered_tracing.py | 20 +++++---
 2 files changed, 47 insertions(+), 34 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index edca008975..afb1543062 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -344,17 +344,17 @@ def run_grainstats(
 
 
 def run_disordered_tracing(
-    image: npt.NDArray,
+    # image: npt.NDArray,
     image_grain_crops: ImageGrainCrops,
     pixel_to_nm_scaling: float,
     filename: str,
     basename: str,
-    core_out_path: Path,
-    tracing_out_path: Path,
+    # core_out_path: Path,
+    # tracing_out_path: Path,
     disordered_tracing_config: dict,
-    plotting_config: dict,
+    # plotting_config: dict,
     grainstats_df: pd.DataFrame = None,
-) -> dict:
+) -> tuple[dict | None, pd.DataFrame, pd.DataFrame]:
     """
     Skeletonise and prune grains, adding results to statistics data frames and optionally plot results.
 
@@ -411,7 +411,7 @@ def run_disordered_tracing(
                 (
                     disordered_traces_cropped_data,
                     _disordered_trace_grainstats,
-                    disordered_tracing_images,
+                    # disordered_tracing_images,
                     disordered_tracing_stats,
                 ) = trace_image_disordered(
                     grain_crops=grain_crop_direction.crops,
@@ -428,20 +428,23 @@ def run_disordered_tracing(
                 # append direction results to dict
                 disordered_traces[direction] = disordered_traces_cropped_data
                 # save plots
-                Images(
-                    image,
-                    masked_array=disordered_tracing_images.pop("pruned_skeleton"),
-                    output_dir=core_out_path,
-                    filename=f"{filename}_{direction}_disordered_trace",
-                    **plotting_config["plot_dict"]["pruned_skeleton"],
-                ).plot_and_save()
-                for plot_name, image_value in disordered_tracing_images.items():
-                    Images(
-                        image,
-                        masked_array=image_value,
-                        output_dir=tracing_out_path / direction,
-                        **plotting_config["plot_dict"][plot_name],
-                    ).plot_and_save()
+
+                # PLOT INDIVIDUALLY INSIDE THE FUNCTION??
+
+                # Images(
+                #     image,
+                #     masked_array=disordered_tracing_images.pop("pruned_skeleton"),
+                #     output_dir=core_out_path,
+                #     filename=f"{filename}_{direction}_disordered_trace",
+                #     **plotting_config["plot_dict"]["pruned_skeleton"],
+                # ).plot_and_save()
+                # for plot_name, image_value in disordered_tracing_images.items():
+                #     Images(
+                #         image,
+                #         masked_array=image_value,
+                #         output_dir=tracing_out_path / direction,
+                #         **plotting_config["plot_dict"][plot_name],
+                #     ).plot_and_save()
             # merge grainstats data with other dataframe
             resultant_grainstats = (
                 pd.merge(grainstats_df, disordered_trace_grainstats, on=["image", "threshold", "grain_number"])
@@ -802,7 +805,9 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
+                    splining_molstats = create_empty_dataframe(
+                        column_set="mol_statistics", index_col="molecule_number"
+                    )
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1101,16 +1106,16 @@ def process_scan(
 
         # Disordered Tracing
         disordered_traces_data, grainstats_df, disordered_tracing_stats = run_disordered_tracing(
-            image=topostats_object["image_flattened"],
-            grain_dict=grain_dict,
+            # image=topostats_object["image_flattened"],
+            image_grain_crops=image_grain_crops,
             pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
             filename=topostats_object["filename"],
             basename=topostats_object["img_path"],
-            core_out_path=core_out_path,
-            tracing_out_path=tracing_out_path,
+            # core_out_path=core_out_path,
+            # tracing_out_path=tracing_out_path,
             disordered_tracing_config=disordered_tracing_config,
             grainstats_df=grainstats_df,
-            plotting_config=plotting_config,
+            # plotting_config=plotting_config,
         )
         topostats_object["disordered_traces"] = disordered_traces_data
 
@@ -1176,7 +1181,9 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
+        disordered_tracing_stats = create_empty_dataframe(
+            column_set="disordered_tracing_statistics", index_col="index"
+        )
         height_profiles = {}
 
     # Get image statistics
diff --git a/topostats/tracing/disordered_tracing.py b/topostats/tracing/disordered_tracing.py
index ecfcadf1db..97ae32de60 100644
--- a/topostats/tracing/disordered_tracing.py
+++ b/topostats/tracing/disordered_tracing.py
@@ -266,7 +266,7 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
     skeletonisation_params: dict,
     pruning_params: dict,
     pad_width: int = 1,
-) -> tuple[dict, pd.DataFrame, dict, pd.DataFrame]:
+) -> tuple[dict[str, dict], pd.DataFrame, pd.DataFrame]:
     """
     Processor function for tracing image.
 
@@ -307,7 +307,7 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
     # img_base = np.zeros_like(image)
     disordered_trace_crop_data = {}
     grainstats_additions = {}
-    # disordered_tracing_stats = pd.DataFrame()
+    disordered_tracing_stats = pd.DataFrame()
 
     # want to get each cropped image, use some anchor coords to match them onto the image,
     #   and compile all the grain images onto a single image
@@ -349,8 +349,14 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
                         np.where(disordered_trace_images["pruned_skeleton"] == 1, grain_crop_image, 0),
                         spacing=pixel_to_nm_scaling,
                     )
-                    skan_df = skan.summarize(skan_skeleton)
-                    skan_df = compile_skan_stats(skan_df, skan_skeleton, grain_crop_image, filename, grain_number)
+                    skan_df = skan.summarize(skel=skan_skeleton)
+                    skan_df = compile_skan_stats(
+                        skan_df=skan_df,
+                        skan_skeleton=skan_skeleton,
+                        image=grain_crop_image,
+                        filename=filename,
+                        grain_number=grain_number,
+                    )
                     total_branch_length = skan_df["branch_distance"].sum() * 1e-9
                 except ValueError:
                     LOGGER.warning(
@@ -376,9 +382,9 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
             #    crop = disordered_trace_images[image_name]
             #    bbox = bboxs[cropped_image_index]
             #    full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
-            disordered_trace_crop_data[f"grain_{cropped_image_index}"] = disordered_trace_images
-            disordered_trace_crop_data[f"grain_{cropped_image_index}"]["bbox"] = grain_crop.bbox
-            disordered_trace_crop_data[f"grain_{cropped_image_index}"]["pad_width"] = grain_crop.padding
+            disordered_trace_crop_data[f"grain_{grain_number}"] = disordered_trace_images
+            disordered_trace_crop_data[f"grain_{grain_number}"]["bbox"] = grain_crop.bbox
+            disordered_trace_crop_data[f"grain_{grain_number}"]["pad_width"] = grain_crop.padding
 
         # when skel too small, pruned to 0's, skan -> ValueError -> skipped
         except Exception as e:  # pylint: disable=broad-exception-caught

From 3a188802abcb1432c39c328bf7618cbe63f731bb Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 26 Nov 2024 10:41:14 +0000
Subject: [PATCH 20/34] [pre-commit.ci] Fixing issues with pre-commit

---
 topostats/processing.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index afb1543062..7fe9e4c5d6 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -805,9 +805,7 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(
-                        column_set="mol_statistics", index_col="molecule_number"
-                    )
+                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1181,9 +1179,7 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(
-            column_set="disordered_tracing_statistics", index_col="index"
-        )
+        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
         height_profiles = {}
 
     # Get image statistics

From 7781c867fd777f34ddf6212f8d315a5d34b9e18d Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Wed, 4 Dec 2024 10:59:33 +0000
Subject: [PATCH 21/34] Add: class index to disordered tracing config

---
 topostats/default_config.yaml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/topostats/default_config.yaml b/topostats/default_config.yaml
index ad7b1e7360..73cacbf349 100644
--- a/topostats/default_config.yaml
+++ b/topostats/default_config.yaml
@@ -65,6 +65,7 @@ grainstats:
   extract_height_profile: true # Extract height profiles along maximum feret of molecules
 disordered_tracing:
   run: true # Options : true, false
+  class_index: 1 # The class index to trace. This is the class index of the grains.
   min_skeleton_size: 10 # Minimum number of pixels in a skeleton for it to be retained.
   pad_width: 1 # Pixels to pad grains by when tracing
   mask_smoothing_params:

From 80b711cac2238c5ac04966f070cac4c916b79847 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 10 Dec 2024 16:40:28 +0000
Subject: [PATCH 22/34] [WIP] Fix: Attempt to fix grain_number double index
 issue

---
 topostats/grainstats.py                 | 3 ++-
 topostats/tracing/disordered_tracing.py | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/topostats/grainstats.py b/topostats/grainstats.py
index f21d2870c9..45b576ec38 100644
--- a/topostats/grainstats.py
+++ b/topostats/grainstats.py
@@ -330,7 +330,8 @@ def calculate_stats(self) -> tuple(pd.DataFrame, dict):
             # Create a dataframe from the list of dictionaries
             grainstats_df = pd.DataFrame(grainstats_rows)
 
-        grainstats_df.index.name = "grain_number"
+        # Change the index column from the arbitrary one to the grain number
+        grainstats_df.set_index("grain_number", inplace=True)
         grainstats_df["image"] = self.image_name
 
         return grainstats_df, all_height_profiles
diff --git a/topostats/tracing/disordered_tracing.py b/topostats/tracing/disordered_tracing.py
index 97ae32de60..344336f68c 100644
--- a/topostats/tracing/disordered_tracing.py
+++ b/topostats/tracing/disordered_tracing.py
@@ -371,7 +371,6 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
                 conv_pruned_skeleton = convolve_skeleton(disordered_trace_images["pruned_skeleton"])
                 grainstats_additions[grain_number] = {
                     "image": filename,
-                    "grain_number": grain_number,
                     "grain_endpoints": np.int64((conv_pruned_skeleton == 2).sum()),
                     "grain_junctions": np.int64((conv_pruned_skeleton == 3).sum()),
                     "total_branch_lengths": total_branch_length,
@@ -396,6 +395,8 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
 
         # convert stats dict to dataframe
         grainstats_additions_df = pd.DataFrame.from_dict(grainstats_additions, orient="index")
+        # Set the name of the index column to be the grain number
+        grainstats_additions_df.index.name = "grain_number"
 
     # return disordered_trace_crop_data, grainstats_additions_df, all_images, disordered_tracing_stats
     return disordered_trace_crop_data, grainstats_additions_df, disordered_tracing_stats

From 4c6f0f3879674a8f0e99432224487289dcdf2864 Mon Sep 17 00:00:00 2001
From: Sylvia Whittle <86117496+SylviaWhittle@users.noreply.github.com>
Date: Wed, 11 Dec 2024 16:24:45 +0000
Subject: [PATCH 23/34] remove raising error on empty direction

Co-authored-by: Neil Shephard <n.shephard@sheffield.ac.uk>
---
 topostats/processing.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/topostats/processing.py b/topostats/processing.py
index 7fe9e4c5d6..f2d8a30ed1 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -406,7 +406,7 @@ def run_disordered_tracing(
                     LOGGER.warning(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
-                    raise ValueError(f"No grains exist for the {direction} direction")
+                    
 
                 (
                     disordered_traces_cropped_data,

From e061001a8c1adf1dafae3d0619721bc1e3cd5d0c Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 17 Dec 2024 10:11:55 +0000
Subject: [PATCH 24/34] Add padding to the grains section of config and remove
 from unet section.

---
 topostats/default_config.yaml | 2 +-
 topostats/validation.py       | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/topostats/default_config.yaml b/topostats/default_config.yaml
index 73cacbf349..df0b8a979d 100644
--- a/topostats/default_config.yaml
+++ b/topostats/default_config.yaml
@@ -30,6 +30,7 @@ filter:
 grains:
   run: true # Options : true, false
   # Thresholding by height
+  padding: 1 # Padding to apply to grains. Needs to be at least 1, more padding may help with unets.
   threshold_method: std_dev # Options : std_dev, otsu, absolute, unet
   otsu_threshold_multiplier: 1.0
   threshold_std_dev:
@@ -47,7 +48,6 @@ grains:
   remove_edge_intersecting_grains: true # Whether or not to remove grains that touch the image border
   unet_config:
     model_path: null # Path to a trained U-Net model
-    grain_crop_padding: 2 # Padding to apply to the grain crop bounding box
     upper_norm_bound: 5.0 # Upper bound for normalisation of input data. This should be slightly higher than the maximum desired / expected height of grains.
     lower_norm_bound: -1.0 # Lower bound for normalisation of input data. This should be slightly lower than the minimum desired / expected height of the background.
   vetting:
diff --git a/topostats/validation.py b/topostats/validation.py
index 0f93f7b8de..2f11f86f39 100644
--- a/topostats/validation.py
+++ b/topostats/validation.py
@@ -116,6 +116,7 @@ def validate_config(config: dict, schema: Schema, config_type: str) -> None:
                 False,
                 error="Invalid value in config for grains.run, valid values are 'True' or 'False'",
             ),
+            "padding": int,
             "smallest_grain_size_nm2": lambda n: n > 0.0,
             "threshold_method": Or(
                 "absolute",
@@ -182,7 +183,6 @@ def validate_config(config: dict, schema: Schema, config_type: str) -> None:
             ),
             "unet_config": {
                 "model_path": Or(None, str),
-                "grain_crop_padding": int,
                 "upper_norm_bound": float,
                 "lower_norm_bound": float,
             },

From 8855707cc885f5325b84393ee6ecf4343e950dc8 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Tue, 17 Dec 2024 12:23:10 +0000
Subject: [PATCH 25/34] [WIP]: Attempt to fix data passing errors between
 processing and disordered tracing

---
 topostats/default_config.yaml           |  2 +-
 topostats/grains.py                     | 10 +++++--
 topostats/processing.py                 | 15 ++++++----
 topostats/tracing/disordered_tracing.py | 38 +++++++++++++------------
 topostats/validation.py                 |  2 +-
 5 files changed, 39 insertions(+), 28 deletions(-)

diff --git a/topostats/default_config.yaml b/topostats/default_config.yaml
index df0b8a979d..3478efba52 100644
--- a/topostats/default_config.yaml
+++ b/topostats/default_config.yaml
@@ -30,7 +30,7 @@ filter:
 grains:
   run: true # Options : true, false
   # Thresholding by height
-  padding: 1 # Padding to apply to grains. Needs to be at least 1, more padding may help with unets.
+  grain_crop_padding: 1 # Padding to apply to grains. Needs to be at least 1, more padding may help with unets.
   threshold_method: std_dev # Options : std_dev, otsu, absolute, unet
   otsu_threshold_multiplier: 1.0
   threshold_std_dev:
diff --git a/topostats/grains.py b/topostats/grains.py
index 6942373116..0bf1ea3255 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -184,6 +184,7 @@ def __init__(
         image: npt.NDArray,
         filename: str,
         pixel_to_nm_scaling: float,
+        grain_crop_padding: int = 1,
         unet_config: dict[str, str | int | float | tuple[int | None, int, int, int] | None] | None = None,
         threshold_method: str | None = None,
         otsu_threshold_multiplier: float | None = None,
@@ -195,7 +196,6 @@ def __init__(
         remove_edge_intersecting_grains: bool = True,
         classes_to_merge: list[tuple[int, int]] | None = None,
         vetting: dict | None = None,
-        grain_crop_padding: int = 1,
     ):
         """
         Initialise the class.
@@ -935,7 +935,9 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
+                vetting_criteria[1:]
+                for vetting_criteria in class_size_thresholds
+                if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1662,7 +1664,9 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(
+                flat_region.bbox
+            )  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(
diff --git a/topostats/processing.py b/topostats/processing.py
index f2d8a30ed1..3a7e4e33f8 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -344,7 +344,7 @@ def run_grainstats(
 
 
 def run_disordered_tracing(
-    # image: npt.NDArray,
+    full_image: npt.NDArray,
     image_grain_crops: ImageGrainCrops,
     pixel_to_nm_scaling: float,
     filename: str,
@@ -406,7 +406,7 @@ def run_disordered_tracing(
                     LOGGER.warning(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
-                    
+                    continue
 
                 (
                     disordered_traces_cropped_data,
@@ -414,6 +414,7 @@ def run_disordered_tracing(
                     # disordered_tracing_images,
                     disordered_tracing_stats,
                 ) = trace_image_disordered(
+                    full_image=full_image,
                     grain_crops=grain_crop_direction.crops,
                     filename=filename,
                     pixel_to_nm_scaling=pixel_to_nm_scaling,
@@ -805,7 +806,9 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
+                    splining_molstats = create_empty_dataframe(
+                        column_set="mol_statistics", index_col="molecule_number"
+                    )
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1104,7 +1107,7 @@ def process_scan(
 
         # Disordered Tracing
         disordered_traces_data, grainstats_df, disordered_tracing_stats = run_disordered_tracing(
-            # image=topostats_object["image_flattened"],
+            full_image=topostats_object["image_flattened"],
             image_grain_crops=image_grain_crops,
             pixel_to_nm_scaling=topostats_object["pixel_to_nm_scaling"],
             filename=topostats_object["filename"],
@@ -1179,7 +1182,9 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
+        disordered_tracing_stats = create_empty_dataframe(
+            column_set="disordered_tracing_statistics", index_col="index"
+        )
         height_profiles = {}
 
     # Get image statistics
diff --git a/topostats/tracing/disordered_tracing.py b/topostats/tracing/disordered_tracing.py
index 344336f68c..592abe96f2 100644
--- a/topostats/tracing/disordered_tracing.py
+++ b/topostats/tracing/disordered_tracing.py
@@ -257,6 +257,7 @@ def smooth_mask(
 
 
 def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-locals
+    full_image: npt.NDArray,
     grain_crops: dict[int, GrainCrop],
     class_index: int,
     filename: str,
@@ -272,7 +273,7 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
 
     Parameters
     ----------
-    image : npt.NDArray
+    full_image : npt.NDArray
         Full image as Numpy Array.
     grains_mask : npt.NDArray
         Full image as Grains that are labelled.
@@ -304,20 +305,20 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
 
     # cropped_images, cropped_masks, bboxs = prep_arrays(image, grains_mask, pad_width)
     # n_grains = len(cropped_images)
-    # img_base = np.zeros_like(image)
+    img_base = np.zeros_like(full_image)
     disordered_trace_crop_data = {}
     grainstats_additions = {}
     disordered_tracing_stats = pd.DataFrame()
 
     # want to get each cropped image, use some anchor coords to match them onto the image,
     #   and compile all the grain images onto a single image
-    # all_images = {
-    #     "smoothed_grain": img_base.copy(),
-    #     "skeleton": img_base.copy(),
-    #     "pruned_skeleton": img_base.copy(),
-    #     "branch_indexes": img_base.copy(),
-    #     "branch_types": img_base.copy(),
-    # }
+    all_images = {
+        "smoothed_grain": img_base.copy(),
+        "skeleton": img_base.copy(),
+        "pruned_skeleton": img_base.copy(),
+        "branch_indexes": img_base.copy(),
+        "branch_types": img_base.copy(),
+    }
 
     # LOGGER.info(f"[{filename}] : Calculating Disordered Tracing statistics for {n_grains} grains...")
 
@@ -329,7 +330,7 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
             grain_crop_class_mask = grain_crop_tensor[:, :, class_index]
             grain_crop_image = grain_crop.image
 
-            disordered_trace_images = disordered_trace_grain(
+            disordered_trace_images: dict | None = disordered_trace_grain(
                 cropped_image=grain_crop_image,
                 cropped_mask=grain_crop_class_mask,
                 pixel_to_nm_scaling=pixel_to_nm_scaling,
@@ -376,14 +377,15 @@ def trace_image_disordered(  # pylint: disable=too-many-arguments,too-many-local
                     "total_branch_lengths": total_branch_length,
                 }
 
-            # remap the cropped images back onto the original
-            # for image_name, full_image in all_images.items():
-            #    crop = disordered_trace_images[image_name]
-            #    bbox = bboxs[cropped_image_index]
-            #    full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
-            disordered_trace_crop_data[f"grain_{grain_number}"] = disordered_trace_images
-            disordered_trace_crop_data[f"grain_{grain_number}"]["bbox"] = grain_crop.bbox
-            disordered_trace_crop_data[f"grain_{grain_number}"]["pad_width"] = grain_crop.padding
+                # remap the cropped images back onto the original, there are many image crops that we want to
+                #  remap back onto the original image so we iterate over them, as passed by the function
+                for image_name, full_image in all_images.items():
+                    crop = disordered_trace_images[image_name]
+                    bbox = grain_crop.bbox
+                    full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop
+                disordered_trace_crop_data[f"grain_{grain_number}"] = disordered_trace_images
+                disordered_trace_crop_data[f"grain_{grain_number}"]["bbox"] = grain_crop.bbox
+                disordered_trace_crop_data[f"grain_{grain_number}"]["pad_width"] = grain_crop.padding
 
         # when skel too small, pruned to 0's, skan -> ValueError -> skipped
         except Exception as e:  # pylint: disable=broad-exception-caught
diff --git a/topostats/validation.py b/topostats/validation.py
index 2f11f86f39..45d3f3bab8 100644
--- a/topostats/validation.py
+++ b/topostats/validation.py
@@ -116,7 +116,7 @@ def validate_config(config: dict, schema: Schema, config_type: str) -> None:
                 False,
                 error="Invalid value in config for grains.run, valid values are 'True' or 'False'",
             ),
-            "padding": int,
+            "grain_crop_padding": int,
             "smallest_grain_size_nm2": lambda n: n > 0.0,
             "threshold_method": Or(
                 "absolute",

From 6547c596478e6b9fdb9d98d9a9400a87f3634d06 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 18 Dec 2024 12:10:42 +0000
Subject: [PATCH 26/34] [pre-commit.ci] Fixing issues with pre-commit

---
 topostats/grains.py     | 8 ++------
 topostats/processing.py | 8 ++------
 2 files changed, 4 insertions(+), 12 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index f8ce199f63..778f753c08 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -940,9 +940,7 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:]
-                for vetting_criteria in class_size_thresholds
-                if vetting_criteria[0] == class_index
+                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1669,9 +1667,7 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(
-                flat_region.bbox
-            )  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(
diff --git a/topostats/processing.py b/topostats/processing.py
index 3f80df390e..d2dd8d3948 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -811,9 +811,7 @@ def run_splining(
                         f"[{filename}] : No grains exist for the {direction} direction. Skipping disordered_tracing for {direction}."
                     )
                     splining_grainstats = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
-                    splining_molstats = create_empty_dataframe(
-                        column_set="mol_statistics", index_col="molecule_number"
-                    )
+                    splining_molstats = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
                     raise ValueError(f"No grains exist for the {direction} direction")
                 # if grains are found
                 (
@@ -1195,9 +1193,7 @@ def process_scan(
     else:
         grainstats_df = create_empty_dataframe(column_set="grainstats", index_col="grain_number")
         molstats_df = create_empty_dataframe(column_set="mol_statistics", index_col="molecule_number")
-        disordered_tracing_stats = create_empty_dataframe(
-            column_set="disordered_tracing_statistics", index_col="index"
-        )
+        disordered_tracing_stats = create_empty_dataframe(column_set="disordered_tracing_statistics", index_col="index")
         height_profiles = {}
 
     # Get image statistics

From 8030c8e5af0437730b381d1a2f98f799e99025ef Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Wed, 18 Dec 2024 15:17:12 +0000
Subject: [PATCH 27/34] Add post init script to validate input data to
 GrainCrop dataclass

---
 topostats/grains.py | 27 ++++++++++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 778f753c08..8fd4b7129d 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -51,7 +51,7 @@ class GrainCrop:
     image : npt.NDArray[np.float32]
         2-D Numpy array of the cropped image.
     mask : npt.NDArray[np.bool_]
-        2-D Numpy array of the cropped mask.
+        3-D Numpy tensor of the cropped mask.
     padding : int
         Padding added to the bounding box of the grain during cropping.
     bbox: tuple[int, int, int, int]
@@ -63,6 +63,31 @@ class GrainCrop:
     padding: int
     bbox: tuple[int, int, int, int]
 
+    def __post_init__(self):
+        """
+        Validate the data makes sense.
+        """
+
+        # If image is not square
+        if self.image.shape[0] != self.image.shape[1]:
+            raise ValueError(f"Image is not square: {self.image.shape}")
+
+        # If first two dimensions of mask are not the same as the image
+        if self.mask.shape[0] != self.image.shape[0] or self.mask.shape[1] != self.image.shape[1]:
+            raise ValueError(f"Mask dimensions do not match image: {self.mask.shape} vs {self.image.shape}")
+
+        if self.padding < 1:
+            raise ValueError(f"Padding must be >= 1, but is {self.padding}")
+
+        if len(self.bbox) != 4:
+            raise ValueError(f"Bounding box must have 4 elements, but has {len(self.bbox)}")
+
+        # if bbox is not square
+        if self.bbox[2] - self.bbox[0] != self.bbox[3] - self.bbox[1]:
+            raise ValueError(
+                f"Bounding box is not square: {self.bbox}, size: {self.bbox[2] - self.bbox[0]} x {self.bbox[3] - self.bbox[1]}"
+            )
+
 
 def validate_full_mask_tensor_shape(array: npt.NDArray[np.bool_]) -> npt.NDArray[np.bool_]:
     """

From b2ec80ddac986e1b83e8d282c9c2afbba54a798e Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Wed, 18 Dec 2024 15:18:49 +0000
Subject: [PATCH 28/34] Fix bug: wrong (old) graincrops returned from unet
 masking

---
 topostats/grains.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 8fd4b7129d..fab3ed8921 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -794,7 +794,6 @@ def improve_grain_segmentation_unet(
         new_graincrops = {}
         for grain_number, graincrop in graincrops.items():
             LOGGER.debug(f"Unet predicting mask for grain {grain_number} of {len(graincrops)}")
-
             # Run the UNet on the region. This is allowed to be a single channel
             # as we can add a background channel afterwards if needed.
             # Remember that this region is cropped from the original image, so it's not
@@ -807,7 +806,6 @@ def improve_grain_segmentation_unet(
                 upper_norm_bound=unet_config["upper_norm_bound"],
                 lower_norm_bound=unet_config["lower_norm_bound"],
             )
-
             assert len(predicted_mask.shape) == 3
             LOGGER.debug(f"Predicted mask shape: {predicted_mask.shape}")
 
@@ -818,7 +816,7 @@ def improve_grain_segmentation_unet(
                 bbox=graincrop.bbox,
             )
 
-        return graincrops
+        return new_graincrops
 
     @staticmethod
     def keep_largest_labelled_region(

From be9f6632fe2acafdef555cd4e1aaef7bc91d65dc Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Wed, 18 Dec 2024 15:21:41 +0000
Subject: [PATCH 29/34] Fix bug: wrong bbox used in construction of GrainCrop
 dataclass causing non square bboxes

---
 topostats/grains.py | 37 ++++++++++++++++++++++++++++++++++---
 1 file changed, 34 insertions(+), 3 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index fab3ed8921..ad6d917535 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -963,7 +963,9 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
+                vetting_criteria[1:]
+                for vetting_criteria in class_size_thresholds
+                if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1690,7 +1692,9 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(
+                flat_region.bbox
+            )  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(
@@ -1711,26 +1715,53 @@ def extract_grains_from_full_image_tensor(
                 image_shape=(full_mask_tensor.shape[0], full_mask_tensor.shape[1]),
             )
 
+            assert (
+                square_flat_bounding_box[0] - square_flat_bounding_box[2]
+                == square_flat_bounding_box[1] - square_flat_bounding_box[3]
+            )
+
+            print(
+                f"square flat bounding box {square_flat_bounding_box} shape {square_flat_bounding_box[2] - square_flat_bounding_box[0], square_flat_bounding_box[3] - square_flat_bounding_box[1]}"
+            )
+
             # Grab image and mask for the cropped region
             grain_cropped_image = image[
                 square_flat_bounding_box[0] : square_flat_bounding_box[2],
                 square_flat_bounding_box[1] : square_flat_bounding_box[3],
             ]
 
+            print(f"grain cropped image shape {grain_cropped_image.shape}")
+
             grain_cropped_tensor = grain_tensor_full_mask[
                 square_flat_bounding_box[0] : square_flat_bounding_box[2],
                 square_flat_bounding_box[1] : square_flat_bounding_box[3],
                 :,
             ]
 
+            print(f"grain cropped tensor shape {grain_cropped_tensor.shape}")
+
             # Update background class to reflect the removal of any non-connected grains
             grain_cropped_tensor = Grains.update_background_class(grain_mask_tensor=grain_cropped_tensor)
 
+            assert grain_cropped_image.shape[0] == grain_cropped_image.shape[1]
+            print(f"grain cropped tensor shape {grain_cropped_tensor.shape}")
+            assert grain_cropped_tensor.shape[0] == grain_cropped_tensor.shape[1]
+            # Check that the bounding box is square
+            bounding_box_shape = (
+                square_flat_bounding_box[2] - square_flat_bounding_box[0],
+                square_flat_bounding_box[3] - square_flat_bounding_box[1],
+            )
+            print(f"bounding box shape {bounding_box_shape}")
+            assert bounding_box_shape[0] == bounding_box_shape[1]
+            # Check bounding box shape is same as image shape and first two dimensions of tensor
+            assert bounding_box_shape == grain_cropped_image.shape
+            assert bounding_box_shape == (grain_cropped_tensor.shape[0], grain_cropped_tensor.shape[1])
+
             graincrops[grain_number] = GrainCrop(
                 image=grain_cropped_image,
                 mask=grain_cropped_tensor,
                 padding=padding,
-                bbox=flat_bounding_box,
+                bbox=square_flat_bounding_box,
             )
 
         return graincrops

From fc938fbb06f57e36089ecb669ca09a5a190e0d95 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Wed, 18 Dec 2024 16:09:06 +0000
Subject: [PATCH 30/34] Fix: padding wrongly subtracted in nodestats and
 ordered_tracing

---
 topostats/tracing/nodestats.py       | 6 ++++--
 topostats/tracing/ordered_tracing.py | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/topostats/tracing/nodestats.py b/topostats/tracing/nodestats.py
index 4c09f75975..7f6f0e392f 100644
--- a/topostats/tracing/nodestats.py
+++ b/topostats/tracing/nodestats.py
@@ -1578,7 +1578,9 @@ def average_height_trace(  # noqa: C901
             trace = order_branch(trace_img, branch_coords[0])
             height_trace = img[trace[:, 0], trace[:, 1]]
             dist = nodeStats.coord_dist_rad(trace, centre)  # self.coord_dist(trace)
-            dist, height_trace = nodeStats.average_uniques(dist, height_trace)  # needs to be paired with coord_dist_rad
+            dist, height_trace = nodeStats.average_uniques(
+                dist, height_trace
+            )  # needs to be paired with coord_dist_rad
             heights.append(height_trace)
             distances.append(
                 dist - dist_zero_point  # - 0
@@ -1924,7 +1926,7 @@ def nodestats_image(
             for image_name, full_image in all_images.items():
                 crop = nodestats_images[image_name]
                 bbox = disordered_tracing_grain_data["bbox"]
-                full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
+                full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop
 
         except Exception as e:  # pylint: disable=broad-exception-caught
             LOGGER.error(
diff --git a/topostats/tracing/ordered_tracing.py b/topostats/tracing/ordered_tracing.py
index 8e38cc2a1b..d746351117 100644
--- a/topostats/tracing/ordered_tracing.py
+++ b/topostats/tracing/ordered_tracing.py
@@ -957,7 +957,7 @@ def ordered_tracing_image(
             for image_name, full_image in ordered_trace_full_images.items():
                 crop = images[image_name]
                 bbox = disordered_trace_data["bbox"]
-                full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop[pad_width:-pad_width, pad_width:-pad_width]
+                full_image[bbox[0] : bbox[2], bbox[1] : bbox[3]] += crop
 
         except Exception as e:  # pylint: disable=broad-exception-caught
             LOGGER.error(

From ba08b36c1aa812cce7a589e8eb2aaa5f37b17a85 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 18 Dec 2024 16:09:59 +0000
Subject: [PATCH 31/34] [pre-commit.ci] Fixing issues with pre-commit

---
 topostats/grains.py            | 9 ++-------
 topostats/tracing/nodestats.py | 4 +---
 2 files changed, 3 insertions(+), 10 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index ad6d917535..0816e0f024 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -67,7 +67,6 @@ def __post_init__(self):
         """
         Validate the data makes sense.
         """
-
         # If image is not square
         if self.image.shape[0] != self.image.shape[1]:
             raise ValueError(f"Image is not square: {self.image.shape}")
@@ -963,9 +962,7 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:]
-                for vetting_criteria in class_size_thresholds
-                if vetting_criteria[0] == class_index
+                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1692,9 +1689,7 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(
-                flat_region.bbox
-            )  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(
diff --git a/topostats/tracing/nodestats.py b/topostats/tracing/nodestats.py
index 7f6f0e392f..a587a35e99 100644
--- a/topostats/tracing/nodestats.py
+++ b/topostats/tracing/nodestats.py
@@ -1578,9 +1578,7 @@ def average_height_trace(  # noqa: C901
             trace = order_branch(trace_img, branch_coords[0])
             height_trace = img[trace[:, 0], trace[:, 1]]
             dist = nodeStats.coord_dist_rad(trace, centre)  # self.coord_dist(trace)
-            dist, height_trace = nodeStats.average_uniques(
-                dist, height_trace
-            )  # needs to be paired with coord_dist_rad
+            dist, height_trace = nodeStats.average_uniques(dist, height_trace)  # needs to be paired with coord_dist_rad
             heights.append(height_trace)
             distances.append(
                 dist - dist_zero_point  # - 0

From 8da84e29e7d88894ea2a0e02cf30b57bdcb36047 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 20 Dec 2024 11:49:48 +0000
Subject: [PATCH 32/34] Add: image grain crops to topostats object and fix
 process scan both regtest

---
 ...test_processing.test_process_scan_both.out |  10 +++++-----
 ...cess_scan_topostats_file_regtest.topostats | Bin 335700 -> 274920 bytes
 topostats/processing.py                       |   2 ++
 3 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/tests/_regtest_outputs/test_processing.test_process_scan_both.out b/tests/_regtest_outputs/test_processing.test_process_scan_both.out
index 79aff8a940..7f5d4b2ca5 100644
--- a/tests/_regtest_outputs/test_processing.test_process_scan_both.out
+++ b/tests/_regtest_outputs/test_processing.test_process_scan_both.out
@@ -1,8 +1,8 @@
                   image_size_x_m  image_size_y_m  image_area_m2  image_size_x_px  image_size_y_px  image_area_px2  grains_number_above  grains_per_m2_above  grains_number_below  grains_per_m2_below  rms_roughness
 image                                                                                                                                                                                                               
 minicircle_small      1.2646e-07      1.2646e-07     1.5993e-14               64               64            4096                    3           1.8758e+14                    1           6.2526e+13     6.8208e-10
-    centre_x   centre_y  grain_number  radius_min  radius_max  radius_mean  radius_median  height_min  height_max  height_median  height_mean      volume       area  area_cartesian_bbox  smallest_bounding_width  smallest_bounding_length  smallest_bounding_area  aspect_ratio threshold  max_feret  min_feret             image  grain_endpoints  grain_junctions  total_branch_lengths  grain_width_mean  total_contour_length  average_end_to_end_distance
-0 7.5100e-08 4.7559e-08             0  3.9431e-09  2.5631e-08   1.6016e-08     1.6680e-08  9.1991e-10  2.6422e-09     1.5338e-09   1.5341e-09  1.0543e-24 6.8721e-16           1.3198e-15               2.0539e-08                5.0379e-08              1.0347e-15    4.0769e-01     above 5.0379e-08 2.0539e-08  minicircle_small       1.0000e+00       1.0000e+00            8.4571e-08        8.2685e-09            6.5881e-08                   8.8370e-09
-1 8.0241e-08 7.8677e-08             1  6.8951e-09  2.7188e-08   1.6272e-08     1.6263e-08  9.0630e-10  2.4586e-09     1.6144e-09   1.6264e-09  1.0352e-24 6.3645e-16           1.5931e-15               2.0174e-08                5.1212e-08              1.0332e-15    3.9394e-01     above 5.1262e-08 2.0174e-08  minicircle_small       0.0000e+00       0.0000e+00            7.3054e-08        7.6154e-09            5.8272e-08                   0.0000e+00
-2 4.0012e-08 7.5644e-08             2  9.9461e-09  2.3654e-08   1.7561e-08     1.8364e-08  9.0641e-10  2.1066e-09     1.5939e-09   1.5493e-09  1.1192e-24 7.2236e-16           1.5462e-15               3.3592e-08                4.1496e-08              1.3940e-15    8.0952e-01     above 4.4405e-08 3.2528e-08  minicircle_small       0.0000e+00       0.0000e+00            1.0447e-07        7.8033e-09            8.7183e-08                   0.0000e+00
-3 3.2366e-08 1.4036e-08             0  7.7690e-10  1.2272e-08   6.4301e-09     6.4170e-09 -3.7937e-10 -2.1207e-10    -2.4477e-10  -2.6816e-10 -3.0364e-26 1.1323e-16           3.0066e-16               7.0841e-09                2.1505e-08              1.5234e-16    3.2941e-01     below 2.2092e-08 7.0841e-09  minicircle_small              NaN              NaN                   NaN               NaN                   NaN                          NaN
+   class_number  subgrain_number  grain_number   centre_x   centre_y  radius_min  radius_max  radius_mean  radius_median  height_min  height_max  height_median  height_mean      volume       area  area_cartesian_bbox  smallest_bounding_width  smallest_bounding_length  smallest_bounding_area  aspect_ratio threshold  max_feret  min_feret             image  grain_endpoints  grain_junctions  total_branch_lengths  grain_width_mean  total_contour_length  average_end_to_end_distance
+0             1                0             0 3.7555e-08 8.9055e-08  3.9431e-09  2.5631e-08   1.6016e-08     1.6680e-08  9.1991e-10  2.6422e-09     1.5338e-09   1.5341e-09  1.0543e-24 6.8721e-16           1.3198e-15               2.0539e-08                5.0379e-08              1.0347e-15    4.0769e-01     above 5.0379e-08 2.0539e-08  minicircle_small       1.0000e+00       1.0000e+00            8.5729e-08        9.5493e-09            6.7087e-08                   9.4020e-09
+1             1                0             1 7.4313e-08 8.8557e-08  6.8951e-09  2.7188e-08   1.6272e-08     1.6263e-08  9.0630e-10  2.4586e-09     1.6144e-09   1.6264e-09  1.0352e-24 6.3645e-16           1.5931e-15               2.0174e-08                5.1212e-08              1.0332e-15    3.9394e-01     above 5.1262e-08 2.0174e-08  minicircle_small       0.0000e+00       0.0000e+00            7.3054e-08        7.6878e-09            5.7989e-08                   0.0000e+00
+2             1                0             2 7.9532e-08 4.0076e-08  9.9461e-09  2.3654e-08   1.7561e-08     1.8364e-08  9.0641e-10  2.1066e-09     1.5939e-09   1.5493e-09  1.1192e-24 7.2236e-16           1.5462e-15               3.3592e-08                4.1496e-08              1.3940e-15    8.0952e-01     above 4.4405e-08 3.2528e-08  minicircle_small       0.0000e+00       0.0000e+00            1.0331e-07        7.8754e-09            8.5725e-08                   0.0000e+00
+3             1                0             0 2.0510e-08 2.9845e-08  7.7690e-10  1.2272e-08   6.4301e-09     6.4170e-09 -3.7937e-10 -2.1207e-10    -2.4477e-10  -2.6816e-10 -3.0364e-26 1.1323e-16           3.0066e-16               7.0841e-09                2.1505e-08              1.5234e-16    3.2941e-01     below 2.2092e-08 7.0841e-09  minicircle_small              NaN              NaN                   NaN               NaN                   NaN                          NaN
diff --git a/tests/resources/process_scan_topostats_file_regtest.topostats b/tests/resources/process_scan_topostats_file_regtest.topostats
index 43a4380fb4e37d4f80ff3960db834945796860b5..9f06e06d6c9f6bd869ebbabc164aff37850b6ef0 100644
GIT binary patch
delta 19151
zcmeHvX?Rpc)^Oh?q|*t|f$U*vAV32FI;;UA_rkD5)DRS**+L_t86^;9izCo+fyn^I
zbgoj-5e%{zV1xvs35z7bFb3H!phHv!TSW~h7)A#XQNB90bl(oiyfZ((?|q;5`R<c@
zPo1h$r%s)!I_KP~+xhTNM8lWX0=v|ZBzyY25K{MrWP+=*tmY>nR3~YZDf4p)O-vCw
zh9<{~B+D5u(hNPGW;Avb(f6;*3->-NrO6h38LgNz+*8vdHi&WQ%wQJ=M|f&J)(R{^
zTA?KppeP6!L_jiL;1{E%fMjaLVNu$nA+6nj8S+V?kDL{fOydI+up$DDMLf+&IBJBz
zp5w7N0XWMsBykU7pclAf06Tc}zoX)Bq#Td#CfO@%BQ@oS8CCTXTBA5!3fHVhtjMNI
zt+oE4ZBbQkAy{iP2WgHFs}X72M&ON;Mzu*HQQEZ7Fu|mGEb(Rf0Kpwgr7+nT9LV5s
z2ARO7jtn<5IG8~jgH{Y(BPcbm14ZLzcN~0e=nkdhcaL5hF@K{{^zFA(PF~xm1pg9x
zBhh`(tt`0}oWAJx9>wYOK3(y*PnCBr9qC=1U8Z=xdG($CkCiJ8!>4`HV@i#pb|_nM
zNvTnKTV7fG$n|4Nuk$^&?@m9eJXw3SvGLExl{d!EJl^|*GfLNqE8`D%>XgR^{Tf&N
zbft3inl?A~w>`@Jo_RfQ#+NBaJwLyD_nwtX!OE~d6e{jb3hRIm{Pxx#<mWmc^5y?!
z|JFw=&f`^wSG+oI`{pxM<IhCfgk7hrYOSZEit>+D-9PjCLrDoARF%5JQlkf*s2bKO
zwa13?(^XH*_|t@?%RZ`d#=rks`?Ke&BBmca^~=j=t2#%#75}Hp$EsRGi!#1LQTR+D
zy5rN;PyV4*zlYT!>+9oI9UZ5J;3BS3Q49aDy=JxA#gaiNsgKY}v&Nb=m(@C;ZV1(O
zV>=2{fh)fBd4XebgW)lmP=A|N8qsQRT6m!(<}5+ad@Dqjo%&qR8k4M2(GU|`xnEFD
z$y-^5p((U*vDi^_MkgB<uJ&qlCrh%glBedRRuJP7GY5KVh`3_TAWzMq!d15L4wSLU
z2eqVCcUv0bXm)Ol)_CisitNkM^Ade}^aa|2CAyZsJXOnITA-bLs$RR$#f$IHwPoEO
z)jB0SrB!tsipy-;#-~#~dKV2PShd#^#(Q+LwsPqlhMFPE$QaT$#g;|N9AOY?9j1%*
ziQKdcsP9mDQml<F4=c5aj0l@UPfb?K%``l~o7vlVD_;;Dp<3>5=bKT}6Mq_INe-pN
zRta?)EJUFfzZHT#r<OI({h2kRESdA0t)S^qF;dZR6Z&40BeiwE78@1g(bsDmHM4eT
zM7l@MA<3%+$qW`?=)`v}zBBO+<p7@|gRxP2LVI(BqFpV(gMnFGE#T@WqXAiIL#O@g
z5kcqkL@zqO1|I9+kq?gocod3Wt@Ht}Ry^6O6+GzGa;JE;?5SSO<@9Rl)4ZBG-K#lf
zcs1KhuV$U)EzIz0jrjJ?_G-Ed*WUw7%Hi=2JSyR_A07wcaY*!fYQF!!X>vD!#|qL`
z))KEZkh%unJRu40E(i=LA~~qImk`+rlB~>-gd9LM{YaD*5=I(rL2a}Qd0O_7C~aP7
z`)+{vu7H-f+wl2U>S~GCB~F*St0>Z)LJNf}Yok!VKEe*joCur)nEILuC@|(43RpH_
zK+@W#))B<ZQ=X9L=Z+DrIxJG6Fv4yqH6&@TSS`pGv^gRFn&~ccR1l>sGe>lU6t2#I
zzj><kG7^f&TA3vh^bAgxj57A>?U|Y=%soOOodvS?vLzZh%tAzAYq|M#>{*M~@p9Up
z8Wyndv5{@19hz$q^9F;VJ$}k8rVb@0_hGZ>Oomhk-7IDfgOm&BWiIyAtS_ARkd)fS
zgg!qlmuWXo$F*yGigI_bPI*6LZsB7*1!AmX^-`!`9b*=gpN7=@56xoT21sE~B-U@F
zrd?b9ky*la9QFK8m}n`<qs(lP+@iZAt@K1DzHO*@fbg)U&rioljHYJUHE-4!t?_Xi
zs_ri&Yp(1REqy@;LRp<CI5hRB4UocxD)Ksn4qEoYdr@OwAwerx;6SciK}My0g+w6H
zyxEg6T#Rqb_~gt$yS5Q~k~%HJmW=&E5r*mb&h<g?EO3+joC54oRGcoxFj+gfym;@*
z70(7&6t@=|1UycGovqty=q&gOC=O>n&E$jBPvk75#9L@A%yjPa6M>u?_8U2Cj=cMY
zOwIk<S9dbk1jjIonZ@8+#DICyqC{GF(lE<k&TxVL%}Sp?G|ZCz5v^2U38(~~MKrZR
zX(>V{OG7eUrJ3z1UUDzR>(VN495PUV(di`-PeGhbHcnDUlQt9_O*+28(QwEe4~LMS
zhi`B+DuxB%@+k^zFlr?zW^SAaFo+gWt9R+>;vNw{%0*o(*`sTDdvz_dT-Q?H);0S(
zx@NA>wFZ3GRpRn}xO_h@KY+^*;_`QK`5|16@46~nUX9CZaQR`(@H}w601q!bz7Ta!
z&1#?RuMUgTN^8jtVSFf$=Sc#W#N%+k1BCqAx!Q2;1q%*iYNN{ZSNz`NKj)B3IEHBM
z&<Ie)Ya(7Dcp=SSyncY-)V%MfqJsVctTG&cJq{}_5#K5JhIBF%s>s?;NI}{Cg><l2
zR?YhXUCelijsZ=>@^c-BHH~35gY#iI+EDsHVJ2>6GOua(=wmSBF)$QfFf!DPHM|Z+
z1<cruVq#vf$uhqIDKKL(n;`{e%)W)0v224GYwx^BrzAIHE<D4^wn$f~5wm{|wZVuL
ze+4NpVx?b03XGUH6hdQ&P|?LDq7{9&L^Q)=w8i@k<!Wk&4H>NMI|zO@5WECEx+-8p
zQ_4ctY;062u_19`Sj~)0Hh~wq0Pez}EI`Oh#dj{gRX@zZW%wRRcmQ1`7$Xd0Fs$CY
zYelxq=eS3Qg#n8JW<S7e3?a|OZb(D41`7wu00PcwcpA!5@ePX)78ZwjI1Rixti)-y
z=<NQajV7%a#Ly(<fILNX>D{433mP(1PL1vGB-Nzayk)#Y@$Fd<(2j=AGZ#bX97tu+
zEDndo+nrLXcbLVd)Nx4;i_QM;QeS?^VR7|lQiL=$gk~+bl~uI1N55SWGO%F;+>DlA
z70hBD&J=|Vlu#*FqlEqNI0%nJ@Ti8zVR+QS<FD|jgU1K(I1Z1K@X+CL1|A>6<8Sb&
zhsP)II1i5t@bJRp3wT_F$Jb)Dr)D^_M&Z7#1?_{{Sa7EJ)v~or5{0d;c!$BFB_3&m
zRJ&kCpY;;j;BAI5Y_P&OX%Hp>!^KFldJAo}PGj3)CT7hx){0>W4$jPItsNIsg@up=
z+GQ7F$U{q8U_rJ%LR--!Xr9m*EDkYJ9(~dOcnRMH)3$yfj#!Avn@^#2K?=q@p!EI%
zwtZRm68nkcnnBpGCPIKS2vzqJ?gdwctt+vtnvUa`bO$EM{U9B3?6@ce5(%hyAU3)<
zWy4vMQ9+UrgIff{90MaDqv9kXh4P9(jY$H20hXCq*sN#Y3!>nL$6yGS9i+j=+Uc{n
zjbU_40c>t*8=AIa3tQcmrory&R+<LO>kWtWdK`#!MZay3hl$t#mzBy#(A6P=Q!2*T
zsAb}2*kRDp!s8^@J%aYkZ;`O<8w6jjI#ED#9uRo8LyoaD&Y-Kvnn6PauZyj!xfM%s
zC}E9rvJi`9oBXjw1s152pLt%ofIJy8XU>=1t%-tJYDg1MXSba7SZ)(;{>X%h!A!cj
z$*z`POiCBS=p0Sx{NjVUU^Ao@+J6a09hHNrF2N7!lR)$~Lg73cig-YnAf#8!en4m|
zh>Z`R?8(9w_s_4@jTpS=v?7(CO0W;zt9(_I5iz9vjB@4bq=zr8*`thdo%*Tk<9n6o
zJCuJ>eDI@`!J%Teg})q9u2ud%<oRa~D%~nKFFrY`Mk#9b^7UcQ9#F0vIeYEIFQ=86
z{|GuiNAM~;f^Q%E$EU}X0rB&uv~GP(nUim=%Rf`A)JASzr&UzAmCqiOBR9ksD~DYd
z)q|~%E9ZLs@Kbo`R^`2?9h(MzwM+T-`hLsvYfF{={YMQci+n@5XrTi9^h(P4`(kU9
zYfnr%{o}F<<<ng!M-2Px3B~e4?03&S^qI2u!`6?BT@ESFb@NVJ@ZG!0nx2+Z1w$&7
zeV2=(POYm}zAC(L-``ibtCZg3>MHGb?N`KsIoo=@`l&Kw@3;-)fA^lUBK+fd7h4@w
zwAX{D3LjrmIwNJ{#WyPz>4my`LvPh6f1f>I#MbdUm0bshzI>?N4y8?T)sEK6g_O|Z
zrLRufcP2$Bed65<@4TtZa}9aA;@mc+HOn8z@+V9kw8840R-W?Iw4rMs+E}VMqSv?^
zgAXb{oZVgD(D{J!`eoPrHLH&(M{a#HepAT_rLFqPP}AWv%GKR(WM3MiE47Iq_nbBP
zuu{l&E~h-RTVJyqDYMdBO*>QcR!T0@r1nW&O>p}0ls_|gdW$|~)34=9;>?28?%h+i
zE5rVN-{ehi9#F2g`{exlg-4apSFd*|w0)_ZIQ&ZOl+dHf%tsfkn3=Of8R>OQ+S7BZ
z^3u`r4PWmpQ+nl8gj}4pP3ac(P2SQu+m$2BW}K-zdr(=M)%ZhFZIyB@v2g10rSB-;
zroE`7zW0qX{GWrK4u0yO60|hBZJF&KN{=6s+DuFO=9aSP*uER>=9{qV+&WcIgzM<5
zsX|X-3u^5YdWJu~X-wt9UDuRx$%FTnZM>nZL}Q&oXW>IM4};kJ13zoQs}3v9FZUkZ
zcCu9Yq9*NX`{75GIJChjL<_6Y+fHGsXu21*nI=RD6H)hR!XJf3^wu;1Z+<_;-kSf{
zZ}#8S;)w?{F0C-+o58TRXxto@QTQyOo0OMfw6oOaSw#mnrn7|wqW69<Gy?_g*N3}>
z^krQ8YT<FCLT*f^hAb1kxLw}hv(Pzno8xDjOZn-u{rphkcftgVYYJtfn^3`ac_TV>
zMTiRXZ#xSnWSu1>Na>HaXq(04Oq6j&SQ=KDNm<oD#%nnj2Dl4u(IJqMZagNS{!fdM
z=(Vdts?W@hE?`V3D||p@U<&<*6AQT1KT!%!m&x$FwP2>NYuGh~RnMoT-UVDB*P>>$
z6KPt=xlrCo0MMM-!Xzy%+VWTtW5bDzY=M%CWLickIGvs0sZY@y?s9Cofpn=1&wLg{
zKS@Nd%ce99Pjb;q7B?fhWL`wHm=#ydD)7er3}D2umxVXA1???j^}~iKfZv!Vh>oeq
zu~ghfhu$bO9s1nm%`_aD$vOoJfQ#zvna%Rtum)CSUe0B&ThR;)$-bP)Ud<|S+5HRx
zWDgaV7(~A{M3^8pPDSc6F1lZjyxrne&9obt_bhMd4$;le(IWdQ&Z1#;Gc2UkRZMg*
z<Hb4p83c%K5v^{6>~fAcLC8cCbHuG=wEdHpH(Q<GT(c!>Nwc#^CHp$sXlg#MR9Dzc
zrN~TNZoGW6dq$jzvET%LmP~qBMUDbmoL<P=vA)zyJ5qKb(}|f?;9UI-n&`xl$93|2
zH}vCjak`YbUZ|L~LVP*b^H5Z%?JIg7QxGMh{<DRSQYlW{z9PVx3?uI(qB*~bdz(3Q
zfwQU1jr4VtoK3CsnaGBBL_=uHOs<rJkGe7QSw6DH{bGyZs1-NdGf0~C=EWc>_rLSp
zqt63r6UX4=ViP~%ZDeoYT=|BB`?8vj0-iasp_P;-XTQj-`;9|bkN&Vez)lxoHMBh^
z&m~dtXa}py**guqV3KiLm^zb=9jD7?ctM7?DW<bWKS!pp>mVDTLtGNmdo3%y7)v>H
zgOv4HGONa!8`&6zFgH;4oUY*p#m&D^=rhQ-0_mK0Q{i*&p~CK_i?o_j<sIh2rq+rT
zRUPF>Q7SnsB4v}5^*faj@g_43m5idPp9eOX$?&`{MNiE?iLu76N3zFYv?Dq{P&P{i
zM@Twlvt*Wv|BBOJ;dSs|9c>mnsU-jY(dMYtqb#NGQ|WSJ{5vBR3olEnd`G)qYw9MQ
zY!le-Q8m}N0?joxoqhaZ2AByCFUuiFHf)(%HMJz2ZZC{=YK>c6<kq+41#I);ea`w8
z4tU{A#lgn4>6r>xD}y50<gWf&Y-YKad?}*c`=ljd_KUQysb8~UHD~`wbF@(pcZT<0
ze3rsjM2F|_KxBR`qJ;;gS!mrh*&1(R_4$Y!_^|~S6Ka@v6>~02n&i>@k@tGn2`AY@
z@#rD=;Ka<*gF#b?NIfm3VznlsbX`h;#~_-rUY2BJJ&jeIds@Ofzgbp8-Hw_q+W|My
z&q%U1=wv3<Y?(fW#sDy_9A218u#qZN{*xr}PO4OQolUA{%qgfnn=O+jYR9v%g$ku+
zx{mIw5O`G57(waCiId+>I})X42G0JAh@v~nOG2%^X{{2R(Rq{9QL4iUya>TN##($!
zn10b@X}pJ0Thm}oH=S+#n}ss{tN<G=n(JsH8;l2G4$p^2Hk-{4gKEtv_mUKY3NGQX
zuw9l0dGs)%anx{0!h7U+fG$g6P-6FI^O49E{@)&9pznqXHYumC1lH<TUnm%#%j_VZ
zx(1DssS6<^OqcqfV2%(gBv$vs`s;+i3W<7uY9HUnQ{Z2nn^^}?-T#YovldRp$_Me9
zenUh*n+V+g9c^(ArIci||0R$9obS*EA8P+EJg~382lj38P_nNQ{_}Gsax&FkK>tk}
z*>q0krQNy2!-ad$>f6#vNuNk2EOxJsGzoPU<*c~WCZ0NQReSW?+;fUv2brhfImMVs
z{yA%(?w_-PbBG;f*zczkNu^EzoHn`e$3Xd@?*IRR#DDUEghwC8=IH<2^NR6!=jz-<
zc~tQHcNRV}VtEhLIb9wW+<V8ld&hn~2+dBHV_R+Uycs)fKU$eC4-c-|IwE6C^lem?
zE|1b2Q_RRbNlxzBDKC4YHtMF5(dyiVz0RA;u_>|pCv3T?j6YpIcl3)lmC?eEz0)IZ
zqDLpmN$xN_qu_%uep_%F@?d`mJ}?B>2Hz0Ww`C9_Yrw-9Y{y`G1|u+pdXWr9VF-Mq
z@eS=junTQq=q8Oo1~dW$0HFdj1`GfKLx9i#ltY66K>&at5I~>-$ZzNXq@l6qH^}e;
z!`}o);-k<I)b;^6#P-9<*I(&b60TAh$zZJd^#*0&#PrUpW}jt44KoE>7vqi^FgNa~
zL7AYPj_uiGU6X^v%*lLN2Y3>49Fnn9P8A<stewW+L2cYF7toy$|0|2cBFVZrcUmi}
zRGJ~7*~4Y}Qtp^dn@OL;E9vu^sU&94LA&-z8GKhJ&dlm@SM29Q&DCx43X(jRa(B+-
z++!Y<kg(oxW?~IZ$7+BpeF6KJ+q96@vOU4S?q0&!ER)dd_{BQT%J}lu{q6A9;bWrY
z8DK`qLA)!IB;h1J`}dISnoJut6+BLPb5@lMug{{a{s=GcDH-zpLN$6ZL;gw9+nOpo
z9fMj6Qhi$!>S7K$7Ur+RyIF)qZ|Mc~MwUc25-VPtE5`_u_sM4TEvcBnsYmenQ4Ure
zTwFeba!5XZfh}%Am+aRWp@Y=}m;Dcc{L&eoy`1)UHy>g#_c6#R)1&cAI6i|*n=L0x
zUfg9*qrBI*>7XxbPWl`?P+0o8WCPE>_)b0zTy_5Dx#}90H^;GHxq)NMQ=B9Iz)UQi
zv3{{!8<y#&9m{)$bE<y6IZh?dNOa!>-oEV5(n_h%^GbDjEstF0^Ac<=WW||SH9jbQ
zHdKjp_B&Qn=IMEyxxP+naf>D!H%v|p&8!+{Ze#-^$Vc3{CvYTl8D3mY^ZAH-j|9@0
z|4MOt9Tk@D79_8sxz2T5*qD4NaH-O-3dk#E4(<s4qL*#5SXo^zv5}Wi;`&zuB{BV(
z;<gP`l5;Tk&?j%o7SXX1XR)>WB1w|9TL4KC%=a!R3}jZx@cNf%K40LKMU?q6IQ$&6
zTAm=4Opz*{SuNugTrgFt@ca>{#Cc7vsgWH!Qh1q%cQ;;=(C%GwmYC}y@{pjdXuVH1
zu>{dJ-zO{Q>R=%m<!8uvb7G;q(U&FZGXzvqCZ9#?UX<^ba-0(S{6+a!$u$k@eX;yr
zm|s17NQySJJR~YE!$XoK*1ZCP5VPmqNRpts#@Bt~^|rndsJdi^JGTU?F1dtu15_8n
zH|)%Cnj__6A>CJmNbo<qnv_)=%8F6iPT46w(p0_1CaHS$EnK}hTbO!Hn~2ObCf??g
zU$nnT<?^dn^3K3tV3op{nx;!ctKnQ=9v7R@UEd0$q}1(Ha`?Q9f7XtbB|iv#hqVPi
zmH-Iu8X&l7h#7(z<UtM(dLU?N@%Jr;uwhv6k0fm9UNP|j{NZk^owXA~)_jPjOqTbw
zhPL^^TQHjQw%lF_LO}=Q&M54FYz~2@;R-^1-<G3M!rQWS5;O}JhzA$@VBd#BI0F~2
zOMVD~ajq7QLpR`RgYMcV_eZal%VzxRpAGwPJ1^~%t!Ujo*-Rgc0o_tATiv`Z2%F&I
zvNYw%ao%QAnK<6G4?GOxd|^o7VzQ!7_sd<-y?5acg<jhWBJRZ^&VO4rw{Pb6L$H=2
zJ$fr1z4byLzAb0FX}FxMenUdG;VO-oM=+S8&RNygG&f-+gNj;uI5;<{ol0S$s^%XY
zv1rKnQ53#EN~N**XbNwQRwtBho~D-lPPIc={>T`$)s|3oLi|_?|2%dm%NVC#1OfDM
z>S45~$K+Mx)nW+Uho`D5*M1t7@jZhJ?^b6Lfsj}s?`iwi^+Q|AcmJx)%v$qYMAUES
zScSaZZ450z;l?=yWy8wrmLA%0Tj9#R<7aQ6VSGjYRAhN!==_m8x0we!ZVnjvU(uUK
z))kuvP83FJHt8WL4=*oby^QXNuJHN_shRz;>T(G!+$=8%Q{Sf6$5p}2oo!#s<}SAa
zM^W3bFOFf=xbyL|p+<tdVpFG*S_U?RSjNo#Bng%==Kv(BQz?3v%D|abGu&H6^SM*0
zevdLIP9<eKI2GLQCngJ<!OM`0(j7FzAF)5$)Cd~6Qx!R@xe+u~Gb5OHkcdZ&U~!F4
zRI&6h6~SNxt<xoD1S@w+71*_O7o@tK5=sgVibRQn@p9om<H6tM*T2OE${!DKy}OPz
zH&Cf{Y@q7us4^Np;1x=aw>(hhV{D*`Sv9U(KO5Q+!$9e^&`&Z@jVDMF2Fi2_l60VK
zrvnEnmEoCZ0|zSg9A!>iSYJB9P3(h&w<sz3AmQC6Un{tF8wz(*CuW<Diu2@yVaB%T
zuzAn&Vaq$qhRv!IVac%BKk^AKCjX5JWEi%(Qa)_;rObD6NxARpi|)&iM_HWpR4P62
z>r%46>*mIi*WpH*<J{PCoEvA3b7RhN@Og02CjD_I$KkyiE@qs2lNWBt;06O;z~SQ7
zo#VmiuKn1vO)i)1XeRg+=5t=aJ|~>F>Vx2lIty*kyBOTf`vC9yaPdBJoa+Gg0N~=z
zjW^BlaCBd}9Q}Ka^2WHI;;P^{R~7fATnvtLF?dh87+ez}sP6&r+XrwAC?Aj=)OQ`+
zE62FGws{pQw_T;M2ZhyLDOCHY7nAe9OwR8@aZ^8a7@%L;817Ja?tqKqNi3hy=MHE>
zKhWKT;UKlc`%mD%u8eRHEAPe_4`%rZDrq-)2!lNs-4OLj!k{OkAEFX=9eVl7DgWNA
zo_mP;Ai4VZ8ZeP7(I)_{kjQ#6h|zI;p!(`7a4|eLsR{02^*BA(EX#MB%#45&wrf%$
zK~0!F_E7%Cw2>_ac_w-dBv-V0;y7Gg8H`~9C9sahGEH<R(4Cwg&kA}nm2_e~x3iww
zSUy)-SH_?Z6EuMhSt3(?cUC{q$Cu8LM8=QvCG>7Q*u8u`rmF1C@ZN0f_@wLQ8!Fmg
zyRV$$JyfazUcZ}<4=wNNn+p_&N$Q3qACnGD!XB(gF|0y&)^H^2IiHN}SdY7~9$T3J
zoiSs$tm-4=>Yx%7;}piRp|Pr)-D~^pZ+Ih`6-21(09_m9gQv$v&suysQe8|Eu?$A4
zK-F)52bGuiCmA;N5@a_-vHUJ-^ttCg`MJx&k+fmc9Y1i78q7*;4DMW%(dM>TAzqf_
z1XQX^aj4&H*(o((&j1@Q_ycVGIlio-#X8>ii&`w_!8>zUu*u2lJ6}fd{gj43t(Ib<
z$yq7>vKoz9Bw9xs{!0R(&SiB}wspcivvDg0ND3Go{_!hzzynJdni)ln7q4dc{)JFf
zGP-()?5P<GzXfc>y}-W(JSLg&p|byLZE_C&4fQK|mv~TNzS{CmOvxp-6Qi<fT!u^^
z8|r51PE7W<lzG7w&fNNsmUm)Guh5;C$5UA~&fLfb`UY=@a3zOOAIRxk49^Xv`7xBt
zR1h0TSIO}Dc&wlkWXCk;dn#$sH7=?12NI_>{mT3=KjM@9%4=BC&*ZGp>8u);z{my~
z=aM$uqU<@H^;e2JZc}xU9T+%66|EOxtVsw~{3A)?5G=0xdSLwqhMRw;^>@Q=Y2p|1
zbc_8{O6kJC_rJByzk_Hkljx?Yzbfx(?$6B*IO(Y!m%Vru?{1X-6NI5>yLC5dtKW0>
z!^ASYq0=Zw3X%=CL*7IDUF7YQ+-_N;>9`DQbS}|B{V$~B_zQZe(S(PHe1Q;el>$dG
znc>bL7)4OP^6OXudFPpGO@RtfxdIIJ*nYhpBsaHf3<?l)f>8cdd4XgSK@j?dO-!&{
zR@2(BJzII6cfa~P>%MU=(qOf+jcizZ>~k`)I20||EH4eSh0)%rZDb~u$tsg~X!g?z
z{%Hvr7FQb?b&d$iioe08I-ghuHUUTTb~A_%(@IJWCrLc5;xpR^s<eRNrIDC~sPw~D
zLDMBK?i7BPxgs@QV3Wu|$fUJLwX8!ulIqae4qB!<%#LoN!y6a44o#mi9U9gtYsK1G
zZbpYU;;9aEI?(PEb>cd7+T<2>SlNl`Fr8K7IyA5e&|z*2h);A_+=V1Tht+)ob=bgg
zb3aT1It;(_rQ)tH<hw1UU1>#Kgo)~=(Tkb6&VI~QhB4UR=Q^otbubqk&Y*+2;s}-6
zyrd?soCl=bC3D=y>$_RSuB;*t%Mw^Dz;PaM_Gk1ws_DmY9)xh5M?-yCc|Tv&WFnV?
z#TM>9tOpK8(3=$`HHlsZFo}8CHh|IhVBy*T7PNF{4Gv&1!N-q=WjxYmnsmS4g%NU9
z@>rhtgyV_6c#(9##CyzDYi9z+GJafRJlNqOd=xA1N(6R4_+TW<u=&C++CV3opUcNK
zj^UBMF{hz?tZ$Ccc(S8!4$(PdW4M`xYh77{90dFx#CHNKkQiZiHW@`m*i+q1!Zu08
zu^A4~di&^U0Lmj{GYi-Hs`M+iG9CY43^>cmaN*yR;nht}?+gott?t}@E!z-#KWalt
bd&A13Hl(h<--cZHJz>p`2pzZ)2;2S-UWv94

delta 32072
zcmeG_4O~=J_VeByWCjKp6a;jXQAtoyM@2=Y-e^{4lvwz+qTxqI3Wi1IZWg$)t-GQz
zKJ}WbE@qh(Dx@o>q@uQ3Vva857Zw(MF1l2d=3@TuxgT$4U=X$1fB*gSeh%l}bI(2Z
zoOACz_q_YwyMFPBxng%nehf;E^l|i9Y)pSruY>DQaoI)5zlx@HI*U{9nApKTWR6S`
z=@%r(*3Dztkf&s|{kccC{vO?@_cPC6b-w81DC-*G&lyB=*p9;?xYIHDx457m@=+7~
z!mzEY<cr5XfP&NlXCStZL1r~6#}}KA%O;#O1_dmw&_gCbEieV)J>4W<)q2bT8KX|p
zN%X?O15g_r8;AVWfqudG-Ds%;RGH^+KKmbPCnCba)w6Qgo5R7NZW@Tu3fl&tP`q~x
zGN@aet(h)W&QPILtEI=xcyBkU6;Y6;xw5V_5RV#*OzNm(ZE(_9<fjH!gg`>C<{bN1
zb_d?1ABX)p?8u?Sp_#zk1eOyxC_37z)dQLPUr~IvEHG_+aFcTUYjs-0k8da|#*SY7
z?&4iayFvPAQ=TbO9yz@GgR8$gtxWy$?d}Dc#meH(S8eY-b6d7jKX}6Fu4BuTgO?wi
zyy?Zam04{s2H*OdeM)AL?tIk7qsq4jF4WfkrCfQ>_(Il<?++`vPY?c`?ZKnUu(&Jj
z3ZB}bWb+2=7JYN;IkZ)Yk&d?ip<uZ(*|OxGy6tlnF24r4oE`qw$r&3tewt1f{Qm1B
zs#hQ0vmo-vhsJkz?&<YIOaU^L@A>od&+XX%=)pa=Uh3GZ{>6iPqPjjD7rEfbo+qc?
zcl*8_$M?*My!EQ_^ap!<2d)Zu@T(Jhj<lMb=-Bs<J&_@s+W+aR!+Q+8Ns~IdI!ra2
zW90JwYA3(8YLeO7w`L$SV!F)W2u@u0nWu>@pR;_zmr&K=%*R!!l67lZ@H`~%n5I_7
z=M_7rN~*JMxNPX<_)SN@x?_w}W*+2R&8UOJ#=1>5byb)LZP|Njet4>Uvb&?KqiT*Y
z$Y1nul<igvKD9;!bJ4_z(DG_<ySwGupsgwGK1S-2=ood$vK{J&(edhsi&r|F;c8&i
zNdWiVlB#z6X0O9(Q+KO>QnRCj97o5K<N}gZ)Y(xnYGCwlasG12!2Xy=Nj8Twh;p1A
zHHA=l6l??zTmQQ3F@3Yr4fQ+OmkD1Rq}KE+yoM@M1`hH&G)xfV1hcfQX=>cZS3k^n
zzo6b~8n8r?tp;`Vz{;j8sgC$TeldJd{3tckgt~6fAiu7BQWENv*3xP<XRui|^u?D0
zr0rV^md!=#dv_SrCB-Kx+<m)^!d`bza5zV)ft#uE-@7AK{jg{uDM=$xv;K!gVErIk
z;li*a0_V8wOuIt#0yGi$8c=)fg?C%wf(qG;ZAax8Y<J4}q{rPb1YP8rjlt&k<?*UD
zH%?78>7kw#+m4aSf(^|Q-v?Zsgp{R7r6i3c1`1NOf_S$ADo1JsAtp8K4GU-~{X<F{
zRr4rqP?t(&WBf?tTq;_FWcRu?1V(72-aOJ@282^LNmbxTJuaNO#VfYhD|V|_?BBd%
zOT1#=^orf(75kP~>>jV!GOyVEUa@a`#UAvEeMhQtlwE3Ovaib;<hOAqv+E^am|V-$
zs*lX_rC93FstoduUXXG=VUQ2?hEyuSowPy9I@cg?oa!iBq3)eKRgN;LdFRSo+tL^_
z8}F;tM>S`gA8wX)4?)en;Rbo$BP6xr2r^2402KI6^K-&?hVt5Cq1W>GTz_ep@iA7M
z+M4*#O5(*&5yxj&Gb{7e+Gmq-Z7(SvSM`<>acV62W9}`5tG1j3H7Cb}&AlX>nx7M|
z<}XN6ZNImxdr!rv1q9j_*i~n46cDL7rxVr8+;lbdY|4tBQaE<@mGa3zd=9CTmz^Z;
znR{w8NhYalzuf^YntE!h$Wderan>u<wM(lUPP6)q3SjHxWxw}@E(z8&;+lr4#&&o6
z=@tvk!0$FD!pHtD5A{p^tDvIZurV3swYOMMsp$nD%EL|IXB*Q3%W!NoOWll#CRQU6
z4>3!B^G@0xTY{x%!+~Zfr&P&9jp=`9R5m=ZwX}x1;rUm{^6GCI%vR@W83J6t@_Sjf
zW)v*e#L^$i+`}HiIaa+<J$lDLil9FFE|N?lPB|ZZH7gI=d;)AteYcqxucS48yxmSn
z3j`Zqb@%Odyu={w*t(?nexv4osGFMS?tWsx#ctFHXqI~F9L<?=w?^n_^-NYkJmd&1
zH6Gi~6CdsJe8kIfI)Y1TNJm(44@tQ2y5K8UIEAApwHC5kfM?4G!L4kFBkinYrIZew
zD@+8FSUhOeLk&>SsKc4+#*5aYDiAh@m@+^``Y1<M&1ud?URka6hIBsSg|cX-7i9jd
z#z>sm9(lpqE#9!^4&Ja;Ua>_Tz0u}Gdc!7Jy<x4LykTp-V)u6TMw{Qo8#Xn{8`gG<
z7wn3z$XQ%ARQknld-9!n)i0~H!|4w<eDtobLkYPN43r$sxhwigso*S$bG72sK2j=i
zxgjh1NJ*;wWHOAiI9$+Eilt!=folF7AN8%X3OH*L+JK8rJvm&hI+emQ&Rh#~>lWgX
z6R_PTB@)gl)Lok^QwY+w^lTdO*ePmKZWeIE=l>}`thsB4b2?qetbe473t6*l3n2yT
znEg+Xf^{tE&ya$3EOilI$E20&qH{;ugj8_X#egn1PEyatFQPXMu#Qzbp(d<jDept7
zsM;W>e*h_1$Ws1I7BcMYjbe;#IfHh%PquWo-|7=AUlpt-uHDwKbpDs3lq_t1Nc}Lc
zfGiv%)MvV+lC@?dy(I{Y8lf(!fpn(Y@#~EQKZ4_nv~pHfM9O#2<?6m$iqzTp8_6O?
zaP+s4u1(~IgUsYF^wd>e7{*-Y?%NB9C*4SJspL0;{L;wJPJWpjLLTr$sLe$bvb-o4
z5Xy>(`^_ORhQNhhpq1cRyTrK@VajA`w_wd~+<t9UxYbBOY^A6QrYj$=>V*sqjhGNB
zXy^@wv_rjT-KSpok^2gxC7RJD&=rre^Fm4PDU9R0T^A(@r<V0#l(h!9`zVz2WbD^Q
znv3*EQv3N%Ds~o-sj9{P;7!S>4)JizbfMv^h=-qd;rQA?e1E*gjV?nx04IEHtN>bY
zIeJ$G;@i2achM4@YXRbnNvW$b&Yf`KqL7qp8FRJWsD_I(slhpWND$+aR=XD>O<Nh7
zkY;KGF5EndBg>;%D<mrzXC+oXM#NS_!pD^@V3i0RU;lYOln@bfWL(=uDrmI6;i3?4
z*vwFGSX*0f*jlgHqA+i?nc?2Bwg_+7TCdomcHU?++k3;>EZ(rSUa>_TywPTM^oF%X
zdc)Rw#THq;(Pnn?hP8F}hOPCAE$ZTpHZ#f_)^>|jP+W%il4|WL1>(Q~$f%AqnbZ_h
zsQTzX=ys_QXOlfsvOQ`bBQOvjO-7Ia1TCR7{!5RMe=c&3OiKi*TTH%KiI%j48P-=H
zJUdzn_Ld@xp8*$Lkb__swumkJzgihZfDlcjj#8Z`6CgmKcUze!;&5#bDH+yz@K0pc
zlE45~@TChv2(^h&CsLSmI$m{VhtP7m)L-c-MbpSZ<P^1FK@8y<r_s?c!o`SqD~LMw
zo+2`28(^<nukJldgjo<yGuYgdL?}cG3$C?EBS=s&9)2XMNns&MeqV{)Vdf?~oVhfx
zNUx`hJqaRWe}ohS7714&1%XAPj%_;{Ah?iHRI6W@!<mB%%Fu7n-as9`Qik5f@Aj9*
zBPYo{oh_}W`;*4pM7DMqwmOzwlo)P$62pa4oNCjD`-e<osJI8|{gZww%SKHvxqpb%
z(cCUSkFA;vt<$;gw?M3tf_tPG$yhj=6`M=WrAfZ>ijJ}YnUW+tYaz;*xuXFG*DV^a
z#FypR23y;1l2M9qPM6$(&LuK4)DE09OESxq10V}`Yi7CbepYk|J4Z`1rJRx>_e(*N
z95MzcJ|Ml6-R{l4A9aW>P(o_%37azExN>RVM%$muHY&sZHG2Aoe||#=NIiP7^Q^6k
zH8Sa9&h~=|CV#ohqyOHmnBRQFn7g1<X}4p2)_bGNl&S$4qjcr(DsP{j^F(CK0cGSj
zK9y6XD&;xftEJzZIjmf2_t4mYfUlLkOU#Fs96z9p3tgY5mXu^G*<<vfYuXnmyHl&}
zr2*y2$KAfY80`PD^3hYab^SiwtnB-Kr(f=(A|=Q+bU<<Fi^^_4CV>B{gvxs(%9L;C
zj6PQPWQpS3{N9j3e?Ov_7Dv=PJM|;wjf#L-a;M$O)Gk#MX4kx>JRR+KG=D&evj3}<
zZI9;dRgy3CU*D~7n_bGuX%!QO2X9l_^qaS#+p06lR4k2=!m_imvbOq_9m<xxLy7+N
zWlDBx-ytvG^|}(hYvA(TtzT0N1#?%8e)D+3_L&c-zk7dGLZ{+6Z+*D^73Ggp20T@A
zVuNx{djE~%L2H#)$F8*xp1Mehy7<7nku|E)k3WIPVNW2GfFbMwgxym7^p{GR68*1R
zZ+)f9!EB|g<I5GB9UmyuPnDhsne(1f&YwVxd@TLwV?%Z+>lZBXyD;c)%EFv2&OdL#
z$|FgwCLCY6Dd8A5ARjUs2!7c)cHNb2%KnM@&utm|s&d!adz06_vP;pmIbC^h={ri;
zs_#23wSJ<M?t9_DSpV#Il$3`be>!E}Yf5sJZS)(_FDt3<Y+Li$>%~g`yb@#eqzy{{
zwx2JYJLOen^^+5iANohBvM8hW+dc>OC|6>aj+;MsyW%_Qd1b`@&y{Z9#XaTw+fv16
zZdg#U^&6$@w|yq``Mh3v?C_hHTBqym`WX+LiO#vA%+0!U{lXO&6>~pZ?pnNftkgrg
zjNcq9MU!cLj)35q>xS=mbo1BBo$>v*7O%aeWaEf&Ql#`g9z;OI;$1(ezJJ)KjQ?co
zI~(3BQa&!b??Riw?<j$I&NwMd`aS;BIB8tO)q8qxdHRk^%J?>y9=xcYR|aog+jqsQ
z|5B3h596diO8>$O#!Do$`42q$8no@u+}&%gYEOvfOw`x5TR7dmy;sG<7yfyty~&4>
zyhn`$ibrDbWV4l99+7l2_0>p&G5*Mq>=|ndzE?h5qx2h<gzV)sA82R!HXMf8Po;lU
zKKGKve!*~l(f4sWJ9bV$HlwBiK0i_FA{RbbGV)2Ott=-$h`UdgX3P3ya6QFNt3<Xo
zt;SwND7bo&CpSa=EMu7q%#?8K0I5`7J)PGcKAqLBOLeQ=Ix?QtmH3^Do#u2p0cE^k
zxfF(1rlAlV+gj=xQ8$Bo?Ho?YL8hbZ2`ygbJ3p>_gO*f_Nwf{YmchUV0|~}y;BC-$
zLQS~vf^2{gPmg1d%CTC(b&Vc!J4)UFAg>g-Rm%Z>BpBD<98)WHb|tbjU`jeL<Yp!H
z5Yp0#y*i@PH3`+xrZ%8Y)Xg;*tI2K|)MOjng*90&WVmsqUCq-Tv<z*sjpk$!<u@sd
zhuIs9CVg-tscy0rvN2)N4Bns>n!P~_UQ6nm;n!Hc2E!`KbG%z_c+fTM@Ec~}ZnehZ
z{xk<crjbj&8`sTUv=$#ut-<_b&qnh{%vLW@!kv!F1KnJ-H&fD{Hxb``D@2@d-3(++
z!7?-EMkIvhu7|!kP3Y_eCC1#0@XaYrTklQh>qFNg?e%Ka6y3G8<qIWfoAFC%+mtmV
zT=u@)Bjb7u(3DtBDcltP8ivq%{NAP-^M+@yn(A~-LM&a7+jQ_&RhUJ>())6k>}#3f
zTGa&`Z?Gvd2t2<XA}^roaF3&MPcJh&7Q?uhbUhL@X#;4ETuqi!8-;pn;uo}yXbvRv
zX8WIJn<jOI=E!y3<v+<b&=Tk-biGk1>Nf>p+;llA_xCbx+?<NX;JI$ap-pl%OMfi`
zfq%K?#DW&wJO4jrAZ?BGSYDdd6~rZ1vvk+830c91+>lM6#p@YFg04B()S$Vx18!Dl
zXqH@!=%_84alkaHf%b38Fwm^~AiB|fGoFUtYJ(HL<65J$k=`H$xf<cG+a3_xLEHa(
z5sIeyt==?mq#0AV35}arUOeJi=CdaKl)eNRu2tV#0K%47uQ~!d(7jn*X|D#jo2gMt
zFuwm|u{6`Ezq>#}7Xh2p={j@~WAoN)*%}(UR(T@}5w1+`TfE8LNM0)E^^7Zi^9ddB
z3kmGizi9$rr$J#A|G91owW<5FNk)B9sJH&~){peXPz#M+J!XaeMmr5gnwRj#g>wsZ
zt~uPrux|*xOaFhM*`D4eNki>K{wB4&#&i^%1-!`?=*J%+bfDb8RBpC=T;ro9)!Ynf
zvCw~og*s!uAyV4_*XMx{=)*%`<F+qZi0k$`jGXc%?wck}X78p{KjeBhrT#1aZVGum
zeB)1eFYsVG2i$_|f+xK2yct(bmqyBY<8bUlQmKF4Y_^;XKQG~bc91*B)6auYbN%y<
zv$eAt3R+7Ykp6;%UvDdC;M&%b+26_;=Mn6;8R#L`E1)fYflX$*{L?3IKk*@+%Qx`d
z|6eUX8lS)eB)wXjqyGudEnO=Ak}M4n)cg$3QJn6N&3<lvH~7fS+wAj@?`Pmg+Tao2
z&pW(-iN3yO-0C*WTIiaI@mp;5UlF0SbYn|5w)nahU*|1OXu-93#}@C{;vN6Ldq;yh
zBG_t-dz|fznB6{>nxDmR>|-gqaf4$%LFQn$Gu7zEh#7yqT@G_UTg_4S9k!1mA0e4)
zw2dBg;0f7L)?q~wGN6jNB>5@bpfS%S8#F27amrk21K!zRN|P-U@ce9PgWCOEN1U`!
z8Y!*Be_1HKOuyhlc~8kD_dF?u$|&wB8P8fOWuWRO<>Ipb5*!7B=Pj3hgQ}mFOWt2D
zEzeP3{InxQ`BnSy&Fh?I4q<|?5~}~ei0O$<OV}<~EPtmWQ6^<{mkU)aM_iP28Dal6
z(j@qxET9G`$Yp{i`F_|74CS^ThGKn4osKfC<8Q%?AmIf>4)HM7uVPsX?z9+%y^_{Z
z&LPLov>e?-?W)76H-2;vPQ6Dmw>!YQMm&Ibb(Oj<JksI3%Es8THK^oF8jP`&7Z6r=
z${Dxk^6?{Tx^);8N_^WqzZ2U?6m(JXEI)j(9_9oboMv3TuJLF{ScC9icFP&F3L7+}
z4TUwl60?D>yp-idLzSmlKK&|#;$y*BUr&K5+p&DTg$Q(n4p{wPS(PY|iw)T77a1ys
zeb0!k_aTEkd>!t7ADWJrBuRe5D%Ua!(Y^;0BBL-!Vu>eR7SAtxlbQSk^%JP>N8uaf
z*q2C2M_h0~Cf|epjZqzJ=IA9j&miGt)>3UyiRGYNBfEXzL%vkc*FJ<7fgkzWS35a?
z;>hFL@$Sdg&=1*N-@nsr8mmM$nZ~mesg&f#h!mu11@UeLw1r>ip<rhDI;<P}mauNK
zr^b<U8!Sf{H8}>Wuu3iraP>ODGw*46pL<Vv+nadi6f~KR#kifUVloVGHWsIUz)&_m
z7Hw<D0+LmV@U+k5N!Tz%3io##{NO&T`5e&Wb3Xb=auCJ{#Cn4CjuTS6V_;aoc@elM
z0FiT5$ibU%Do={T8AUUqXgVvZDn`G<9RtxbkTq*v%SVaxAEdTxFp!q!0ZAHsstlZo
z6gzzA=g$u32=Uc9{WN>&*(6+BL&6Y3mI!~y7u?6?;qZriMGYzPE>ODbkptNf52v4I
zld^dFO*i{i+a|cik^{W(-XNsMX+bE8G#n|zVWmFZAGmMSaB}phjvoN1GjU^!IqWjN
zaEEjaPuLC<k{-)ac!DP3T$a+4*EcocR5Z(3)`jd;v}>Gdr~-A!rykUe&(5Mpa+%IE
zs?^rFFkCJ)&X<{EgU`JO8j05A)CqP<&TtRI5UW*@597_)8i-DkldhPh>I0DlTY`|?
z;WPmz7Uu;a6E=et2|5mNhLZ9iv=X%()uh9}lB1d&WnaJpx^Ow~DqeaIdd7GFu|eRz
zZlNghL0OOcP)0;ugEiCBPporFEM_e~1%n*oRLQiwlRlCzMaa^MESR5qKW~m5C6~^Z
zS~4Qx?9RsH)^$9bU-<+!=@Mn?ApPBzWio!RxK6d<r>VG|RY`B#H){Te-^CMOM$Saw
z2Ugm&Qd&KeWr=qhp{4#?Pg5&)u9Uo<7FEt8GjiM!_zM(g)${1x(m{{A{9Aq^V08yJ
zdP*Y+WQ$K>KQ9BSGLj9b_F24waG-7sipx4k^n3~T2#_2=oi~QJ0RiD#E$8B}a2Egn
zJQv4Nwu6Pg<q>2wuv_E0I7l*|qaOf1=URt&6%Jv;$^8&7FQ;z7xxab$GSUPYSEixu
z>dv?kDCB(vXOM^D!2_htsBoN^64M9k@#~Y(B<fvJCCSa2g4nHU-FU_(WFkFZJ#h-l
zrca~UWz$)9+YEhi*~>0BM&G(>Ih+e3M1Ds>Qsj88aLplGS&H~^m%Dc1XfSCeZ&~ob
zlhx*|=)DE0$cWgMo=3KT@g{V9IdMCjs|^X0;Gm0N?Wsp_GE;^xKcoVej739z*~t{q
z`07|R$oF%4Y(;;pj6)HvSO!lXhX(t86*ldw@dZENzm7vgS8S8Q@uhJnCVF=3t246I
zAC$?hPJFm^{AK0v*od7YU%IT^b!^+TVb5Py_OIGk99j1RzH2<{lkH2~6CC6U$7%$S
zpAUyH$OFj_U=aC1C<=GH0O4$MfUP-f!(j-4Ko6(Y0)$if0EUqt@P#{eXeWaFpj`s9
zX{ZheIOhWh;52Ti01y}e1cm^i0f5jT(1QQ~K_Gx2AV3fpASeJ36ay&E5doptIU*C;
z!UocT6#&A4B2Wkrq=W|HegYsc0VvWyBTgN*k~Cb}e=p!D3&l$&pw{@%tH|a%d;6m^
zhRp9;a`jatNv)=?#qw89-0F4oMz*_T&@h-*ND2Zp_n!y2*9Sq}>qA-dC#3k$9n!ps
zU>Bk++aZPhr1md0*zIv!++qD^+dUO2tYSRezh0Ga>O;tks+qeipUKK<(y8|>eUxP@
z=jhdmZ924%<KA7w7jC|n#uw!DOcXLuuMQp25o`{ZWil)uaI=PDAB1EVc4mrBm(qZl
zyXP)@8K<V9$M82RP&gd3g?Ii272vNfqd{_tRC44I6d)mMiiAmOE88QBeOS-;$nvJ?
zOUO%cJtdTtCf1LN8G1Zq4$7F7+n`am^?`~jwu$1nb%myGhM{y~J1eCDotwqTb{5^`
zLkX#5?ddP!oSEc=uow@ir)LNWO&NN;Xg<n7rH{B#+dQa~9~0E|b2xS3Y)VabuBW1S
zLfmr|+8<Dn#k#F-E`e-!BkL(qHRMS>+Z)agZ`I$R8&B0C{dAtGDin5mHIY)Jak`56
zurUS)p`q$4G)cPW0X-eYI{1#0ULL#C1e4G4g<LX0=OHPfDd2RbqcBwZ?cmwc2Ew_$
z9!V(<ak}9heLy(a^{$-9Y*7C@%?9(c#SpPPrN`59QAYc1oJz>!R;!o?%mLN&87F<V
z;AG5wmZ9AD&!e(QP(qGZn^nr`1y7IqL{KrrOzRjEK^MpK$?F+j46$XUo^+)=U%8vu
zi0zox<S-Hw$i;3H%5xct&T?O};R3r<vCooeum(+|Tc)=k;I)JR9$ef6&2xAPu#M*{
zR3gB4*!m{%AFChtnDrsgLI2TNAN5B_()h8*@IQHSO6U2!MO;|kkLk~oI?HkD^JpeI
zG)iA`<awCzMM?USv-yw;6=CAeiD*NB2N$NjK94c<xbH4pq{si>L#9~VLJ)~Q_(@!>
zNrzIFko2eQAiwd%&5(oqNLy)SRt4-w|G+0#qWjUx(R$qT1#|_?6O+x&8p8I8%eXr3
zJA`#AwS5MfpYLI+Wy?VVYO0r3Xu|Y3u+r01H9T+ni>Ik-){;8Q;mMT|Rv);0=}}fI
zQ&_*o>u}{xRNA1_sg;)0>_V+1)YikK9<4Ex=JJL;Od9(aZqjWSM5iV#U#)3@D%X%S
zH)+oOdT!GE`}NeM!9rs3QH?Klq@niotSn=}{ajHscy1UT`2wm!CrLf(;5?=Zb8vYZ
zgoDdkugAN#p(g_B_(TfXKsYfaDqi7C-Fqf(gC0LsLcBpauOKFb8-rq4KIurW|Eov$
zSzZNg+4}aHrg)FULe`LAmBsVRHgf@3SLw3eAa!W>t$08hi7pWh?|nd85B*6Ix@3mH
zISw)x!aO?3gD2C^w1u(v=%j6%c_$^j4g%6nN-ol5MbnE(ns<_6jOe6rk|N{&7Ri{q
zl`%ELLCVYNmnxDEesvRy#2Ysu^NKSl9H%Cb|Egvl$=xmD7gmP_YM`vhLv#lh1aELb
z)CU*DQgA^u2$wJJa~g%=nA3>-Pb#O7HCX%~I}rK7Wyb3YQJ7@HUlWim*oxAIP>3a<
zh)+ZwLNB-k9zrv?L|zP4K_|rJbok=yC>ZZ0!jco8$~K{}2%)cl5a+=~Rwy25hHpms
z*%>s@jJJ1x?L^`opAE9J8-^hqCUAJG{UygM5fgUa#-U;-ha2X_PO>w2)^0D_=bO`K
z=}-n+53?^R9rE~qgB;#<yFDG+slVMmvS|GTd+~64474}n4uvH?yMrU{w6}W6-#)Va
zUG{3Yjd9*(--jRXn*4_m_5u)K@ZAjVyxacFqBEw+kK)N&&>LZIg`6H-^7R$vf#<GX
z!h5eOtGjgmvR&O(99e{3l@K1CK#r~yh7@)Xh8KC!VInU&ROH1t5%v`0MC8S?D#n!<
zF1({!)7uW{0MV1dIQSImjq?kU6@O7k#&A_38EQ`Wv+@)QmjZC{DHJ(H$N_gVa0xlZ
zfAS&p5P6}e$cq6Z>WhH|cSUd+u#q$~NT7upLQ0{g$O|)w@fU!HlHngAv=-6`tN3B*
z44Fmwn~0KRbjHIU171lcDmI@%hHk9}#*W!EB>lG6HyHZcF894n)Q@Ds>9qfBPVIE=
z5fn3qJ7BjML`0Y&`-B6o-_0HH{5?#J)nx=S2bHp)p*+Hp(lYLVlX(T<HQgAv9B}R~
z4+s3x0q|th0dIR-b80;*4zh*>t9qWdyvqgP3f3PXb=dmZi=*gKkj7<yV}zOb%ue){
z=Q)u^eF-lt`R6X=E1_@?$L!G-(|Z|j*26I;yvrT)fxm;`)G=2c)^tEM<s{7=^9ao`
zk0+e~0ZFisSFE1BVD%j8m$ie*Rz>&KE+E4scAw;M8gS-0WDj`wu_=$telTOo<Fg(o
z-LmN%Qq_;+=h0h;v@v>h@wsC4zTkGo^b?E^yTPPEbwve3)qOx9RSadF5kYmx2_94z
zaw=hZw_4DPEIdESIMtsPoQ(57Vkq~Z+H!^m)!R6|;92QjpY)rUHZT3im<YOZp0DBS
zjR>mi1KErbbWuSpKRlTDXBJdfS29+5kc;)=fhvZw^6<fal?2sMAL;QuS4B`=%4-S6
z9$ef6P5#DH(D^*Sx{mbucWf7S{xOe{Zu%I828{Ifu{EBxlX*Vt0wFm<-#pq|NA5N)
zpXhOi?~v+nw#v5H+4zoOFejM9p7tf5#FKx4_JeaQ<e$wBe-F+var!>?YA9af&+~ol
zN6O*i<vINt(8qB;L7cw7eGFXLV}k92XiNP$!#4KUUxRA^7oaVN16hSJ5u9!yXJp~@
z13B!dl`}q8t^Q+FgH9TM$Yvj$(^(UM3DQ}UkZIn9=i@l!`5d7!r_brb%Xwbco#(k_
zbs;?8mBV%bv&pYBB*@j(-gU;gi8YH7Z<8+PjT-S};%!ZCOR);N8T1OYv+OMmOx%@Q
zEZEL8kJc>Ax~?0sOg404TVAn;W@6UYvAhp0oIy|Ckd>RT7sq$v^?URBU9>#wf~W@c
zw{Ska?5nb0fJ>MpmRA&#3W1_|h3;B8(>TW7k6xm_U@Y=coL=BNb3RcVTJ70P!p@qe
ztj{}XmSFivu0dOF!uEX3gm3}EIKvP->k`q|0=DL}BgB4Z*rcpucZgXKo^7Z3umF2r
zubp4KxTu|d{US*C+X)pZ>R|u!ML_u3p<37^t5)6{$a{6RQ?g?Lnr6doriCAmc=R?M
zzplOWNX8$N*vuu5<xoi6Q<I#{Z_($W+!DRHCBnG?!cSPZ>1`Zu)ykP%omu_tCB1n;
z6lc_n=ezQGEAlbCK`{us^Dc?y^bKKe-ZedVxsCG|)3OicpZ!UEcV5tk_jxzYxGxtV
zhQofGt}7R$ABVSadg1Y+v~ud<#&pr@>*&>4<4^nz%s;Cp0GpSQT;e`_dUen=U<`!I
z?!_5XM?#8(1w>bL=k-H*ePO{iJP#J8^S?DmL~#r7HGnuKrh&+JVstrXPT!H|g$obo
z^r70ASw*LJM1wJ9;d!&YUpAe;Ry!jUqqsBgJAclo3$NhE>2A?x0h<-V5BT$X-MD~e
zK8t$T-x~DUD=`m0sm%)3-aAT3{72*k3bqNGyi1SI&P0a;+}2GI(9~R@0Zo96VRe_6
zUBK^#>SrU%KVk8uAzEB5*>erxtWtK7?Wr(*sD!2yrgV40Z7COlL&$e}9b6uPLunli
z95T@GuRYfJl<&25e*Ps^PolMTx|ewVz-3tHm)Gilhvrx7afYNnj|xeGD&<BbluG!w
zhf>v-nNsN&HKlUv-IUr^EtD!HOsP^HffNh<qVzuejRh0g51?Bvc@hhw3K4-UW=i*A
zC=b=kkq@32sLyCWomUVFx-o!Qk;OsFu6WpITPqTyu{fy0z?h5b^*nF!hcFKB-H(Q;
z?b^koN}|Tw!Nf>4pK&9(Q3@k1^Y?*A2)yS}g1-+w9IDSi2YjJD7M+`XTwP!HnXpnl
z;YzLK7DfuL^FxG{rU!E?RR=Sv^&y&-4upEKQd9`HQYo(>)O2IuvQkZekNblcQ<!F@
zSblpbsweY&R(r4#eTAqOA7HL~iDPx>ZlmRMHvX#c;b8sesQybaZ(JONo;B?~thcN^
cx1J?jvtC>WH|y!YBG#i1a!g<8>4TjA0cPw)_W%F@

diff --git a/topostats/processing.py b/topostats/processing.py
index d2dd8d3948..b5f048c84d 100644
--- a/topostats/processing.py
+++ b/topostats/processing.py
@@ -1103,6 +1103,8 @@ def process_scan(
         grains_config=grains_config,
     )
 
+    topostats_object["grain_crops"] = image_grain_crops
+
     if image_grain_crops.above is not None or image_grain_crops.below is not None:
         # Grainstats :
         grainstats_df, height_profiles = run_grainstats(

From cec53015af9ede6d23515a4fd1237080a08089f0 Mon Sep 17 00:00:00 2001
From: SylviaWhittle <sylviwhittle@gmail.com>
Date: Fri, 20 Dec 2024 15:03:38 +0000
Subject: [PATCH 33/34] Add more rigorous GrainCrop setters and getters

---
 topostats/grains.py | 79 +++++++++++++++++++++++++++++++++------------
 1 file changed, 59 insertions(+), 20 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 0816e0f024..4309ae99d2 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -58,34 +58,69 @@ class GrainCrop:
         Bounding box of the crop including padding.
     """
 
-    image: npt.NDArray[np.float32]
-    mask: npt.NDArray[np.bool_]
-    padding: int
-    bbox: tuple[int, int, int, int]
+    _image: npt.NDArray[np.float32]
+    _mask: npt.NDArray[np.bool_]
+    _padding: int
+    _bbox: tuple[int, int, int, int]
 
     def __post_init__(self):
         """
         Validate the data makes sense.
         """
-        # If image is not square
-        if self.image.shape[0] != self.image.shape[1]:
-            raise ValueError(f"Image is not square: {self.image.shape}")
+        self.image = self._image
+        self.mask = self._mask
+        self.padding = self._padding
+        self.bbox = self._bbox
 
-        # If first two dimensions of mask are not the same as the image
-        if self.mask.shape[0] != self.image.shape[0] or self.mask.shape[1] != self.image.shape[1]:
-            raise ValueError(f"Mask dimensions do not match image: {self.mask.shape} vs {self.image.shape}")
+    @property
+    def image(self) -> npt.NDArray[np.float32]:
+        """Return the image."""
+        return self._image
 
-        if self.padding < 1:
-            raise ValueError(f"Padding must be >= 1, but is {self.padding}")
+    @image.setter
+    def image(self, value: npt.NDArray[np.float32]):
+        if value.shape[0] != value.shape[1]:
+            raise ValueError(f"Image is not square: {value.shape}")
+        self._image = value
 
-        if len(self.bbox) != 4:
-            raise ValueError(f"Bounding box must have 4 elements, but has {len(self.bbox)}")
+    @property
+    def mask(self) -> npt.NDArray[np.bool_]:
+        """Return the mask."""
+        return self._mask
 
-        # if bbox is not square
-        if self.bbox[2] - self.bbox[0] != self.bbox[3] - self.bbox[1]:
+    @mask.setter
+    def mask(self, value: npt.NDArray[np.bool_]):
+        if value.shape[0] != self.image.shape[0] or value.shape[1] != self.image.shape[1]:
+            raise ValueError(f"Mask dimensions do not match image: {value.shape} vs {self.image.shape}")
+        self._mask = value
+
+    @property
+    def padding(self) -> int:
+        """Return the padding."""
+        return self._padding
+
+    @padding.setter
+    def padding(self, value: int):
+        if not isinstance(value, int):
+            raise ValueError(f"Padding must be an integer, but is {value}")
+        if value < 1:
+            raise ValueError(f"Padding must be >= 1, but is {value}")
+        self._padding = value
+
+    @property
+    def bbox(self) -> tuple[int, int, int, int]:
+        """Return the bounding box."""
+        return self._bbox
+
+    @bbox.setter
+    def bbox(self, value: tuple[int, int, int, int]):
+        if len(value) != 4:
+            raise ValueError(f"Bounding box must have 4 elements, but has {len(value)}")
+        if value[2] - value[0] != value[3] - value[1]:
             raise ValueError(
-                f"Bounding box is not square: {self.bbox}, size: {self.bbox[2] - self.bbox[0]} x {self.bbox[3] - self.bbox[1]}"
+                f"Bounding box is not square: {value}, size: {value[2] - value[0]} x {value[3] - value[1]}"
             )
+        self._bbox = value
 
 
 def validate_full_mask_tensor_shape(array: npt.NDArray[np.bool_]) -> npt.NDArray[np.bool_]:
@@ -121,7 +156,7 @@ class GrainCropsDirection:
     """
 
     crops: dict[int, GrainCrop]
-    full_mask_tensor: npt.NDArray[np.bool_]
+    _full_mask_tensor: npt.NDArray[np.bool_]
 
     def __post_init__(self):
         """
@@ -962,7 +997,9 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
+                vetting_criteria[1:]
+                for vetting_criteria in class_size_thresholds
+                if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1689,7 +1726,9 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(
+                flat_region.bbox
+            )  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(

From 7f75723c08d03f831bdbe6c5bddc1b8a6ae740ac Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 20 Dec 2024 15:30:29 +0000
Subject: [PATCH 34/34] [pre-commit.ci] Fixing issues with pre-commit

---
 topostats/grains.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/topostats/grains.py b/topostats/grains.py
index 4309ae99d2..464663930a 100644
--- a/topostats/grains.py
+++ b/topostats/grains.py
@@ -997,9 +997,7 @@ def vet_class_sizes_single_grain(
                 continue
 
             lower_threshold, upper_threshold = [
-                vetting_criteria[1:]
-                for vetting_criteria in class_size_thresholds
-                if vetting_criteria[0] == class_index
+                vetting_criteria[1:] for vetting_criteria in class_size_thresholds if vetting_criteria[0] == class_index
             ][0]
 
             if lower_threshold is not None:
@@ -1726,9 +1724,7 @@ def extract_grains_from_full_image_tensor(
 
             # Crop the tensor
             # Get the bounding box for the region
-            flat_bounding_box: tuple[int, int, int, int] = tuple(
-                flat_region.bbox
-            )  # min_row, min_col, max_row, max_col
+            flat_bounding_box: tuple[int, int, int, int] = tuple(flat_region.bbox)  # min_row, min_col, max_row, max_col
 
             # Pad the mask
             padded_flat_bounding_box = pad_bounding_box(