diff --git a/CHANGELOG.md b/CHANGELOG.md index b215a03..c9bb2cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,49 @@ All notable changes to spatialtissuepy are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +Substrate / microenvironment access for PhysiCell simulations — the main +functional gap called out in the v0.2.0 defect report. Additive; nothing +existing changes behavior. + +### Added + +- **Environmental substrate sampling on `PhysiCellTimeStep`.** The + microenvironment `.mat` (oxygen and other diffusible fields) is now reachable + through the high-level API, having previously been parseable only by calling + `parse_microenvironment_mat` directly: + - `substrates` → `{name: (n_voxels,) array}` of concentration fields + - `voxel_positions` → `(n_voxels, 3)` voxel-center coordinates + - `substrate_at(name, x, y, z=None)` → concentration in the voxel nearest each + query point, via a cached `scipy.spatial.cKDTree`. Verified to match a + brute-force `argmin` nearest-voxel lookup exactly on every cell of the + bundled example frame. + - `substrate_names` on both `PhysiCellTimeStep` and `PhysiCellSimulation`. +- **Per-cell internalized substrates.** `PhysiCellTimeStep.internalized_substrates()` + returns a DataFrame of PhysiCell's `internalized_total_substrates` field, one + column per substrate (named), one row per cell, in `positions` order. This is + the amount each cell has actually taken up — accounting for PhysiCell's uptake + dynamics — as opposed to the environmental concentration at the cell's + location that `substrate_at` gives. (Populated only when the simulation + enabled `track_internalized_substrates_in_each_agent`; otherwise present but + zero, which the method documents.) +- `parse_physicell_xml` now records the microenvironment field-data filename at + `metadata.extra['microenvironment_file']`, read from the XML's + `microenvironment/domain/data/filename` element so relocated or renamed + outputs still resolve. The reader falls back to the default + `output{index:08d}_microenvironment0.mat` naming. +- Regression tests (`tests/test_physicell_microenvironment.py`, 22 tests) + covering nearest-voxel equivalence, tie handling, 2-D/3-D and scalar/array + queries, the internalized-vs-raw-matrix mapping, empty-microenvironment + frames, file resolution, and `discover_physicell_timesteps` arity. + +### Notes + +- `discover_physicell_timesteps` keeps its 3-tuple return type unchanged; the + microenvironment file is resolved lazily inside `PhysiCellTimeStep`, so no + caller that unpacks the tuple is affected. + ## [0.3.0] - 2026-07-23 PhysiCell reader correctness fixes, reported against v0.2.0 by the diff --git a/spatialtissuepy/synthetic/physicell/parser.py b/spatialtissuepy/synthetic/physicell/parser.py index 59fc12f..58e9e85 100644 --- a/spatialtissuepy/synthetic/physicell/parser.py +++ b/spatialtissuepy/synthetic/physicell/parser.py @@ -134,6 +134,17 @@ def parse_physicell_xml(xml_path: Path) -> PhysiCellMetadata: name = var.get('name', 'unknown') substrate_names.append(name) + # Extract the microenvironment field-data filename. PhysiCell records it at + # microenvironment/domain/data/filename; this is distinct from the mesh file + # under mesh/voxels/filename. Reading it lets the reader locate the + # substrate .mat even when it does not follow the default naming convention. + microenvironment_file = None + me_elem = root.find('.//microenvironment') + if me_elem is not None: + data_fn = me_elem.find('.//data/filename') + if data_fn is not None and data_fn.text: + microenvironment_file = data_fn.text.strip() + # Extract cell type names and IDs cell_type_names = [] cell_type_ids = [] @@ -170,6 +181,7 @@ def parse_physicell_xml(xml_path: Path) -> PhysiCellMetadata: extra = { 'custom_labels': custom_labels, 'xml_path': str(xml_path), + 'microenvironment_file': microenvironment_file, } return PhysiCellMetadata( diff --git a/spatialtissuepy/synthetic/physicell/reader.py b/spatialtissuepy/synthetic/physicell/reader.py index 6ca8132..839b3b4 100644 --- a/spatialtissuepy/synthetic/physicell/reader.py +++ b/spatialtissuepy/synthetic/physicell/reader.py @@ -11,7 +11,7 @@ import xml.etree.ElementTree as ET from dataclasses import dataclass, field from pathlib import Path -from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union +from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple, Union import numpy as np import pandas as pd @@ -24,9 +24,54 @@ PhysiCellMetadata, get_cell_type_mapping, parse_cells_mat, + parse_microenvironment_mat, parse_physicell_xml, ) + +def _find_microenvironment_mat( + xml_path: Path, + microenvironment_file: Optional[str] = None, +) -> Optional[Path]: + """ + Locate the microenvironment ``.mat`` file for a frame. + + Resolution order: + + 1. The filename declared in the frame's XML (``microenvironment/domain/ + data/filename``), if given -- this survives renamed or relocated output. + 2. The default naming convention ``output{index:08d}_microenvironment0.mat`` + next to the XML. + + Parameters + ---------- + xml_path : Path + Path to the frame's ``output*.xml``. + microenvironment_file : str, optional + Filename from the XML, as parsed into + ``metadata.extra['microenvironment_file']``. + + Returns + ------- + Path or None + Path to the microenvironment ``.mat`` if found, else ``None``. + """ + xml_path = Path(xml_path) + + if microenvironment_file: + candidate = xml_path.parent / microenvironment_file + if candidate.exists(): + return candidate + + match = re.search(r'output(\d+)', xml_path.stem) + if match: + index = int(match.group(1)) + candidate = xml_path.parent / f'output{index:08d}_microenvironment0.mat' + if candidate.exists(): + return candidate + + return None + # ----------------------------------------------------------------------------- # PhysiCell TimeStep # ----------------------------------------------------------------------------- @@ -69,8 +114,13 @@ class PhysiCellTimeStep(ABMTimeStep): cells_mat_path: Path = None cell_type_mapping: Dict[int, str] = field(default_factory=dict) include_dead_cells: bool = False + microenvironment_mat_path: Optional[Path] = None _cell_data: Optional[Dict[str, np.ndarray]] = field(default=None, repr=False) _physicell_metadata: Optional[PhysiCellMetadata] = field(default=None, repr=False) + _microenvironment: Optional[Dict[str, Any]] = field(default=None, repr=False) + _me_loaded: bool = field(default=False, repr=False) + _voxel_tree: Any = field(default=None, repr=False) + _voxel_tree_ndim: Optional[int] = field(default=None, repr=False) def _load_cell_data(self) -> Dict[str, np.ndarray]: """Load cell data from MAT file (cached).""" @@ -99,6 +149,200 @@ def _load_metadata(self) -> PhysiCellMetadata: self._physicell_metadata = parse_physicell_xml(self.source_path) return self._physicell_metadata + def _load_microenvironment(self) -> Dict[str, Any]: + """ + Load the substrate field data for this frame (cached). + + Resolves the microenvironment ``.mat`` from the XML (or the default + naming convention) if a path was not supplied. A frame with no + microenvironment file yields empty substrates rather than raising, so + callers can probe ``substrates`` unconditionally. + """ + if self._me_loaded: + return self._microenvironment + + empty = { + 'voxel_positions': np.empty((0, 3)), + 'concentrations': {}, + 'raw_data': None, + } + + path = self.microenvironment_mat_path + if path is None: + try: + me_file = self._load_metadata().extra.get('microenvironment_file') + except (OSError, ET.ParseError): + me_file = None + path = _find_microenvironment_mat(self.source_path, me_file) + self.microenvironment_mat_path = path + + if path is None: + self._microenvironment = empty + else: + try: + substrate_names = self._load_metadata().substrate_names + except (OSError, ET.ParseError): + substrate_names = None + self._microenvironment = parse_microenvironment_mat( + path, substrate_names + ) + + self._me_loaded = True + return self._microenvironment + + @property + def substrate_names(self) -> List[str]: + """Names of the diffusible substrates, in matrix order.""" + return list(self._load_metadata().substrate_names) + + @property + def substrates(self) -> Dict[str, np.ndarray]: + """ + Substrate concentration fields, keyed by name. + + Each value is a ``(n_voxels,)`` array aligned with + :attr:`voxel_positions`. Empty if the frame has no microenvironment + file. + """ + return self._load_microenvironment()['concentrations'] + + @property + def voxel_positions(self) -> np.ndarray: + """``(n_voxels, 3)`` array of voxel-center coordinates.""" + return self._load_microenvironment()['voxel_positions'] + + def substrate_at( + self, + name: str, + x: Union[float, np.ndarray], + y: Union[float, np.ndarray], + z: Optional[Union[float, np.ndarray]] = None, + ) -> np.ndarray: + """ + Sample a substrate's environmental concentration at spatial points. + + Returns the concentration in the voxel nearest each query point, found + with a KD-tree over the voxel centers. This is the concentration of the + substrate *in the space around* a location; for the amount a cell has + actually taken up, use :meth:`internalized_substrates`. + + Parameters + ---------- + name : str + Substrate name; must be one of :attr:`substrate_names`. + x, y : float or array-like + Query coordinates. Scalars or equal-length arrays. + z : float or array-like, optional + Query z. If given, nearest-voxel search is 3-D; otherwise it uses + the x-y plane (appropriate for 2-D simulations). + + Returns + ------- + np.ndarray + Concentration at each query point, matching the input shape. + """ + me = self._load_microenvironment() + concentrations = me['concentrations'] + if name not in concentrations: + available = ', '.join(concentrations) or '(none)' + raise ValueError( + f"Unknown substrate {name!r}. Available: {available}" + ) + + voxels = me['voxel_positions'] + if voxels.shape[0] == 0: + raise ValueError( + f"No microenvironment voxels available for {self.source_path}" + ) + + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + if z is None: + query = np.column_stack([np.ravel(x), np.ravel(y)]) + ndim = 2 + else: + z = np.asarray(z, dtype=float) + query = np.column_stack([np.ravel(x), np.ravel(y), np.ravel(z)]) + ndim = 3 + + tree = self._voxel_kdtree(ndim) + _, idx = tree.query(query) + result = concentrations[name][idx] + + return result if x.ndim else result[0] + + def _voxel_kdtree(self, ndim: int): + """Build (and cache) a KD-tree over voxel centers for `ndim` axes.""" + if self._voxel_tree is None or self._voxel_tree_ndim != ndim: + from scipy.spatial import cKDTree + + voxels = self._load_microenvironment()['voxel_positions'] + self._voxel_tree = cKDTree(voxels[:, :ndim]) + self._voxel_tree_ndim = ndim + return self._voxel_tree + + def internalized_substrates( + self, + include_dead_cells: Optional[bool] = None, + ) -> pd.DataFrame: + """ + Per-cell internalized substrate amounts, one column per substrate. + + These come from PhysiCell's ``internalized_total_substrates`` field on + each cell, so they reflect the framework's actual uptake dynamics rather + than the environmental concentration at the cell's location. Use this + instead of :meth:`substrate_at` when you need what a cell has taken up, + not what surrounds it. + + Rows are in the same order as :attr:`positions`. + + .. note:: + PhysiCell only accumulates this field when the simulation enables + ``track_internalized_substrates_in_each_agent``. If that option was + off, the column is present but uniformly zero -- that is the model's + recorded value, not a parsing artifact. For the environmental + concentration, which is always available, use :meth:`substrate_at`. + + Parameters + ---------- + include_dead_cells : bool, optional + Whether to include dead cells. Defaults to the instance attribute. + + Returns + ------- + pd.DataFrame + Columns named by substrate; one row per cell. + + Raises + ------ + ValueError + If the frame's cell data does not record internalized substrates + (the model did not write the ``internalized_total_substrates`` + field, or the XML labels were unavailable). + """ + data = self._load_cell_data() + columns = data.get('columns') or {} + names = self.substrate_names + + field_keys = [f'internalized_total_substrates_{i}' for i in range(len(names))] + if not names or not all(k in columns for k in field_keys): + raise ValueError( + "This frame does not record internalized substrates. It " + "requires the PhysiCell 'internalized_total_substrates' field " + "and readable XML ; neither is optional here." + ) + + if include_dead_cells is None: + include_dead_cells = self.include_dead_cells + if include_dead_cells: + mask = np.ones(len(data['cell_types']), dtype=bool) + else: + mask = data['dead_flags'] == 0 + + return pd.DataFrame( + {name: columns[key][mask] for name, key in zip(names, field_keys)} + ) + @property def n_cells(self) -> int: """Number of cells at this time step.""" @@ -390,6 +634,18 @@ def n_timesteps(self) -> int: """Number of time steps.""" return len(self._timestep_files) + @property + def substrate_names(self) -> List[str]: + """ + Names of the diffusible substrates in this simulation. + + Read from the first frame's XML; empty if there are no time steps. + """ + if not self._timestep_files: + return [] + first_xml = self._timestep_files[0][1] + return list(parse_physicell_xml(first_xml).substrate_names) + @property def times(self) -> np.ndarray: """Array of simulation times.""" diff --git a/tests/test_physicell_microenvironment.py b/tests/test_physicell_microenvironment.py new file mode 100644 index 0000000..a60abc1 --- /dev/null +++ b/tests/test_physicell_microenvironment.py @@ -0,0 +1,237 @@ +""" +Tests for PhysiCell microenvironment / substrate access (v0.4.0). + +Covers two distinct notions of "how much substrate": + +- environmental concentration at a location, sampled from the voxel field via + ``PhysiCellTimeStep.substrate_at`` (nearest-voxel, KD-tree); +- the amount a cell has actually internalized, read per-cell from PhysiCell's + ``internalized_total_substrates`` field via + ``PhysiCellTimeStep.internalized_substrates``. + +The bundled example simulation is used for the equivalence checks; small +synthetic inputs cover ties, missing files, and error paths. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.io import loadmat, savemat + +from spatialtissuepy.synthetic.physicell import ( + PhysiCellTimeStep, + discover_physicell_timesteps, + read_physicell_simulation, + read_physicell_timestep, +) +from spatialtissuepy.synthetic.physicell.reader import _find_microenvironment_mat + + +@pytest.fixture +def example_physicell_dir(): + path = ( + Path(__file__).parent.parent + / 'examples' / 'sample_data' / 'example_physicell_sim' + ) + if not path.exists(): + pytest.skip(f"Example PhysiCell data not found at {path}") + return path + + +@pytest.fixture +def frame_xml(example_physicell_dir): + xmls = sorted(example_physicell_dir.glob('output*.xml')) + if not xmls: + pytest.skip("No frames in example simulation") + return xmls[len(xmls) // 2] + + +def _injected_timestep(voxels, concentrations): + """A PhysiCellTimeStep with a microenvironment injected, no files read.""" + ts = PhysiCellTimeStep(time=0.0, time_index=0, source_path=Path('none.xml')) + ts._microenvironment = { + 'voxel_positions': np.asarray(voxels, dtype=float), + 'concentrations': {k: np.asarray(v, dtype=float) + for k, v in concentrations.items()}, + 'raw_data': None, + } + ts._me_loaded = True + return ts + + +class TestSubstrateFields: + def test_substrate_names(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert 'oxygen' in ts.substrate_names + + def test_substrates_shape(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + oxygen = ts.substrates['oxygen'] + assert oxygen.ndim == 1 + assert oxygen.shape[0] == ts.voxel_positions.shape[0] + + def test_voxel_positions_are_3d(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert ts.voxel_positions.shape[1] == 3 + + +class TestSubstrateAt: + def test_matches_bruteforce_argmin(self, frame_xml): + """The report's acceptance criterion: exact match to nearest-voxel.""" + ts = read_physicell_timestep(frame_xml) + coords = ts.to_spatial_data().coordinates + xs, ys = coords[:, 0], coords[:, 1] + + got = ts.substrate_at('oxygen', xs, ys) + + vox = ts.voxel_positions[:, :2] + field = ts.substrates['oxygen'] + brute = np.array([ + field[np.argmin((vox[:, 0] - x) ** 2 + (vox[:, 1] - y) ** 2)] + for x, y in zip(xs, ys) + ]) + np.testing.assert_array_equal(got, brute) + + def test_scalar_query_returns_scalar(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + val = ts.substrate_at('oxygen', 0.0, 0.0) + assert np.isscalar(val) or val.ndim == 0 + + def test_array_query_preserves_length(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + out = ts.substrate_at('oxygen', np.array([0.0, 10.0, 20.0]), + np.array([0.0, 10.0, 20.0])) + assert out.shape == (3,) + + def test_unknown_substrate_raises(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + with pytest.raises(ValueError, match="Unknown substrate"): + ts.substrate_at('unobtanium', 0.0, 0.0) + + def test_tie_is_deterministic_and_valid(self): + """A point equidistant between two voxels returns one of their values. + + KD-tree tie-breaking need not match np.argmin, so this asserts the + result is one of the tied voxels and is stable, rather than pinning a + specific choice. + """ + voxels = [[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]] + ts = _injected_timestep(voxels, {'oxygen': [1.0, 2.0]}) + + first = ts.substrate_at('oxygen', 5.0, 0.0) + second = ts.substrate_at('oxygen', 5.0, 0.0) + + assert first in (1.0, 2.0) + assert first == second # deterministic + + def test_3d_query(self): + voxels = [[0, 0, 0], [0, 0, 10]] + ts = _injected_timestep(voxels, {'o2': [5.0, 9.0]}) + assert ts.substrate_at('o2', 0.0, 0.0, 1.0) == 5.0 + assert ts.substrate_at('o2', 0.0, 0.0, 9.0) == 9.0 + + +class TestInternalizedSubstrates: + def test_columns_match_raw_matrix(self, frame_xml, example_physicell_dir): + """internalized_total_substrates maps to rows 92+ in matrix order.""" + ts = read_physicell_timestep(frame_xml, include_dead_cells=True) + matrix = loadmat(str(ts.cells_mat_path))['cells'] + if matrix.shape[0] < 97: + pytest.skip("Frame lacks internalized_total_substrates rows") + + inte = ts.internalized_substrates(include_dead_cells=True) + for i, name in enumerate(ts.substrate_names): + np.testing.assert_array_equal(inte[name].to_numpy(), matrix[92 + i, :]) + + def test_columns_named_by_substrate(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + inte = ts.internalized_substrates() + assert list(inte.columns) == ts.substrate_names + + def test_row_count_follows_dead_filter(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert len(ts.internalized_substrates()) == ts.n_cells + assert len(ts.internalized_substrates(include_dead_cells=True)) == \ + ts.n_cells_total + + def test_row_order_matches_positions(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + inte = ts.internalized_substrates() + assert len(inte) == ts.positions.shape[0] + + def test_raises_without_field(self, tmp_path): + """A frame whose cell data lacks the field raises informatively.""" + # Wider than tall so orientation is unambiguous; no , so the + # parser exposes no 'columns' and no substrate names are declared. + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((5, 20))}) + (tmp_path / 'output00000000.xml').write_text( + "" + ) + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + with pytest.raises(ValueError, match="internalized"): + ts.internalized_substrates() + + +class TestFileResolution: + def test_resolves_from_xml_filename(self, frame_xml, example_physicell_dir): + expected = example_physicell_dir / f'{frame_xml.stem}_microenvironment0.mat' + got = _find_microenvironment_mat( + frame_xml, f'{frame_xml.stem}_microenvironment0.mat' + ) + assert got == expected + + def test_falls_back_to_naming_convention(self, frame_xml): + # No filename hint -> should still find it by convention. + got = _find_microenvironment_mat(frame_xml, None) + assert got is not None and got.exists() + + def test_returns_none_when_absent(self, tmp_path): + (tmp_path / 'output00000000.xml').write_text("") + assert _find_microenvironment_mat( + tmp_path / 'output00000000.xml', None + ) is None + + +class TestEmptyMicroenvironment: + def test_missing_file_yields_empty(self, tmp_path): + """No microenvironment file -> empty substrates, no exception.""" + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((30, 4))}) + (tmp_path / 'output00000000.xml').write_text( + "" + ) + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + + assert ts.substrates == {} + assert ts.voxel_positions.shape == (0, 3) + + def test_substrate_at_on_empty_raises(self, tmp_path): + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((30, 4))}) + (tmp_path / 'output00000000.xml').write_text("") + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + with pytest.raises(ValueError): + ts.substrate_at('oxygen', 0.0, 0.0) + + +class TestSimulationIntegration: + def test_substrate_names(self, example_physicell_dir): + sim = read_physicell_simulation(example_physicell_dir) + assert 'oxygen' in sim.substrate_names + + def test_timesteps_have_substrate_access(self, example_physicell_dir): + sim = read_physicell_simulation(example_physicell_dir) + ts = sim.get_timestep(0) + assert ts.voxel_positions.shape[0] > 0 + + +class TestBackwardCompatibility: + def test_discover_still_returns_3_tuples(self, example_physicell_dir): + """Adding microenvironment support must not change discover's shape.""" + steps = discover_physicell_timesteps(example_physicell_dir) + assert steps and all(len(s) == 3 for s in steps) + index, xml_path, mat_path = steps[0] # unpacks as before + assert xml_path.suffix == '.xml' + assert mat_path.suffix == '.mat'