TopoStats classes internally and for writing/reading HDF5

[ns-rse]

Closes #1102
Closes #1143

This PR (draft for now) is the logical extension of the GrainCrops GrainCropDirections and ImageGrainCrops introduced by @SylviaWhittle in #1022 and switches to using the TopoStats class @ns-rse introduced in #1145 for handling images and the derived datasets (arrays) such that the unit of interest is individual grains.

It is at the moment far from complete as the checklist shows below but because of the large amount of changes and reorganisation I was keen to share it in stages. The tests for each commit pass (thanks pytest-testmon 😀 ) but until all steps are complete the integration tests (tests_processing.py and tests_run_modules.py won't pass, I'm working on them as I go through each class).

Its perhaps worth reading the commit messages for the individual commits for a little more information on the re-organisation that has been done so far.

• LoadScans
• Filters
• Grains
• GrainStats
• DisorderedTracing
• NodeStats
• OrderedTracing
• Splining
• Curvature

Of note...

Shared Methods

Some methods from Grains were used @staticmethod from GrainCrops and have been moved out to utils. I've setup some skeleton tests for these but they fail (get 2x5 arrays back when I would have expected 5x5 arrays from flattening 5x5x3 arrays). Couldn't see any existing tests for these.

Documentation

I intend to document the class structures and in turn the HDF5 format these are written to.

Syrupy

I've closed #1143 as the key test which used .pkl's and required manually updating has been addressed. I think we could still switch all tests to use syrupy (see #1152).

AFMReader

These changes will also require modifications to AFMReader and perhaps moving LoadScans over but I'm wary of introducing a circular dependency and have already discussed with @SylviaWhittle this. We felt that perhaps AFMReader should only load files and return dictionaries. Re-constructing these to TopoStats / ImageGrainCrops / GrainCropsDirection / GrainCrop should be the domain of TopoStats. I don't think this should be a problem for the Napari plug-in as it could import and use whatever it needs from either.

---

Before submitting a Pull Request please check the following.

• Existing tests pass.
• Documentation has been updated and builds. Remember to update as required...
• Pre-commit checks pass.
• New functions/methods have typehints and docstrings.
• New functions/methods have tests which check the intended behaviour is correct.

[SylviaWhittle]

This is looking very good!

I've checked this out locally and just followed the threads of where data goes, and it looks very nice.

One thing, is to maybe move things out of / don't add more to utils.py given that IIRC we are wanting to eventually eliminate it? Perhaps a grain_handling.py or something?

I'll have a further look next week with Laura but looking great so far 👍

[SylviaWhittle]

Oops meant to merely comment, not approve, sorry!

[ns-rse]

Now have run_nodestats working with TopoStats objects, haven't yet fixed all the unit tests against this but pushing on with ordered_tracing, splining and curvature and will address all failing unit tests once these are done.

[ns-rse]

I'm revisiting Nodestats and have some questions about the various nested dictionaries I was wondering if anyone might
know the answers to ( @Laura Wiggins I think you worked on this most with Max so maybe you have some insight.

Matched and Unmatched Branches

1. There are Matched and Unmatched branches. At the moment I've created a class MatchedBranch to replace the former
   but have left the later as a dictionary unmatched_branches. As far as I can tell these have some overlapping
   attributes/values, they both have angles but MatchedBranches (nee matched_branches) also includes
   ordered_coords, heights, fwhm. For simplicity I think it would make sense to have a Branch class and not
   populate the additional attributes where the branch is not matched (and perhaps have a bool attribute to indicate
   whether it is matched or not).

2. The MatchedBranch class mirrors the structure of the existing matched_branch which included
   matched_branch["fwhm"]. Reading the typehint and docstring of nodeStats.analyse_node_branches() the fwhm is
   "full-width half maximum of the branches" and is meant to be npt.NDArray[np.number] but looking at the tests its
   actually a dictionary with the following (which comes from tests/resources/tracing/nodestats/catenane_node_0_matched_branches_analyse_node_branches.pkl)

matched_branch["fwhm"]={
    'fwhm': np.float64(40.314950876985606),
    'half_maxs': [np.float64(-18.432397834711107),
                  np.float64(21.882553042274495),
                  np.float64(2.821882638708631e-09)],
    'peaks': [np.int64(93),
              np.float64(0.0),
              np.float64(3.819342748549839e-09)]
}

My queston here stems from the fact that the top-level key fwhm doesn't describe what the dictionary contains and so
I'm wondering if the nesting here can simply be removed and so the actual fwhm (the number 40.314950876985606 in
the above example) becomes an attribute of a [Matched]Branch along with half_maxs and peaks? Obviously these would
be empty/None attribute in branches that aren't matched?

[ns-rse]

Answers after chatting with @SylivaWhittle

1. Keep MatchedBranches and have UnmatchedBranches as a separate class.

2. Everything nested within fwhm is required (i.e. fwhm / half_maxs / peaks) and it is fine to bump these up to
   be attributes of MatchedBranches as long as there is clear nomenclature that all three pertain to "full-width
   half-maximum".

Also discussed...

Having classes with attributes isn't too dissimilar to dictionaries with key/value pairs, but with the added benefit
that the types are clearly defined (no more 🦆 typing in the code base). Because the code base is growing it is
becoming harder for any one to keep it all in their head and reason through all steps and so people will be working on
one small feature that has to fit into a bigger picture. Ensuring datatypes are consistent as they are passed
around/through the analysis pipeline makes development considerably easier as errors manifest themselves earlier in the
development process and can be corrected (I've found this to be the case so far with some types being wrong in the tests
as I've eased Pydantic dataclass validation into the current branch I'm working on).

[ns-rse]

I have now added config as an attribute to TopoStats class which means the configuration which is loaded when
loading files follows the TopoStats object around as it passes between the processing classes
(i.e. Filters/Grains/GrainStats/DisorderedTracing/NodeStats/OrderedTracing/Splining/Curvature).
This makes life a lot easier as we don't have to pull apart the config dictionary (loaded from
topostats/default_config.yaml or user provided option) to pass into each of the processing classes as it sits adjacent
to the image data itself.
We will have to remove the options and populate the class attributes from the TopoStats.config attribute.

For example we currently have...

class Filters:
    def __init__(
        self,
        topostats_object: TopoStats,
        row_alignment_quantile: float = 0.5,
        threshold_method: str = "otsu",
        otsu_threshold_multiplier: float = 1.7,
        threshold_std_dev: dict | None = None,
        threshold_absolute: dict | None = None,
        gaussian_size: float = None,
        gaussian_mode: str = "nearest",
        remove_scars: dict = None,
    ):
        self.topostats_object = topostats_object
        self.image = topostats_object.image_original
        self.filename = topostats_object.filename
        self.pixel_to_nm_scaling = topostats_object.pixel_to_nm_scaling
        self.gaussian_size = gaussian_size
        self.gaussian_mode = gaussian_mode
        self.row_alignment_quantile = row_alignment_quantile
        self.threshold_method = threshold_method
        self.otsu_threshold_multiplier = otsu_threshold_multiplier
        # Convert to lists since the thresholding function expects lists of thresholds but
        # we don't want to use more than one value for the filters step.
        if threshold_std_dev is None:
            threshold_std_dev = {"above": 1.0, "below": 1.0}
        else:
            self.threshold_std_dev = {
                "above": [threshold_std_dev["above"]],
                "below": [threshold_std_dev["below"]],
            }
        if threshold_absolute is None:
            threshold_absolute = {"above": 1.0, "below": 10.0}
        else:
            self.threshold_absolute = {
                "above": [threshold_absolute["above"]],
                "below": [threshold_absolute["below"]],
            }
        self.remove_scars_config = remove_scars
        self.images = {
            "pixels": self.image,
            "initial_median_flatten": None,
            "initial_tilt_removal": None,
            "initial_quadratic_removal": None,
            "initial_scar_removal": None,
            "initial_zero_average_background": None,
            "masked_median_flatten": None,
            "masked_tilt_removal": None,
            "masked_quadratic_removal": None,
            "secondary_scar_removal": None,
            "scar_mask": None,
            "mask": None,
            "final_zero_average_background": None,
            "gaussian_filtered": None,
        }
        self.thresholds = None
        self.medians = {"rows": None, "cols": None}
        self.results = {
            "diff": None,
            "median_row_height": None,
            "x_gradient": None,
            "y_gradient": None,
            "threshold": None,
        }

This would become...

EDIT 2025-11-10 : On reflection I think it could be foolish to remove the class options up-front. Instead we should populate them from the topostats_object.config attribute if they are not provided.

class Filters:
    def __init__(
        self,
        topostats_object: TopoStats,
        row_alignment_quantile: float = 0.5,
        threshold_method: str = "otsu",
        otsu_threshold_multiplier: float = 1.7,
        threshold_std_dev: dict | None = None,
        threshold_absolute: dict | None = None,
        gaussian_size: float = None,
        gaussian_mode: str = "nearest",
        remove_scars: dict = None,
    ):
        self.topostats_object = topostats_object
        filter_config = self.topostats_object["config"]["filter"]
        self.image = topostats_object.image_original
        self.filename = topostats_object.filename
        self.pixel_to_nm_scaling = topostats_object.pixel_to_nm_scaling
        self.gaussian_size = filter_config["gaussian_size"] if gaussian_size is None else gaussian_size
        self.gaussian_mode = filter["gaussian_mode"] if gaussian_mode is None else gaussian_mode
        self.row_alignment_quantile = filter["row_alignment_quantile"] if row_alignment_quantile is None else row_alignment_quantile
        self.threshold_method = filter["threshold_method"]  if threshold_method is None else threshold_method
        self.otsu_threshold_multiplier = filter["otsu_threshold_multiplier"]  if otsu_threshold_multiplier is None else otsu_threshold_multiplier
        # Convert to lists since the thresholding function expects lists of thresholds but
        # we don't want to use more than one value for the filters step.
        if threshold_std_dev is None:
            self.threshold_std_dev = {
                "above": [filter_config['threshold_std_dev["above"]']],
                "below": [filter_config['threshold_std_dev["below"]']],
            }
        else:
            self.threshold_std_dev = threshold_std_dev
        if threshold_absolute is None:
            self.threshold_absolute = {
                "above": [filter_config['threshold_absolute["above"]']],
                "below": [filter_config['threshold_absolute["below"]']],
            }
        else:
            self.threshold_absolute = threshold_absolute
        self.remove_scars_config = filter_config["remove_scars"] if remove_scars is None else remove_scars
        self.images = {
            "pixels": self.image,
            "initial_median_flatten": None,
            "initial_tilt_removal": None,
            "initial_quadratic_removal": None,
            "initial_scar_removal": None,
            "initial_zero_average_background": None,
            "masked_median_flatten": None,
            "masked_tilt_removal": None,
            "masked_quadratic_removal": None,
            "secondary_scar_removal": None,
            "scar_mask": None,
            "mask": None,
            "final_zero_average_background": None,
            "gaussian_filtered": None,
        }
        self.thresholds = None
        self.medians = {"rows": None, "cols": None}
        self.results = {
            "diff": None,
            "median_row_height": None,
            "x_gradient": None,
            "y_gradient": None,
            "threshold": None,
        }

[ns-rse]

Copying this from Slack chat so its recorded...

2025-12-05

I've been making some decent progress on my refactoring this week and have been working on making sure all the plots are correctly generated at each stage.

However, I've got some questions about the plots for disordered tracing and would appreciate some feedback/thoughts.

Nodestats

On the main branch I run processing against good 'ol minicircle.spm and enable plotting just for nodes. This is performed by code in processing.py starting on line 593 there are "whole image plots" made of which there are....

• minicircle_above_nodes.png
• 25-convolved_skeleton.png
• 26-node_centres.ong

Astute readers will notice that these are essentially identical images and do not have nodes, convolved skeletons nor node centres plotted on them, suggesting we have problems (on the main branch) with this plotting.

Questions

• Does anyone use these whole image plots from NodeStats?
• Does anyone use whole image plots from other stages or is it individual grain plots that people actually look at?

I would expect the individual plots of grains to be of greater interest as the spatial relationship of grains within an image is not something I'm yet to hear people being interested in (perhaps a minor exception is excluding grains that touch the edge of images but that is because they can't be properly analysed).

Ordered traces also suffer from the same afflicition, the whole image plot is missing skeletons of the back bone.

Proposed Solution

Rather than spend time trying to fix this I am inclined to remove the generation of these plots at this stage in the refactoring and saying adieu to these plots

1. It reduces the amount of code, particularly redundant code that isn't required.
2. These are intermediary files, at the end we do still get whole image plots of spliend traces produced correctly (and of curvature).

Would this cause problems for anyone?

2025-12-10

Further I've discovered the "whole image" plots made on the main branch during ordered tracing suffer from a similar problem as the above, although at first it isn't apparent.

These are the three plots produced with...

plotting:
  image_set:
    - ordered_tracing

..and if you zoom in on the grains it is just possible to make out a skeleton.

Thus for the time being I'll make sure these are generated but still wonder what the utility of them actually is given they are barely legible and suspect the cropped plots are the one that are more useful.

It is also worth noting that the Mask label at the top right of the plots is incorrect.