extended moments features

CHANGELOG

@@ -1,3 +1,6 @@
+0.12.2
+ - enh: feature computation for area_um_raw, deform_raw, per_ratio,
+   and per_um_raw
 0.12.1
  - fix: avoid NaN-valued features due to invalid contours (#9)
 0.12.0

dcnum/feat/feat_moments/mt_legacy.py

@@ -11,45 +11,61 @@ def moments_based_features(mask, pixel_size):
     size = mask.shape[0]
 
     empty = np.full(size, np.nan, dtype=np.float64)
-    deform = np.copy(empty)
-    size_x = np.copy(empty)
-    size_y = np.copy(empty)
-    pos_x = np.copy(empty)
-    pos_y = np.copy(empty)
+
+    # features from raw contour
     area_msd = np.copy(empty)
     area_ratio = np.copy(empty)
-    area_um = np.copy(empty)
+    area_um_raw = np.copy(empty)
     aspect = np.copy(empty)
+    deform_raw = np.copy(empty)
+    inert_ratio_prnc = np.copy(empty)
+    inert_ratio_raw = np.copy(empty)
+    per_ratio = np.copy(empty)
+    per_um_raw = np.copy(empty)
+    size_x = np.copy(empty)
+    size_y = np.copy(empty)
     tilt = np.copy(empty)
+
+    # features from convex hull
+    area_um = np.copy(empty)
+    deform = np.copy(empty)
     inert_ratio_cvx = np.copy(empty)
-    inert_ratio_raw = np.copy(empty)
-    inert_ratio_prnc = np.copy(empty)
+    pos_x = np.copy(empty)
+    pos_y = np.copy(empty)
+
     # The following valid-array is not a real feature, but only
     # used to figure out which events need to be removed due
     # to invalid computed features, often due to invalid contours.
     valid = np.full(size, False)
 
     for ii in range(size):
+        # raw contour
         cont_raw = contour_single_opencv(mask[ii])
         if len(cont_raw.shape) < 2:
             continue
         if cv2.contourArea(cont_raw) == 0:
             continue
+
         mu_raw = cv2.moments(cont_raw)
+        arc_raw = np.float64(cv2.arcLength(cont_raw, True))
+        area_raw = np.float64(mu_raw["m00"])
 
         # convex hull
         cont_cvx = np.squeeze(cv2.convexHull(cont_raw))
         mu_cvx = cv2.moments(cont_cvx)
+        arc_cvx = np.float64(cv2.arcLength(cont_cvx, True))
+        # circ
+        circ = 2 * np.sqrt(np.pi * mu_cvx["m00"]) / arc_cvx
+        deform[ii] = 1 - circ
 
         if mu_cvx["m00"] == 0 or mu_raw["m00"] == 0:
-            # Contour size too small
+            # contour size too small
             continue
 
-        arc = cv2.arcLength(cont_cvx, True)
-
+        # bounding box
         x, y, w, h = cv2.boundingRect(cont_raw)
 
-        # tilt
+        # tilt angle in radians
         oii = 0.5 * np.arctan2(2 * mu_raw['mu11'],
                                mu_raw['mu02'] - mu_raw['mu20'])
         # +PI/2 because relative to channel axis
@@ -60,18 +76,18 @@ def moments_based_features(mask, pixel_size):
             tilti -= np.pi
         tilt[ii] = np.abs(tilti)
 
-        # circ
-        circ = 2 * np.sqrt(np.pi * mu_cvx["m00"]) / arc
-        deform[ii] = 1 - circ
-
         size_x[ii] = w * pixel_size
         size_y[ii] = h * pixel_size
         pos_x[ii] = mu_cvx["m10"] / mu_cvx["m00"] * pixel_size
         pos_y[ii] = mu_cvx["m01"] / mu_cvx["m00"] * pixel_size
         area_msd[ii] = mu_raw["m00"]
+        area_um_raw[ii] = area_raw * pixel_size**2
         area_ratio[ii] = mu_cvx["m00"] / mu_raw["m00"]
-        area_um[ii] = mu_cvx["m00"] * pixel_size ** 2
+        area_um[ii] = mu_cvx["m00"] * pixel_size**2
         aspect[ii] = w / h
+        per_um_raw[ii] = arc_raw * pixel_size
+        deform_raw[ii] = 1 - 2 * np.sqrt(np.pi * area_raw) / arc_raw
+        per_ratio[ii] = arc_raw / arc_cvx
 
         # inert_ratio_cvx
         if mu_cvx['mu02'] > 0:  # defaults to zero
@@ -95,6 +111,8 @@ def moments_based_features(mask, pixel_size):
         root_prnc = mprnc["mu20"] / mprnc["mu02"]
         if root_prnc > 0:  # defaults to zero
             inert_ratio_prnc[ii] = np.sqrt(root_prnc)
+
+        # specify validity
         valid[ii] = True
 
     return {
@@ -111,5 +129,9 @@ def moments_based_features(mask, pixel_size):
         "inert_ratio_cvx": inert_ratio_cvx,
         "inert_ratio_raw": inert_ratio_raw,
         "inert_ratio_prnc": inert_ratio_prnc,
+        "area_um_raw": area_um_raw,
+        "per_um_raw": per_um_raw,
+        "deform_raw": deform_raw,
+        "per_ratio": per_ratio,
         "valid": valid,
     }

tests/test_feat_moments_based_extended.py

@@ -0,0 +1,58 @@
+import pathlib
+
+import h5py
+import numpy as np
+
+from dcnum.feat import feat_moments
+
+from helper_methods import retrieve_data
+
+data_path = pathlib.Path(__file__).parent / "data"
+
+
+def test_moments_based_features():
+    # This file has new cell features belonging to
+    # fmt-hdf5_cytoshot_full-features_2023.zip
+    path = retrieve_data(data_path /
+                         "fmt-hdf5_cytoshot_extended-moments-features.zip")
+
+    feats = [
+            "area_um_raw",
+            "per_um_raw",
+            "deform_raw",
+            "per_ratio",
+    ]
+
+    # Make data available
+    with h5py.File(path) as h5:
+        data = feat_moments.moments_based_features(
+            mask=h5["events/mask"][:],
+            pixel_size=0.2645
+        )
+        for feat in feats:
+            rtol = 1e-5
+            atol = 1e-8
+            assert np.allclose(h5["events"][feat][:],
+                               data[feat],
+                               rtol=rtol,
+                               atol=atol), f"Feature {feat} mismatch!"
+
+
+def test_mask_2d():
+    masks = np.array([
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 1, 1, 0, 0],
+        [0, 0, 1, 1, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0],
+    ], dtype=bool)[np.newaxis]
+    data = feat_moments.moments_based_features(
+                mask=masks,
+                pixel_size=0.2645
+            )
+    assert data["deform_raw"].shape == (1,)
+    # This is the deformation of a square (compared to circle)
+    assert np.allclose(data["deform_raw"][0], 0.11377307454724206)
+    # Without moments-based computation, this would be 4*pxsize=0.066125
+    assert np.allclose(data["area_um_raw"][0], 0.06996025)