Update convert model units to infer model resolution

src/arcade_collection/output/convert_model_units.py

@@ -1,27 +1,76 @@
+import re
 from typing import Optional, Union
 
 import pandas as pd
 
 
 def convert_model_units(
-    data: pd.DataFrame, ds: float, dt: float, regions: Optional[Union[list[str], str]] = None
+    data: pd.DataFrame,
+    ds: Optional[float],
+    dt: Optional[float],
+    regions: Optional[Union[list[str], str]] = None,
 ) -> None:
-    data["time"] = round(dt * data["TICK"], 2)
+    """
+    Converts data from simulation units to true units.
 
-    if "NUM_VOXELS" in data.columns:
-        data["volume"] = ds * ds * ds * data["NUM_VOXELS"]
+    Simulations use spatial unit of voxels and temporal unit of ticks. Spatial
+    resolution (microns/voxel) and temporal resolution (hours/tick) are used to
+    convert data to true units. If spatial or temporal resolution is not given,
+    they will be estimated from the ``KEY`` column of the data.
 
-    if "MAX_Z" in data.columns and "MIN_Z" in data.columns:
-        data["height"] = ds * (data["MAX_Z"] - data["MIN_Z"] + 1)
+    The following columns are added to the data:
 
-    if "CX" in data.columns:
-        data["cx"] = ds * data["CX"]
+    =============  ===================  =============================
+    Target column  Source column(s)      Conversion
+    =============  ===================  =============================
+    ``time``       ``TICK``             ``dt * TICK``
+    ``volume``     ``NUM_VOXELS``       ``ds * ds * ds * NUM_VOXELS``
+    ``height``     ``MAX_Z`` ``MIN_Z``  ``ds * (MAX_Z - MIN_Z + 1)``
+    ``cx``         ``CENTER_X``         ``ds * CENTER_X``
+    ``cy``         ``CENTER_Y``         ``ds * CENTER_Y``
+    ``cz``         ``CENTER_Z``         ``ds * CENTER_Z``
+    =============  ===================  =============================
 
-    if "CY" in data.columns:
-        data["cy"] = ds * data["CY"]
+    For each region (other than ``DEFAULT``), the following columns are added to the data:
 
-    if "CZ" in data.columns:
-        data["cz"] = ds * data["CZ"]
+    =================  =================================  ==========================================
+    Target column      Source column(s)                   Conversion
+    =================  =================================  ==========================================
+    ``volume.REGION``  ``NUM_VOXELS.REGION``              ``ds * ds * ds * NUM_VOXELS.REGION``
+    ``height.REGION``  ``MAX_Z.REGION`` ``MIN_Z.REGION``  ``ds * (MAX_Z.REGION - MIN_Z.REGION + 1)``
+    =================  =================================  ==========================================
+
+    The following property columns are rescaled:
+
+    =====================  =====================  ==========================
+    Target column          Source column(s)        Conversion
+    =====================  =====================  ==========================
+    ``area``               ``area``               ``ds * ds * area``
+    ``perimeter``          ``perimeter``          ``ds * perimeter``
+    ``axis_major_length``  ``axis_major_length``  ``ds * axis_major_length``
+    ``axis_minor_length``  ``axis_minor_length``  ``ds * axis_minor_length``
+    =====================  =====================  ==========================
+
+    Parameters
+    ----------
+    data
+        Simulation data.
+    ds
+        Spatial resolution in microns/voxel, use None to estimate from keys.
+    dt
+        Temporal resolution in hours/tick, use None to estimate from keys.
+    regions
+        List of regions.
+    """
+
+    if dt is None:
+        dt = data["KEY"].apply(estimate_temporal_resolution)
+
+    if ds is None:
+        ds = data["KEY"].apply(estimate_spatial_resolution)
+
+    convert_temporal_units(data, dt)
+    convert_spatial_units(data, ds)
 
     if regions is None:
         return
@@ -33,5 +82,146 @@ def convert_model_units(
         if region == "DEFAULT":
             continue
 
-        data[f"volume.{region}"] = ds * ds * ds * data[f"NUM_VOXELS.{region}"]
-        data[f"height.{region}"] = ds * (data[f"MAX_Z.{region}"] - data[f"MIN_Z.{region}"])
+        convert_spatial_units(data, ds, region)
+
+
+def convert_temporal_units(data: pd.DataFrame, dt: float) -> None:
+    """
+    Converts temporal data from simulation units to true units.
+
+    Simulations use temporal unit of ticks. Temporal resolution (hours/tick) is
+    used to convert data to true units.
+
+    The following temporal columns are converted:
+
+    =============  ================  =============
+    Target column  Source column(s)  Conversion
+    =============  ================  =============
+    ``time``       ``TICK``          ``dt * TICK``
+    =============  ================  =============
+
+    Parameters
+    ----------
+    data
+        Simulation data.
+    dt
+        Temporal resolution in hours/tick.
+    """
+
+    if "TICK" in data.columns:
+        data["time"] = round(dt * data["TICK"], 2)
+
+
+def convert_spatial_units(data: pd.DataFrame, ds: float, region: Optional[str] = None) -> None:
+    """
+    Converts spatial data from simulation units to true units.
+
+    Simulations use spatial unit of voxels. Spatial resolution (microns/voxel)
+    is used to convert data to true units.
+
+    The following spatial columns are converted:
+
+    =====================  =====================  =============================
+    Target column          Source column(s)        Conversion
+    =====================  =====================  =============================
+    ``volume``             ``NUM_VOXELS``         ``ds * ds * ds * NUM_VOXELS``
+    ``height``             ``MAX_Z`` ``MIN_Z``    ``ds * (MAX_Z - MIN_Z + 1)``
+    ``cx``                 ``CENTER_X``           ``ds * CENTER_X``
+    ``cy``                 ``CENTER_Y``           ``ds * CENTER_Y``
+    ``cz``                 ``CENTER_Z``           ``ds * CENTER_Z``
+    ``area``               ``area``               ``ds * ds * area``
+    ``perimeter``          ``perimeter``          ``ds * perimeter``
+    ``axis_major_length``  ``axis_major_length``  ``ds * axis_major_length``
+    ``axis_minor_length``  ``axis_minor_length``  ``ds * axis_minor_length``
+    =====================  =====================  =============================
+
+    Note that the centroid columns (``cx``, ``cy``, and ``cz``) are only
+    converted for the entire cell (``region == None``).
+
+    Parameters
+    ----------
+    data
+        Simulation data.
+    ds
+        Spatial resolution in microns/voxel.
+    region
+        Name of region.
+    """
+
+    suffix = "" if region is None else f".{region}"
+
+    if f"NUM_VOXELS{suffix}" in data.columns:
+        data[f"volume{suffix}"] = ds * ds * ds * data[f"NUM_VOXELS{suffix}"]
+
+    if f"MAX_Z{suffix}" in data.columns and f"MIN_Z{suffix}" in data.columns:
+        data[f"height{suffix}"] = ds * (data[f"MAX_Z{suffix}"] - data[f"MIN_Z{suffix}"] + 1)
+
+    if "CENTER_X" in data.columns and region is None:
+        data["cx"] = ds * data["CENTER_X"]
+
+    if "CENTER_Y" in data.columns and region is None:
+        data["cy"] = ds * data["CENTER_Y"]
+
+    if "CENTER_Z" in data.columns and region is None:
+        data["cz"] = ds * data["CENTER_Z"]
+
+    property_conversions = [
+        ("area", ds * ds),
+        ("perimeter", ds),
+        ("axis_major_length", ds),
+        ("axis_minor_length", ds),
+    ]
+
+    for name, conversion in property_conversions:
+        column = f"{name}{suffix}"
+
+        if column not in data.columns:
+            continue
+
+        data[column] = data[column] * conversion
+
+
+def estimate_temporal_resolution(key: str) -> float:
+    """
+    Estimates temporal resolution based on condition key.
+
+    If the key contains ``DT##``, where ``##`` denotes the temporal resolution
+    in minutes/tick, temporal resolution is estimated from ``##``. Otherwise,
+    the default temporal resolution is 1 hours/tick.
+
+    Parameters
+    ----------
+    key
+        Condition key.
+
+    Returns
+    -------
+    :
+        Temporal resolution (hours/tick).
+    """
+
+    matches = [re.fullmatch(r"DT([0-9]+)", k) for k in key.split("_")]
+    return next((float(match.group(1)) / 60 for match in matches if match is not None), 1.0)
+
+
+def estimate_spatial_resolution(key: str) -> float:
+    """
+    Estimates spatial resolution based on condition key.
+
+    If the key contains ``DS##``, where ``##`` denotes the spatial resolution
+    in micron/voxel, spatial resolution is estimated from ``##``. Otherwise,
+    the default spatial resolution is 1 micron/voxel.
+
+    Parameters
+    ----------
+    key
+        Condition key.
+
+    Returns
+    -------
+    :
+        Spatial resolution (micron/voxel).
+    """
+
+    matches = [re.fullmatch(r"DS([0-9]+)", k) for k in key.split("_")]
+    return next((float(match.group(1)) for match in matches if match is not None), 1.0)

tests/arcade_collection/output/test_convert_model_units.py

@@ -0,0 +1,53 @@
+import random
+import unittest
+
+from arcade_collection.output.convert_model_units import (
+    estimate_spatial_resolution,
+    estimate_temporal_resolution,
+)
+
+
+class TestExtractVoxelContours(unittest.TestCase):
+    def test_estimate_temporal_resolution_missing_temporal_key(self):
+        self.assertEqual(1, estimate_temporal_resolution(""))
+        self.assertEqual(1, estimate_temporal_resolution("A"))
+        self.assertEqual(1, estimate_temporal_resolution("A_B"))
+        self.assertEqual(1, estimate_temporal_resolution("A_B_C"))
+
+    def test_estimate_temporal_resolution_valid_temporal_key(self):
+        dt = int(random.random() * 10)
+        dt_key = f"DT{dt:03d}"
+
+        self.assertEqual(dt / 60, estimate_temporal_resolution(f"{dt_key}_B_C"))
+        self.assertEqual(dt / 60, estimate_temporal_resolution(f"A_{dt_key}_C"))
+        self.assertEqual(dt / 60, estimate_temporal_resolution(f"A_B_{dt_key}"))
+
+    def test_estimate_temporal_resolution_invalid_temporal_key(self):
+        dt = int(random.random() * 10)
+        dt_key = f"DT{dt:03d}x"
+
+        self.assertEqual(1, estimate_temporal_resolution(f"{dt_key}_B_C"))
+        self.assertEqual(1, estimate_temporal_resolution(f"A_{dt_key}_C"))
+        self.assertEqual(1, estimate_temporal_resolution(f"A_B_{dt_key}"))
+
+    def test_estimate_spatial_resolution_missing_spatial_key(self):
+        self.assertEqual(1, estimate_spatial_resolution(""))
+        self.assertEqual(1, estimate_spatial_resolution("A"))
+        self.assertEqual(1, estimate_spatial_resolution("A_B"))
+        self.assertEqual(1, estimate_spatial_resolution("A_B_C"))
+
+    def test_estimate_spatial_resolution_valid_spatial_key(self):
+        ds = int(random.random() * 10)
+        ds_key = f"DS{ds:03d}"
+
+        self.assertEqual(ds, estimate_spatial_resolution(f"{ds_key}_B_C"))
+        self.assertEqual(ds, estimate_spatial_resolution(f"A_{ds_key}_C"))
+        self.assertEqual(ds, estimate_spatial_resolution(f"A_B_{ds_key}"))
+
+    def test_estimate_spatial_resolution_invalid_spatiall_key(self):
+        ds = int(random.random() * 10)
+        ds_key = f"DS{ds:03d}x"
+
+        self.assertEqual(1, estimate_spatial_resolution(f"{ds_key}_B_C"))
+        self.assertEqual(1, estimate_spatial_resolution(f"A_{ds_key}_C"))
+        self.assertEqual(1, estimate_spatial_resolution(f"A_B_{ds_key}"))