Adds topological features into better tracing

.gitignore

@@ -57,3 +57,8 @@ topostats/_version.py
 
 # default output directory, often common from testing
 output/
+
+# Include all files in tests and all subdirectories except processed and __pycache__
+!tests/**
+tests/resources/processed/
+__pycache__/

pyproject.toml

@@ -36,14 +36,15 @@ keywords = [
 ]
 requires-python = ">=3.9"
 dependencies = [
+  "AFMReader",
   "h5py",
   "igor2",
   "keras",
   "matplotlib",
   "numpy",
   "pandas",
-  "pySPM",
   "pyfiglet",
+  "pySPM",
   "pyyaml",
   "ruamel.yaml",
   "schema",
@@ -53,7 +54,7 @@ dependencies = [
   "skan",
   "snoop",
   "tifffile",
-  "AFMReader",
+  "topoly",
   "tqdm",
   "tensorflow",
 ]

tests/resources/tracing/disordered_tracing/catenanes_disordered_tracing_grainstats.csv

@@ -0,0 +1,3 @@
+,image,grain_number,grain_endpoints,grain_junctions,total_branch_lengths
+0,test_image,0,0,14,575.5249825836854
+1,test_image,1,0,12,574.7857998943139

tests/resources/tracing/disordered_tracing/catenanes_disordered_tracing_stats.csv

@@ -0,0 +1,25 @@
+,image,grain_number,branch_distance,branch_type,connected_segments,mean_pixel_value,stdev_pixel_value,min_value,median_value,middle_value
+0,test_image,0,29.857677897827887,2,"[1, 2, 6]",2.6683043052890536,0.3576668425125195,1.993486285613394,2.6464969869174926,2.6743964509318303
+1,test_image,0,28.644860587199464,2,"[0, 2, 6]",2.728953130422273,0.21545195325432523,2.5506513793095915,2.6838064728094184,2.6579579145357006
+2,test_image,0,1.8682724368761408,2,"[0, 1, 3, 4]",4.067541277649083,0.16905835685823303,3.7890382980796193,4.138536286707056,4.13958834095104
+3,test_image,0,146.60562184380709,2,"[2, 4, 9]",2.6693913722865377,0.2947637179915079,1.2200007753773205,2.7080175948240655,2.9115862857598005
+4,test_image,0,26.170179495009112,2,"[2, 3, 9]",2.834959559756504,0.3436409841063358,2.535374801513602,2.7607877845412387,2.7651337225310604
+5,test_image,0,224.58506419260596,2,"[6, 7, 10, 11]",2.700995491724554,0.2653435468736685,1.787397162651758,2.7330374920770115,2.2573671026451536
+6,test_image,0,6.9504086553142095,2,"[0, 1, 5, 7]",3.626373814193468,0.6406710041895474,2.588846104614013,3.84195778773634,3.84195778773634
+7,test_image,0,31.203269242513674,2,"[5, 6, 8, 11]",2.843560036090906,0.5572924823368127,2.1307434183593315,2.705164120987648,2.75295411300335
+8,test_image,0,31.286996805637532,2,"[7, 9, 10, 11]",2.7443198761048193,0.3822328907364185,2.4149214282007074,2.653495727717319,2.4706619619066386
+9,test_image,0,5.570136218438069,2,"[3, 4, 8, 10]",2.7308714736136963,0.049010037214155706,2.6525934702850984,2.7378802570467666,2.749457125278597
+10,test_image,0,38.18835899001824,2,"[5, 8, 9, 11]",2.7864423529784395,0.15142352544214144,2.445066671627662,2.7710339081422832,2.6967078592376783
+11,test_image,0,4.59413621843807,2,"[5, 7, 8, 10]",4.126603901043527,0.37133617585304246,3.351114087306773,4.238726491399509,4.238726491399509
+0,test_image,1,37.70035899001824,2,"[1, 2, 3, 6]",2.7732793700561382,0.13491385165238218,2.445066671627662,2.7690167332032654,2.7141994425994813
+1,test_image,1,223.89492797416787,2,"[0, 2, 4, 9]",2.69655183205219,0.26110756906553556,1.8203021919659572,2.7314109425842963,2.2743843560760055
+2,test_image,1,5.08213621843807,2,"[0, 1, 3, 4]",4.029037066510506,0.4638075153095662,3.0770475331177405,4.18502108746829,4.18502108746829
+3,test_image,1,31.001133024075603,2,"[0, 2, 4, 6]",2.746190739871742,0.4020468170580353,2.443926567140926,2.6515637708741715,2.443926567140926
+4,test_image,1,31.405405460951744,2,"[1, 2, 3, 9]",2.8422583900088516,0.5622599874985189,2.1307434183593315,2.694237316611716,2.75295411300335
+5,test_image,1,147.49789428068328,2,"[6, 7, 11]",2.6728103603722104,0.2940846131940484,1.2200007753773205,2.712208521357603,2.9467266502403127
+6,test_image,1,5.77227243687614,2,"[0, 3, 5, 7]",2.7225654392347494,0.051361792194696554,2.6525934702850984,2.7127873939436142,2.7180535973290754
+7,test_image,1,25.682179495009112,2,"[5, 6, 11]",2.85129816853438,0.36309882085275597,2.535374801513602,2.7617733887534093,2.770465264733052
+8,test_image,1,27.870996805637535,2,"[9, 10, 11]",2.7169147220384864,0.1798385572282301,2.5506513793095915,2.6863743944825202,2.6579579145357006
+9,test_image,1,6.9504086553142095,2,"[1, 4, 8, 10]",3.679401125994582,0.5950886464826225,2.6485713386465375,3.84195778773634,3.84195778773634
+10,test_image,1,29.369677897827888,2,"[8, 9, 11]",2.640305654101202,0.3106025836163093,1.993486285613394,2.648283332587627,2.682776277975024
+11,test_image,1,2.558408655314211,2,"[5, 7, 8, 10]",3.974846970762266,0.23923528279875325,3.6040697432150033,4.075122442799481,4.075122442799481

tests/resources/tracing/disordered_tracing/rep_int_disordered_tracing_grainstats.csv

@@ -0,0 +1,2 @@
+,image,grain_number,grain_endpoints,grain_junctions,total_branch_lengths
+0,test_image,0,0,13,968.5225788725928

tests/resources/tracing/disordered_tracing/rep_int_disordered_tracing_stats.csv

@@ -0,0 +1,13 @@
+,image,grain_number,branch_distance,branch_type,connected_segments,mean_pixel_value,stdev_pixel_value,min_value,median_value,middle_value
+0,test_image,0,172.69207377569276,2,"[1, 2, 3, 8]",2.147391874570661,0.18580255309380492,1.4117732114082249,2.147881683634531,2.193944711921703
+1,test_image,0,338.03541679389866,2,"[0, 2, 7, 11]",2.1261598928571943,0.1841696602430932,1.127883594233125,2.1329128300921694,2.236342991231205
+2,test_image,0,0.6901362184380704,2,"[0, 1, 3, 4]",3.5129004230159246,0.2441665496774352,3.2687338733384856,3.5129004230159246,3.5129004230159246
+3,test_image,0,75.41508335877963,2,"[0, 2, 4, 8]",2.1781050277619003,0.22952699769811186,1.4117732114082249,2.178401374885392,2.2734312397952836
+4,test_image,0,31.893405460951744,2,"[2, 3, 5, 6]",2.1378544765518392,0.338975206263122,1.7555674668849865,2.094437231758162,2.2880045741720902
+5,test_image,0,51.51744873752279,2,"[4, 6, 7]",2.233921716584725,0.23610837345626215,1.7265179765777217,2.206985179983061,2.0220519901178635
+6,test_image,0,1.1781362184380704,2,"[4, 5, 7]",3.4335547526985297,0.29158462263860985,3.0268131230204918,3.578160584306969,3.6956905507681284
+7,test_image,0,153.8275289019402,2,"[1, 5, 6, 11]",2.1708118960130536,0.23399433389250326,1.32021443851052,2.1773699140297205,2.1908233129216366
+8,test_image,0,58.753721174398926,2,"[0, 3, 10, 11]",2.0020290183250613,0.3370066648426065,0.8859157914742134,2.08256912431822,2.247454035474802
+9,test_image,0,38.962222771580166,3,[10],2.2735993718579453,0.17804635219989257,2.0015543670262534,2.2595377460435833,2.2032782158757636
+10,test_image,0,0.488,2,"[8, 9, 11]",2.875087133879444,0.009881348191510873,2.865205785688025,2.875087133879444,2.875087133879444
+11,test_image,0,45.069405460951735,2,"[1, 7, 8, 10]",2.1854926308595277,0.22637773949873832,1.289281969607514,2.227074000531621,2.227074000531621

tests/resources/tracing/nodestats/catenanes_nodestats_grainstats.csv

@@ -0,0 +1,3 @@
+,image,grain_number,num_crossings,avg_crossing_confidence,min_crossing_confidence
+grain_0,test_image,0,4,0.4013589828832889,0.2129989376767838
+grain_1,test_image,1,4,0.3441057054647598,0.17063184531586506

tests/resources/tracing/nodestats/rep_int_nodestats_grainstats_no_pair_odd_branches.csv

@@ -0,0 +1,2 @@
+,image,grain_number,num_crossings,avg_crossing_confidence,min_crossing_confidence
+grain_0,test_image,0,5,0.07082753253520613,0.01059564637975774

tests/resources/tracing/nodestats/rep_int_nodestats_grainstats_pair_odd_branches.csv

@@ -0,0 +1,2 @@
+,image,grain_number,num_crossings,avg_crossing_confidence,min_crossing_confidence
+grain_0,test_image,0,5,0.07082753253520613,0.01059564637975774

tests/resources/tracing/ordered_tracing/catenanes_ordered_tracing_grainstats.csv

@@ -0,0 +1,3 @@
+,image,grain_number,num_mols,writhe_string
+grain_0,catenane,0,2,++++
+grain_1,catenane,1,2,+-++

tests/resources/tracing/ordered_tracing/catenanes_ordered_tracing_molstats.csv

@@ -0,0 +1,5 @@
+,image,grain_number,molecule_number,circular,topology,topology_flip,processing
+0,catenane,0,0,True,4^2_1,2^2_1,nodestats
+1,catenane,0,1,True,4^2_1,2^2_1,nodestats
+2,catenane,1,0,True,2^2_1,0_1U0_1,nodestats
+3,catenane,1,1,True,2^2_1,0_1U0_1,nodestats

tests/resources/tracing/ordered_tracing/rep_int_ordered_tracing_grainstats.csv

@@ -0,0 +1,2 @@
+,image,grain_number,num_mols,writhe_string
+grain_0,replication_intermediate,0,3,---

tests/resources/tracing/ordered_tracing/rep_int_ordered_tracing_molstats.csv

@@ -0,0 +1,4 @@
+,image,grain_number,molecule_number,circular,topology,topology_flip,processing
+0,replication_intermediate,0,0,False,linear,linear,nodestats
+1,replication_intermediate,0,1,False,linear,linear,nodestats
+2,replication_intermediate,0,2,False,linear,linear,nodestats

tests/resources/tracing/splining/catenanes_splining_grainstats.csv

@@ -0,0 +1,3 @@
+,image,grain_number,total_contour_length,average_end_to_end_distance
+grain_0,catenane,0,1113.3149766322022,0.0
+grain_1,catenane,1,1113.4528311181873,0.0

tests/resources/tracing/splining/catenanes_splining_molstats.csv

@@ -0,0 +1,5 @@
+,image,grain_number,molecule_number,contour_length,end_to_end_distance
+0,catenane,0,0,846.6004494307274,0
+1,catenane,0,1,266.7145272014749,0
+2,catenane,1,0,846.5829590713623,0
+3,catenane,1,1,266.869872046825,0

tests/resources/tracing/splining/rep_int_splining_grainstats.csv

@@ -0,0 +1,2 @@
+,image,grain_number,total_contour_length,average_end_to_end_distance
+grain_0,replication_intermediate,0,1773.7493902268593,165.97964434071082

tests/resources/tracing/splining/rep_int_splining_molstats.csv

@@ -0,0 +1,4 @@
+,image,grain_number,molecule_number,contour_length,end_to_end_distance
+0,replication_intermediate,0,0,748.2001882005334,167.88686666919483
+1,replication_intermediate,0,1,766.7436650901011,167.0029939851379
+2,replication_intermediate,0,2,258.80553693622494,163.04907236779974

tests/test_io.py

@@ -15,6 +15,7 @@
 from topostats.io import (
     LoadScans,
     convert_basename_to_relative_paths,
+    dict_almost_equal,
     dict_to_hdf5,
     find_files,
     get_date_time,
@@ -59,48 +60,6 @@
 # pylint: disable=too-many-lines
 
 
-def dict_almost_equal(dict1, dict2, abs_tol=1e-9):
-    """Recursively check if two dictionaries are almost equal with a given absolute tolerance.
-
-    Parameters
-    ----------
-    dict1: dict
-        First dictionary to compare.
-    dict2: dict
-        Second dictionary to compare.
-    abs_tol: float
-        Absolute tolerance to check for equality.
-
-    Returns
-    -------
-    bool
-        True if the dictionaries are almost equal, False otherwise.
-    """
-    if dict1.keys() != dict2.keys():
-        return False
-
-    LOGGER.info("Comparing dictionaries")
-
-    for key in dict1:
-        LOGGER.info(f"Comparing key {key}")
-        if isinstance(dict1[key], dict) and isinstance(dict2[key], dict):
-            if not dict_almost_equal(dict1[key], dict2[key], abs_tol=abs_tol):
-                return False
-        elif isinstance(dict1[key], np.ndarray) and isinstance(dict2[key], np.ndarray):
-            if not np.allclose(dict1[key], dict2[key], atol=abs_tol):
-                LOGGER.info(f"Key {key} type: {type(dict1[key])} not equal: {dict1[key]} != {dict2[key]}")
-                return False
-        elif isinstance(dict1[key], float) and isinstance(dict2[key], float):
-            if not np.isclose(dict1[key], dict2[key], atol=abs_tol):
-                LOGGER.info(f"Key {key} type: {type(dict1[key])} not equal: {dict1[key]} != {dict2[key]}")
-                return False
-        elif dict1[key] != dict2[key]:
-            LOGGER.info(f"Key {key} not equal: {dict1[key]} != {dict2[key]}")
-            return False
-
-    return True
-
-
 def test_get_date_time() -> None:
     """Test the fetching of a formatted date and time string."""
     assert datetime.strptime(get_date_time(), "%Y-%m-%d %H:%M:%S")
@@ -281,6 +240,13 @@ def test_load_array() -> None:
             False,
             id="float not equal",
         ),
+        pytest.param(
+            {"a": np.nan},
+            {"a": np.nan},
+            0.0001,
+            True,
+            id="nan equal",
+        ),
     ],
 )
 def test_dict_almost_equal(dict1: dict, dict2: dict, tolerance: float, expected: bool) -> None:
@@ -769,7 +735,7 @@ def test_dict_to_hdf5_all_together_group_path_non_standard(tmp_path: Path) -> No
         assert list(f.keys()) == list(expected.keys())
         assert f["d"]["a"][()] == expected["d"]["a"]
         np.testing.assert_array_equal(f["d"]["b"][()], expected["d"]["b"])
-        assert f["d"]["c"][()].decode("utf-8") == expected["d"]["c"]
+        assert f["d"]["c"][()].decode("utf-8") == expected["d"]["c"]  # pylint: disable=no-member
         assert f["d"]["d"]["e"][()] == expected["d"]["d"]["e"]
         np.testing.assert_array_equal(f["d"]["d"]["f"][()], expected["d"]["d"]["f"])
         assert f["d"]["d"]["g"][()].decode("utf-8") == expected["d"]["d"]["g"]

tests/tracing/conftest.py

@@ -1,6 +1,5 @@
 """Fixtures for the tracing tests."""
 
-import pickle
 from pathlib import Path
 
 import numpy as np
@@ -42,7 +41,7 @@
 @pytest.fixture()
 def test_dnatracing() -> dnaTrace:
     """Instantiate a dnaTrace object."""
-    return dnaTrace(image=FULL_IMAGE, grain=GRAINS, filename="Test", pixel_to_nm_scaling=1.0)
+    return dnaTrace(image=FULL_IMAGE, mask=GRAINS, filename="Test", pixel_to_nm_scaling=1.0)
 
 
 @pytest.fixture()
@@ -165,18 +164,6 @@ def catenane_image() -> npt.NDArray[np.number]:
     return np.load(RESOURCES / "catenane_image.npy")
 
 
-@pytest.fixture()
-def catenane_skeleton() -> npt.NDArray[np.bool_]:
-    """Skeleton of the catenane test image."""
-    return np.load(RESOURCES / "catenane_skeleton.npy")
-
-
-@pytest.fixture()
-def catenane_smoothed_mask() -> npt.NDArray[np.bool_]:
-    """Smoothed mask of the catenane test image."""
-    return np.load(RESOURCES / "catenane_smoothed_mask.npy")
-
-
 @pytest.fixture()
 def catenane_node_centre_mask() -> npt.NDArray[np.int32]:
     """
@@ -202,50 +189,30 @@ def catenane_connected_nodes() -> npt.NDArray[np.int32]:
 @pytest.fixture()
 def nodestats_catenane(
     catenane_image: npt.NDArray[np.number],
-    catenane_smoothed_mask: npt.NDArray[np.bool_],
-    catenane_skeleton: npt.NDArray[np.bool_],
-    catenane_node_centre_mask: npt.NDArray[np.int32],
-    catenane_connected_nodes: npt.NDArray[np.int32],
 ) -> nodeStats:
     """Fixture for the nodeStats object for a catenated molecule, to be used in analyse_nodes."""
+    catenane_smoothed_mask: npt.NDArray[np.bool_] = np.load(RESOURCES / "catenane_smoothed_mask.npy")
+    catenane_skeleton: npt.NDArray[np.bool_] = np.load(RESOURCES / "catenane_skeleton.npy")
+    catenane_node_centre_mask = np.load(RESOURCES / "catenane_node_centre_mask.npy")
+    catenane_connected_nodes = np.load(RESOURCES / "catenane_connected_nodes.npy")
+
     # Create a nodestats object
     nodestats = nodeStats(
         filename="test_catenane",
         image=catenane_image,
         mask=catenane_smoothed_mask,
         smoothed_mask=catenane_smoothed_mask,
         skeleton=catenane_skeleton,
-        px_2_nm=np.float64(0.18124609375),
+        pixel_to_nm_scaling=np.float64(0.18124609375),
         n_grain=1,
         node_joining_length=7,
         node_extend_dist=14.0,
         branch_pairing_length=20.0,
+        pair_odd_branches=True,
     )
 
     nodestats.node_centre_mask = catenane_node_centre_mask
     nodestats.connected_nodes = catenane_connected_nodes
     nodestats.skeleton = catenane_skeleton
 
     return nodestats
-
-
-# pylint: disable=unspecified-encoding
-@pytest.fixture()
-def nodestats_catenane_node_dict() -> dict:
-    """Node dictionary for the catenane test image."""
-    with Path.open(RESOURCES / "catenane_node_dict.pkl", "rb") as file:
-        return pickle.load(file)  # noqa: S301 - Pickles unsafe but we don't care
-
-
-# pylint: disable=unspecified-encoding
-@pytest.fixture()
-def nodestats_catenane_image_dict() -> dict:
-    """Image dictionary for the catenane test image."""
-    with Path.open(RESOURCES / "catenane_image_dict.pkl", "rb") as file:
-        return pickle.load(file)  # noqa: S301 - Pickles unsafe but we don't care
-
-
-@pytest.fixture()
-def nodestats_catenane_all_connected_nodes() -> npt.NDArray[np.int32]:
-    """All connected nodes for the catenane test image."""
-    return np.load(RESOURCES / "catenane_all_connected_nodes.npy")