Shape Features

requirements.txt

@@ -3,11 +3,12 @@ matplotlib~=3.7.0
 nibabel~=5.0.1
 notebook~=6.5.1
 numpy~=1.24.1
-pandas~=1.5.1
+pandas~=2.0.2
 pathos~=0.3.0
-renalsegmentor~=1.3.6
-scikit-image~=0.19.1
-scipy~=1.10.0
+renalsegmentor~=1.3.7
+scikit-image~=0.20.0
+scipy~=1.9.1
 scikit-learn~=1.2.1
 tabulate~=0.9.0
 tqdm~=4.64.1
+trimesh~=3.22.0

ukat/segmentation/__init__.py

@@ -1 +1 @@
-from . import whole_kidney
+from . import shape_features, whole_kidney

ukat/segmentation/shape_features.py

@@ -0,0 +1,161 @@
+import numpy as np
+import pandas as pd
+import trimesh
+
+from nibabel.affines import voxel_sizes
+from skimage.measure import label, marching_cubes, regionprops
+from trimesh import smoothing
+
+
+class ShapeFeatures:
+    def __init__(self, pixel_array, affine=None, zoom=None, kidneys=True,
+                 region_labels=None):
+        """
+        Calculate shape features for a mask.
+
+        Parameters
+        ----------
+        pixel_array : np.ndarray(dtype=np.uint8)
+            An array containing the mask, this array can either be a binary
+            array of 0s and 1s where 0 represents background tissue or an
+            array of integers where each integer represents a different tissue.
+            If a binary mask is provided, the features of each isolated
+            region will be calculated.
+        affine : np.ndarray
+            A matrix giving the relationship between voxel coordinates and
+            world coordinates.
+        zoom : tuple of float, shape (ndim,)
+            A tuple of floats giving the voxel size in mm.
+        kidneys : bool, optional
+            Default True
+            If true, it will be assumed that the two regions in the mask are
+            the left and right kidneys.
+        region_labels : list, optional
+            A list of string labels corresponding to each integer label in
+            pixel_array
+        """
+        if (affine is None) and (zoom is None):
+            raise ValueError('Affine or zoom must be specified')
+
+        self.pixel_array = pixel_array
+        self.affine = affine
+        if zoom is None:
+            zoom = tuple(voxel_sizes(self.affine))
+        self.zoom = zoom
+        self.n_labels = len(np.unique(pixel_array[pixel_array > 0]))
+        if self.n_labels == 1:
+            self.labels = label(self.pixel_array)
+            self.n_labels = len(np.unique(self.labels[self.labels > 0]))
+        else:
+            self.labels = self.pixel_array
+
+        if kidneys:
+            if self.n_labels != 2:
+                raise ValueError('Expected two labels (L and R) if '
+                                 'kidney=True')
+            if region_labels is not None:
+                self.region_labels = region_labels
+            else:
+                self.region_labels = ['L', 'R']
+        elif region_labels is not None:
+            if len(region_labels) != self.n_labels:
+                raise ValueError('The number of labels must match the number '
+                                 'of regions in the mask')
+            self.region_labels = region_labels
+        else:
+            self.region_labels = np.arange(1, self.n_labels + 1)
+
+        self.features = self.get_features()
+
+    def get_features(self):
+        """
+        Calculate shape features for each region in the mask.
+
+        Returns
+        -------
+        props_df : pd.DataFrame
+            A dataframe containing the calculated shape features for each
+            region in the mask.
+        """
+        properties = ['volume', 'surface_area', 'volume_bbox', 'volume_convex',
+                      'volume_filled', 'n_vox', 'long_axis',
+                      'short_axis', 'compactness', 'euler_number', 'solidity']
+        props_df = pd.DataFrame(index=self.region_labels, columns=properties)
+
+        for region, label in zip(self.region_labels,
+                                 np.unique(self.labels[self.labels > 0])):
+            props_df.loc[region] = self._get_region_props(self.labels == label)
+        return props_df
+
+    def save_features_csv(self, path):
+        """
+        Save the calculated shape features to a csv file.
+
+        Parameters
+        ----------
+        path : str
+            The path to save the csv file to.
+        """
+        self.features.to_csv(path)
+
+    def _get_region_props(self, region):
+        """
+        Calculate shape features for a single region.
+
+        Parameters
+        ----------
+        region : np.ndarray(dtype=np.bool)
+            A binary array of 0s and 1s where 0 represents background tissue.
+
+        Returns
+        -------
+        props_dict : dict
+            A dictionary containing the calculated shape features for the
+            region.
+        """
+        props = regionprops(region.astype(np.uint8), spacing=self.zoom)[0]
+        mesh = self._get_smoothed_mesh(region)
+
+        props_dict = {}
+        props_dict.update({'volume': props['area'] / 1000})  # mm^3 to mL
+        props_dict.update({'surface_area': mesh.area / 100})  # mm^2 to cm^2
+        props_dict.update(
+            {'volume_bbox': props['bbox_area'] / 1000})  # mm^3 to mL
+        props_dict.update(
+            {'volume_convex': props['convex_area'] / 1000})  # mm^3 to mL
+        props_dict.update(
+            {'volume_filled': props['filled_area'] / 1000})  # mm^3 to mL
+        props_dict.update({'n_vox': props['num_pixels']})
+        props_dict.update(
+            {'long_axis': props['major_axis_length'] / 10})  # mm to cm
+        props_dict.update(
+            {'short_axis': props['minor_axis_length'] / 10})  # mm to cm
+        compactness = (props_dict['volume'] / props_dict['surface_area']) / \
+                      (props['equivalent_diameter_area'] / 6)
+        props_dict.update({'compactness': compactness})
+        props_dict.update({'euler_number': props['euler_number']})
+        props_dict.update({'solidity': props['solidity']})
+
+        return props_dict
+
+    def _get_smoothed_mesh(self, region):
+        """
+        Generate a smoothed mesh from a binary region. Parameters have been
+        optimised for kidneys.
+
+        Parameters
+        ----------
+        region : np.ndarray(dtype=np.bool)
+            A binary array of 0s and 1s where 0 represents background tissue.
+
+        Returns
+        -------
+        mesh : trimesh.Trimesh
+            A smoothed mesh representation of the region.
+        """
+        verts, faces, _, _ = marching_cubes(region.astype(np.uint8),
+                                            spacing=self.zoom, level=0.5,
+                                            step_size=1.0)
+        mesh = trimesh.Trimesh(vertices=verts, faces=faces)
+        mesh = smoothing.filter_laplacian(mesh, lamb=1, iterations=20)
+        return mesh

ukat/segmentation/tests/__init__.py

@@ -1 +1 @@
-from . import test_whole_kidney
+from . import test_shape_features, test_whole_kidney

ukat/segmentation/tests/test_shape_features.py

@@ -0,0 +1,189 @@
+import csv
+import os
+import shutil
+
+import numpy as np
+import numpy.testing as npt
+import pandas as pd
+import pytest
+
+from ukat.data import fetch
+from ukat.segmentation.shape_features import ShapeFeatures
+from ukat.segmentation.whole_kidney import Segmentation
+from ukat.utils import arraystats
+
+
+class TestShapeFeatures:
+    image, affine = fetch.t2w_volume_philips()
+    zoom = (1.5000001667213594, 1.5000000306853905, 5.49999480801915)
+    segmentation = Segmentation(image, affine)
+    mask = segmentation.get_mask()
+    kidneys = segmentation.get_kidneys()
+
+    mask_two_body = np.zeros((50, 50, 10))
+    mask_two_body[5:15, 5:45, :] = 1
+    mask_two_body[20:30, 5:45, :] = 1
+
+    mask_three_body = np.zeros((50, 50, 10))
+    mask_three_body[5:15, 5:45, :] = 1
+    mask_three_body[20:30, 5:45, :] = 1
+    mask_three_body[35:45, 5:45, :] = 1
+
+    simulated_affine = np.eye(4)
+
+    def test_get_smoothed_mesh(self):
+        expected = [175.680058, 104.798285, 3.267293, 291.491921]
+        shape_features = ShapeFeatures(self.kidneys, self.affine)
+        mesh = shape_features._get_smoothed_mesh(self.kidneys == 1)
+
+        assert mesh.is_watertight
+        vertex_stats = arraystats.ArrayStats(mesh.vertices).calculate()
+        npt.assert_allclose([vertex_stats["mean"], vertex_stats["std"],
+                             vertex_stats["min"], vertex_stats["max"]],
+                            expected, rtol=1e-6, atol=1e-4)
+
+    def test_get_region_props(self):
+        shape_features = ShapeFeatures(self.kidneys, self.affine)
+        props_dict = shape_features._get_region_props(self.kidneys == 1)
+        assert props_dict == pytest.approx({'volume': 118.19352898042803,
+                                            'surface_area': 148.05689835989392,
+                                            'volume_bbox': 360.8547068442343,
+                                            'volume_convex': 170.52736146479933,
+                                            'volume_filled': 118.19352898042803,
+                                            'n_vox': 9551,
+                                            'long_axis': 11.793750181329315,
+                                            'short_axis': 4.347012606736469,
+                                            'compactness': 0.07866555492167773,
+                                            'euler_number': 2,
+                                            'solidity': 0.6931059506531204},
+                                           rel=1e-20, abs=1e-4)
+
+    def test_shape_features_labels_affine(self):
+        shape_features = ShapeFeatures(self.kidneys, self.affine)
+        features_df = shape_features.get_features()
+        gold_df = pd.DataFrame(index=['L', 'R'],
+                               columns=['volume', 'surface_area',
+                                        'volume_bbox', 'volume_convex',
+                                        'volume_filled', 'n_vox',
+                                        'long_axis', 'short_axis',
+                                        'compactness', 'euler_number',
+                                        'solidity'],
+                               data=[[118.19352898042803, 148.05689835989392,
+                                      360.8547068442343,
+                                      170.52736146479933, 118.19352898042803,
+                                      9551.0,
+                                      11.793750181329315, 4.347012606736469,
+                                      0.07866555492167773, 2.0,
+                                      0.6931059506531204],
+                                     [121.80702604484902, 154.5164344861936,
+                                      263.7357857429465,
+                                      150.36850284171092, 121.80702604484902,
+                                      9843.0,
+                                      12.317681104489647, 3.6430077535636998,
+                                      0.0769055483268716, 2.0,
+                                      0.8100567854497571]])
+        pd.testing.assert_frame_equal(features_df, gold_df, check_dtype=False)
+
+    def test_shape_features_labels_zoom(self):
+        shape_features = ShapeFeatures(self.kidneys, zoom=self.zoom)
+        features_df = shape_features.get_features()
+        gold_df = pd.DataFrame(index=['L', 'R'],
+                               columns=['volume', 'surface_area',
+                                        'volume_bbox', 'volume_convex',
+                                        'volume_filled', 'n_vox',
+                                        'long_axis', 'short_axis',
+                                        'compactness', 'euler_number',
+                                        'solidity'],
+                               data=[[118.19352898042803, 148.05689835989392,
+                                      360.8547068442343,
+                                      170.52736146479933, 118.19352898042803,
+                                      9551.0,
+                                      11.793750181329315, 4.347012606736469,
+                                      0.07866555492167773, 2.0,
+                                      0.6931059506531204],
+                                     [121.80702604484902, 154.5164344861936,
+                                      263.7357857429465,
+                                      150.36850284171092, 121.80702604484902,
+                                      9843.0,
+                                      12.317681104489647, 3.6430077535636998,
+                                      0.0769055483268716, 2.0,
+                                      0.8100567854497571]])
+        pd.testing.assert_frame_equal(features_df, gold_df, check_dtype=False)
+
+    def test_shape_features_mask(self):
+        shape_features = ShapeFeatures(self.mask_two_body,
+                                       self.simulated_affine)
+        features_df = shape_features.get_features()
+        gold_df = pd.DataFrame(index=['L', 'R'],
+                               columns=['volume', 'surface_area',
+                                        'volume_bbox', 'volume_convex',
+                                        'volume_filled', 'n_vox',
+                                        'long_axis', 'short_axis',
+                                        'compactness', 'euler_number',
+                                        'solidity'],
+                               data=[[4.0, 6.354534260930347, 4.0, 4.0, 4.0,
+                                      4000.0, 5.162363799656123,
+                                      1.284523257866513, 0.1917669347647048,
+                                      1.0, 1.0],
+                                     [4.0, 6.413147961900247, 4.0, 4.0, 4.0,
+                                      4000.0, 5.162363799656123,
+                                      1.284523257866513, 0.1900142588811934,
+                                      1.0, 1.0]])
+        pd.testing.assert_frame_equal(features_df, gold_df, check_dtype=False)
+
+    def test_no_zoom_no_affine(self):
+        with pytest.raises(ValueError):
+            ShapeFeatures(self.kidneys)
+
+    def test_region_labels_length_doesnt_match(self):
+
+        with pytest.raises(ValueError):
+            ShapeFeatures(self.mask_three_body, self.simulated_affine,
+                          region_labels=['L', 'R'],
+                          kidneys=False)
+
+    def test_custom_labels(self):
+        shape_features = ShapeFeatures(self.mask, self.affine,
+                                       region_labels=['L', 'R', 'O'],
+                                       kidneys=False)
+        assert np.all(shape_features.get_features().index == ['L', 'R', 'O'])
+
+    def test_no_labels_not_kidneys(self):
+        shape_features = ShapeFeatures(self.mask_three_body,
+                                       self.simulated_affine,
+                                       kidneys=False)
+        assert np.all(shape_features.get_features().index == [1, 2, 3])
+
+    def test_three_regions_kidney(self):
+        with pytest.raises(ValueError):
+            ShapeFeatures(self.mask_three_body, self.simulated_affine,
+                          kidneys=True)
+
+    def test_save_csv(self):
+        shape_features = ShapeFeatures(self.kidneys, self.affine)
+        expected_cols = ['', 'volume', 'surface_area', 'volume_bbox',
+                         'volume_convex', 'volume_filled', 'n_vox',
+                         'long_axis', 'short_axis', 'compactness',
+                         'euler_number', 'solidity']
+
+        if os.path.exists('test_output'):
+            shutil.rmtree('test_output')
+        os.makedirs('test_output', exist_ok=True)
+
+        shape_features.save_features_csv('test_output/features.csv')
+        output_files = os.listdir('test_output')
+
+        assert 'features.csv' in output_files
+
+        with open('test_output/features.csv', 'r') as csv_file:
+            reader = csv.reader(csv_file)
+            list_from_csv = [row for row in reader]
+        assert list_from_csv[0] == expected_cols
+        assert len(list_from_csv) == 3
+        assert len(list_from_csv[1]) == 12
+        assert list_from_csv[1][0] == 'L'
+        assert list_from_csv[2][0] == 'R'
+
+        for f in os.listdir('test_output'):
+            os.remove(os.path.join('test_output', f))
+        shutil.rmtree('test_output')