Merge pull request #853 from ANTsX/labelstats

BUG: label_geometry_measures with uint32 labels

Author: cookpa

ants/label/label_geometry_measures.py

@@ -43,6 +43,8 @@ def label_geometry_measures(label_image, intensity_image=None):
         if not np.all(label_image.numpy() == label_image_int.numpy()):
             raise ValueError('Input label values must be representable as uint32.')
         label_image = label_image.clone('unsigned int')
+    else:
+        label_image_int = label_image
     veccer = [label_image.dimension, label_image_int, intensity_image, outcsv]
     veccer_processed = process_arguments(veccer)
     libfn = get_lib_fn('LabelGeometryMeasures')

ants/label/label_stats.py

@@ -1,5 +1,6 @@
 __all__ = ["label_stats"]
 
+import numpy as np
 import pandas as pd
 
 from ants.internal import get_lib_fn
@@ -8,7 +9,7 @@
 @image_method
 def label_stats(image, label_image):
     """
-    Get label statistics from image
+    Get label statistics from an image. The labels must be representable as uint32.
 
     ANTsR function: `labelStats`
 
@@ -22,7 +23,7 @@ def label_stats(image, label_image):
 
     Returns
     -------
-    ndarray ?
+    pandas.DataFrame
 
     Example
     -------
@@ -34,10 +35,16 @@ def label_stats(image, label_image):
     >>> stats = ants.label_stats(image, segs1['segmentation'])
     """
     image_float = image.clone("float")
-    label_image_int = label_image.clone("unsigned int")
+
+    label_image_int = label_image.clone('unsigned int')
+
+    if label_image.pixeltype != 'unsigned int':
+        if not np.all(label_image.numpy() == label_image_int.numpy()):
+            raise ValueError('Input label values must be representable as uint32.')
 
     libfn = get_lib_fn("labelStats%iD" % image.dimension)
     df = libfn(image_float.pointer, label_image_int.pointer)
     df = pd.DataFrame(df)
     df.sort_values(by=["LabelValue"], inplace=True)
+    df = df.reset_index(drop=True)
     return df

src/labelStats.cxx

@@ -40,12 +40,12 @@ nb::dict labelStatsHelper(
 
   long nlabs = labelStatisticsImageFilter->GetNumberOfLabels();
 
-  std::vector<double> labelvals(nlabs);
+  std::vector<LabelType> labelvals(nlabs);
   std::vector<double> means(nlabs);
   std::vector<double> mins(nlabs);
   std::vector<double> maxes(nlabs);
   std::vector<double> variances(nlabs);
-  std::vector<double> counts(nlabs);
+  std::vector<LabelType> counts(nlabs);
   std::vector<double> volumes(nlabs);
   std::vector<double> x(nlabs);
   std::vector<double> y(nlabs);
@@ -137,7 +137,7 @@ nb::dict labelStatsHelper(
 template <unsigned int Dimension>
 nb::dict labelStats(AntsImage<itk::Image<float, Dimension>> & py_image,
                     AntsImage<itk::Image<unsigned int, Dimension>> & py_labelImage)
-{ 
+{
   typedef itk::Image<float, Dimension> FloatImageType;
   typedef itk::Image<unsigned int, Dimension> IntImageType;
   typedef typename FloatImageType::Pointer FloatImagePointerType;

tests/test_utils.py

@@ -996,7 +996,7 @@ def test_channels_first(self):
         img4 = ants.merge_channels([image,image2], channels_first=True)
         self.assertTrue(np.allclose(img3.numpy()[:,:,0], img4.numpy()[0,:,:]))
         self.assertTrue(np.allclose(img3.numpy()[:,:,1], img4.numpy()[1,:,:]))
-        
+
     def test_polar_decomposition(self):
         # Helper functions for creating known matrices
         def make_known_rotation(theta_deg):
@@ -1006,16 +1006,16 @@ def make_known_rotation(theta_deg):
                         [np.sin(theta),  np.cos(theta)]]
             return R
         def make_known_scaling_matrix(sx, sy, sz):
-            return np.diag([sx, sy, sz])        
+            return np.diag([sx, sy, sz])
 
-        # 1. Setup: Create a matrix from a known P and Z.        
+        # 1. Setup: Create a matrix from a known P and Z.
         # The key is to multiply them in the order P @ Z.
         P_known = make_known_scaling_matrix(2.5, 1.0, 1.5)
-        Z_known = make_known_rotation(42) # Use Z for "orthogonal" part        
+        Z_known = make_known_rotation(42) # Use Z for "orthogonal" part
         # Construct X using the left decomposition structure: X = P @ Z
         X = P_known @ Z_known
 
-        # 2. Decompose the matrix using the function       
+        # 2. Decompose the matrix using the function
         result = ants.polar_decomposition(X)
         P_est = result["P"]
         Z_est = result["Z"]
@@ -1028,7 +1028,7 @@ def make_known_scaling_matrix(sx, sy, sz):
         nptest.assert_almost_equal(P_known, P_est)
         # c. Check if the estimated orthogonal part is close to the known one.
         nptest.assert_almost_equal(Z_known, Z_est)
-        
+
     def test_convergence_monitoring(self):
         f = [1 / i for i in range(1, 100)]
         convergence = ants.convergence_monitoring(f, window_size=10)