diff --git a/.gitignore b/.gitignore index d157ddaae..f44cfbecf 100644 --- a/.gitignore +++ b/.gitignore @@ -54,3 +54,29 @@ _version.py # benchmarking and profiling .asv/ profile.speedscope.json + +# other +node_modules/ + +.mypy_cache +.ruff_cache + +# Kilo and plans +.kilo/ +plans/ + +# test coverage +.coverage +htmlcov/ + +# local hatch config +hatch.toml + +# Kilo and plans +.kilo/ +plans/ + +# test coverage +.coverage +htmlcov/ + diff --git a/docs/api/models_utils.md b/docs/api/models_utils.md index c2a74803b..071cd8ba7 100644 --- a/docs/api/models_utils.md +++ b/docs/api/models_utils.md @@ -7,7 +7,6 @@ .. autofunction:: validate_element .. autodata:: SpatialElement .. autofunction:: get_axes_names -.. autofunction:: get_spatial_axes .. autofunction:: points_geopandas_to_dask_dataframe .. autofunction:: points_dask_dataframe_to_geopandas .. autofunction:: get_channel_names diff --git a/docs/api/transformations_utils.md b/docs/api/transformations_utils.md index dafa228b6..307b4aee4 100644 --- a/docs/api/transformations_utils.md +++ b/docs/api/transformations_utils.md @@ -10,4 +10,5 @@ .. autofunction:: get_transformation_between_landmarks .. autofunction:: align_elements_using_landmarks .. autofunction:: remove_transformations_to_coordinate_system +.. autofunction:: get_spatial_axes ``` diff --git a/pyproject.toml b/pyproject.toml index 4cc02bcb4..13ab89cc6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,8 @@ dependencies = [ "networkx", "numba>=0.55", "numpy", - "ome-zarr>=0.16", + "ome-zarr>=0.19.2", + "ome-zarr-models>=1.8", "pandas", "pooch", "pyarrow", @@ -51,7 +52,7 @@ dependencies = [ "spatial-image>=1.2.3", "typing-extensions>=4.8", "universal-pathlib>=0.2.6", - "xarray>=2024.10", + "xarray>=2025.4", "xarray-spatial>=0.3.5", "zarr>=3", ] diff --git a/src/spatialdata/__init__.py b/src/spatialdata/__init__.py index 96e7032dd..82292e65f 100644 --- a/src/spatialdata/__init__.py +++ b/src/spatialdata/__init__.py @@ -4,7 +4,10 @@ from importlib.metadata import version from typing import TYPE_CHECKING, Any -import spatialdata.models._accessor +from spatialdata.models._accessor import register_attrs_accessor_dask_df_series + +# Import registers Dask .attrs accessor, enabling metadata handling on Dask collections +register_attrs_accessor_dask_df_series() __version__ = version("spatialdata") @@ -59,7 +62,6 @@ "polygon_query": "spatialdata._core.query.spatial_query", # _core.spatialdata "SpatialData": "spatialdata._core.spatialdata", - # _io._utils "get_dask_backing_files": "spatialdata._io._utils", # _io.format "SpatialDataFormatType": "spatialdata._io.format", diff --git a/src/spatialdata/_core/_deepcopy.py b/src/spatialdata/_core/_deepcopy.py index 2a68ab7a2..ec3f3e7ae 100644 --- a/src/spatialdata/_core/_deepcopy.py +++ b/src/spatialdata/_core/_deepcopy.py @@ -11,14 +11,14 @@ from xarray import DataArray, DataTree from spatialdata._core.spatialdata import SpatialData -from spatialdata.models._utils import SpatialElement -from spatialdata.models.models import ( +from spatialdata.models import ( Image2DModel, Image3DModel, Labels2DModel, Labels3DModel, PointsModel, RasterSchema, + SpatialElement, get_model, ) diff --git a/src/spatialdata/_core/centroids.py b/src/spatialdata/_core/centroids.py index 9979b9f20..65c361b3d 100644 --- a/src/spatialdata/_core/centroids.py +++ b/src/spatialdata/_core/centroids.py @@ -11,9 +11,7 @@ from xarray import DataArray, DataTree from spatialdata._core.operations.transform import transform -from spatialdata.models import get_axes_names -from spatialdata.models._utils import SpatialElement -from spatialdata.models.models import Labels2DModel, Labels3DModel, PointsModel, get_model +from spatialdata.models import Labels2DModel, Labels3DModel, PointsModel, SpatialElement, get_axes_names, get_model from spatialdata.transformations.operations import get_transformation from spatialdata.transformations.transformations import BaseTransformation diff --git a/src/spatialdata/_core/data_extent.py b/src/spatialdata/_core/data_extent.py index 32f05834b..5efdcfdb8 100644 --- a/src/spatialdata/_core/data_extent.py +++ b/src/spatialdata/_core/data_extent.py @@ -14,9 +14,7 @@ from spatialdata._core.operations.transform import transform from spatialdata._core.spatialdata import SpatialData -from spatialdata.models import get_axes_names -from spatialdata.models._utils import SpatialElement -from spatialdata.models.models import PointsModel +from spatialdata.models import PointsModel, SpatialElement, get_axes_names from spatialdata.transformations.operations import get_transformation BoundingBoxDescription = dict[str, tuple[float, float]] diff --git a/src/spatialdata/_core/operations/_utils.py b/src/spatialdata/_core/operations/_utils.py index 600dbed75..cd097e5d8 100644 --- a/src/spatialdata/_core/operations/_utils.py +++ b/src/spatialdata/_core/operations/_utils.py @@ -5,7 +5,8 @@ from anndata import AnnData from xarray import DataArray, DataTree -from spatialdata.models import SpatialElement, get_axes_names, get_spatial_axes +from spatialdata.models import SpatialElement, get_axes_names +from spatialdata.transformations import get_spatial_axes if TYPE_CHECKING: from spatialdata._core.spatialdata import SpatialData diff --git a/src/spatialdata/_core/operations/map.py b/src/spatialdata/_core/operations/map.py index b0738dbff..4dc860133 100644 --- a/src/spatialdata/_core/operations/map.py +++ b/src/spatialdata/_core/operations/map.py @@ -13,7 +13,7 @@ from xarray import DataArray, DataTree from spatialdata._types import IntArrayLike -from spatialdata.models._utils import get_axes_names, get_channel_names, get_raster_model_from_data_dims +from spatialdata.models import get_axes_names, get_channel_names, get_raster_model_from_data_dims from spatialdata.transformations import get_transformation __all__ = ["map_raster", "relabel_sequential"] diff --git a/src/spatialdata/_core/operations/rasterize.py b/src/spatialdata/_core/operations/rasterize.py index d47f4bfca..eac8fef10 100644 --- a/src/spatialdata/_core/operations/rasterize.py +++ b/src/spatialdata/_core/operations/rasterize.py @@ -31,7 +31,7 @@ get_axes_names, get_model, ) -from spatialdata.models._utils import get_spatial_axes +from spatialdata.transformations import get_spatial_axes from spatialdata.transformations._utils import _get_scale, compute_coordinates from spatialdata.transformations.operations import get_transformation, remove_transformation, set_transformation from spatialdata.transformations.transformations import ( diff --git a/src/spatialdata/_core/operations/transform.py b/src/spatialdata/_core/operations/transform.py index 410e92fc3..421715bef 100644 --- a/src/spatialdata/_core/operations/transform.py +++ b/src/spatialdata/_core/operations/transform.py @@ -20,9 +20,14 @@ from spatialdata._core.spatialdata import SpatialData from spatialdata._types import ArrayLike from spatialdata._utils import disable_dask_tune_optimization -from spatialdata.models import SpatialElement, get_axes_names, get_model -from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM, get_channel_names -from spatialdata.models.models import RasterSchema +from spatialdata.literals import DEFAULT_COORDINATE_SYSTEM, TRANSFORM_KEY +from spatialdata.models import ( + RasterSchema, + SpatialElement, + get_axes_names, + get_channel_names, + get_model, +) from spatialdata.transformations._utils import _get_scale, compute_coordinates, scale_radii if TYPE_CHECKING: @@ -356,7 +361,6 @@ def _( schema = get_model(data) assert issubclass(schema, RasterSchema) from spatialdata.models import Image2DModel, Image3DModel, Labels2DModel, Labels3DModel - from spatialdata.models._utils import TRANSFORM_KEY from spatialdata.transformations import get_transformation, set_transformation from spatialdata.transformations.transformations import Identity, Sequence @@ -441,7 +445,6 @@ def _( to_coordinate_system: str | None = None, ) -> DaskDataFrame: from spatialdata.models import PointsModel - from spatialdata.models._utils import TRANSFORM_KEY from spatialdata.transformations import Identity, get_transformation transformation, to_coordinate_system = _validate_target_coordinate_systems( @@ -519,7 +522,6 @@ def _( to_coordinate_system: str | None = None, ) -> GeoDataFrame: from spatialdata.models import ShapesModel - from spatialdata.models._utils import TRANSFORM_KEY from spatialdata.transformations import Identity, get_transformation axes = get_axes_names(data) diff --git a/src/spatialdata/_core/operations/vectorize.py b/src/spatialdata/_core/operations/vectorize.py index 431bdbd17..14a16a072 100644 --- a/src/spatialdata/_core/operations/vectorize.py +++ b/src/spatialdata/_core/operations/vectorize.py @@ -30,8 +30,8 @@ SpatialElement, get_axes_names, get_model, + points_dask_dataframe_to_geopandas, ) -from spatialdata.models._utils import points_dask_dataframe_to_geopandas from spatialdata.transformations.operations import get_transformation from spatialdata.transformations.transformations import Identity diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index f0360c310..d31b977fb 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -21,6 +21,7 @@ from spatialdata._types import ArrayLike, ListOrNDArrayFloating, Raster_T from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( + ATTRS_KEY, PointsModel, ShapesModel, SpatialElement, @@ -28,8 +29,8 @@ points_dask_dataframe_to_geopandas, points_geopandas_to_dask_dataframe, ) -from spatialdata.models._utils import ValidAxis_t, get_spatial_axes -from spatialdata.models.models import ATTRS_KEY +from spatialdata.transformations import get_spatial_axes +from spatialdata.transformations._graph.types import ValidAxis_t from spatialdata.transformations.operations import set_transformation from spatialdata.transformations.transformations import ( Affine, diff --git a/src/spatialdata/_core/spatialdata.py b/src/spatialdata/_core/spatialdata.py index 9f5c91c7a..3d162b77b 100644 --- a/src/spatialdata/_core/spatialdata.py +++ b/src/spatialdata/_core/spatialdata.py @@ -4,7 +4,7 @@ import json import os import warnings -from collections.abc import Generator, Mapping +from collections.abc import Generator, Iterable, Mapping from itertools import chain from pathlib import Path from typing import TYPE_CHECKING, Any, Literal, cast @@ -22,6 +22,7 @@ from zarr.errors import GroupNotFoundError from spatialdata._core._elements import Images, Labels, Points, Shapes, Tables +from spatialdata._core.transformation_manager import TransformationsFacade from spatialdata._core.validation import ( check_all_keys_case_insensitively_unique, check_target_region_column_symmetry, @@ -39,16 +40,15 @@ Labels3DModel, PointsModel, ShapesModel, - TableModel, - get_model, - get_table_keys, -) -from spatialdata.models._utils import ( SpatialElement, + TableModel, convert_region_column_to_categorical, get_axes_names, + get_model, + get_table_keys, set_channel_names, ) +from spatialdata.transformations import BaseTransformation, get_transformation if TYPE_CHECKING: from pandas._typing import DtypeObj @@ -132,6 +132,7 @@ def __init__( self._shapes: Shapes = Shapes(shared_keys=self._shared_keys) self._tables: Tables = Tables(shared_keys=self._shared_keys) self.attrs = attrs if attrs else {} # type: ignore[assignment] + self.transformations = TransformationsFacade() element_names = list(chain.from_iterable([e.keys() for e in [images, labels, points, shapes] if e is not None])) @@ -172,6 +173,59 @@ def __init__( self.validate_table_in_spatialdata(table) self.tables[k] = table + # to be deleted later, when cleaning up parts that rely on transformation info in attrs + for element_set in (images, labels, shapes, points): + if element_set is not None: + for element_key, element in element_set.items(): + self.transformations._add_new_element_to_graph_creating_default_cs( + new_element=element, new_element_key=element_key + ) + + all_element_dict = (images or {}) | (labels or {}) | (shapes or {}) | (points or {}) + elements_grouped_by_transformations = self.group_elements_by_target_cs(all_element_dict.keys()) + + for target_cs_name, element_name_transformation_config in elements_grouped_by_transformations.items(): + self.transformations.set_transformations_from_elements_to_common_cs( + target_cs_name=target_cs_name, + element_name_transformation_config=element_name_transformation_config, + ) + + def group_elements_by_target_cs( + self, spatial_element_names: Iterable[str] + ) -> dict[str, dict[str, BaseTransformation]]: + """ + Groups SpatialElement objects by target coordinate systems + + Based on the transformation information in their .attrs + + Parameters + ---------- + spatial_element_names + names of SpatialElement objects + + Returns + ------- + dictionary with target coordinate systems as keys and element-wise transformation information as values. + These values are themselves dictionaries with element_names as keys and transformation objects mapping + element to the corresponding target coordinate system as values + """ + spatial_elements_grouped_by_transformations: dict[str, dict[str, BaseTransformation]] = {} + for element_name in spatial_element_names: + element = self.get(element_name) + if (element is not None) and isinstance(element, SpatialElement): + transformations_config = get_transformation(element, get_all=True) + for target_cs_name, transformation_object in transformations_config.items(): + if target_cs_name in spatial_elements_grouped_by_transformations: + spatial_elements_grouped_by_transformations[target_cs_name][element_name] = ( + transformation_object + ) + else: + spatial_elements_grouped_by_transformations[target_cs_name] = { + element_name: transformation_object + } + + return spatial_elements_grouped_by_transformations + def validate_table_in_spatialdata(self, table: AnnData) -> None: """ Validate the presence of the annotation target of a SpatialData table in the SpatialData object. @@ -2368,6 +2422,12 @@ def _find_element(self, element_name: str) -> tuple[str, str, SpatialElement | A return found[0] + def assert_element_does_not_exist(self, element_name: str) -> None: + """Raises a key error if an element with name `element_name` already exists.""" + existing_element_names = [x[1] for x in self.gen_elements()] + if element_name in existing_element_names: + raise KeyError(f"Element with name '{element_name}' already exists") + @classmethod def init_from_elements( cls, @@ -2524,6 +2584,8 @@ def __setitem__(self, key: str, value: SpatialElement | AnnData) -> None: """ Add the element to the SpatialData object. + Also create a "default" coordinate system for the element and link the element to it. + Parameters ---------- key @@ -2535,19 +2597,49 @@ def __setitem__(self, key: str, value: SpatialElement | AnnData) -> None: schema = get_model(value) if schema in (Image2DModel, Image3DModel): assert isinstance(value, DataArray | DataTree) + self.images[key] = value + # this will throw KeyError if we try to assign an element of different type, but with the same name as an + # existing spatial element; will overwrite if it is the same type + self.transformations._remove_element_from_graph_if_exists(key) + self.transformations._add_new_element_to_graph_creating_default_cs(new_element=value, new_element_key=key) + # TODO add transformations from element attrs to graph + elif schema in (Labels2DModel, Labels3DModel): assert isinstance(value, DataArray | DataTree) + self.labels[key] = value + # this will throw KeyError if we try to assign an element of different type, but with the same name as an + # existing spatial element; will overwrite if it is the same type + self.transformations._remove_element_from_graph_if_exists(key) + self.transformations._add_new_element_to_graph_creating_default_cs(new_element=value, new_element_key=key) + # TODO add transformations from element attrs to graph + elif schema == PointsModel: assert isinstance(value, DaskDataFrame) + self.points[key] = value + # this will throw KeyError if we try to assign an element of different type, but with the same name as an + # existing spatial element; will overwrite if it is the same type + self.transformations._remove_element_from_graph_if_exists(key) + self.transformations._add_new_element_to_graph_creating_default_cs(new_element=value, new_element_key=key) + # TODO add transformations from element attrs to graph + elif schema == ShapesModel: assert isinstance(value, GeoDataFrame) + self.shapes[key] = value + # this will throw KeyError if we try to assign an element of different type, but with the same name as an + # existing spatial element; will overwrite if it is the same type + self.transformations._remove_element_from_graph_if_exists(key) + self.transformations._add_new_element_to_graph_creating_default_cs(new_element=value, new_element_key=key) + # TODO add transformations from element attrs to graph + elif schema == TableModel: assert isinstance(value, AnnData) + self.validate_table_in_spatialdata(value) self.tables[key] = value + # transformation data is not stored in AnnData.attrs else: raise TypeError(f"Unknown element type with schema: {schema!r}.") @@ -2562,6 +2654,7 @@ def __delitem__(self, key: str) -> None: """ element_type, _, _ = self._find_element(key) getattr(self, element_type).__delitem__(key) + self.transformations._remove_element_from_graph_if_exists(key) @property def attrs(self) -> dict[str, JSONValue]: diff --git a/src/spatialdata/_core/transformation_manager/__init__.py b/src/spatialdata/_core/transformation_manager/__init__.py new file mode 100644 index 000000000..e4af58490 --- /dev/null +++ b/src/spatialdata/_core/transformation_manager/__init__.py @@ -0,0 +1,5 @@ +from __future__ import annotations + +from spatialdata._core.transformation_manager.transformations_facade import TransformationsFacade + +__all__ = ["TransformationsFacade"] diff --git a/src/spatialdata/_core/transformation_manager/transformations_facade.py b/src/spatialdata/_core/transformation_manager/transformations_facade.py new file mode 100644 index 000000000..4ca9fa5a2 --- /dev/null +++ b/src/spatialdata/_core/transformation_manager/transformations_facade.py @@ -0,0 +1,243 @@ +from __future__ import annotations + +from collections.abc import Sequence +from itertools import chain + +from spatialdata.models import SpatialElement +from spatialdata.transformations import BaseTransformation +from spatialdata.transformations._graph import ( + BaseTransformationEdge, + CoordinateSystemNotFoundError, + CoordSystem, + ElementNotRegisteredToAnyCoordinateSystemError, + ProjectAxisEdge, + TransformationGraph, + get_instrinsic_coordinate_system, +) +from spatialdata.transformations._graph.utils import get_intrinsic_coordinate_system_name, get_union_of_axes + + +class TransformationsFacade: + def __init__(self) -> None: + """Initialize a TransformationsFacade with an empty transformation graph.""" + self._graph = TransformationGraph() + + def _get_element_transformation_edges( + self, element_name: str, target_coordinate_system_names: Sequence[str] | None = None + ) -> dict[str, BaseTransformationEdge]: + """ + Get the NGFF transformation edges associated with an element. + + Parameters + ---------- + element_name + The name of the element. + target_coordinate_system_names + Optional sequence of target coordinate system names to filter by. + + Returns + ------- + dict, with edge keys as keys and transformations from the element's coordinate system as values. + """ + return self._graph.get_outgoing_edges(element_name, target_coordinate_system_names) + + def _add_new_element_to_graph_creating_default_cs(self, new_element: SpatialElement, new_element_key: str) -> None: + """ + Add a new element to the transformation graph with a default coordinate system. + + Parameters + ---------- + new_element + The spatial element to add to the graph. + new_element_key + The key/name to use for the new element in the graph. + """ + element_cs = get_instrinsic_coordinate_system( + element_to_be_added_to_graph=new_element, + coordinate_system_name=get_intrinsic_coordinate_system_name(new_element_key), + ) + + self._graph.add_coordinate_system(element_cs) + self._graph.add_element(element_name=new_element_key, coordinate_system=element_cs) + + def _remove_element_from_graph_if_exists(self, element_name: str) -> None: + """ + Remove element from the graph if it exists, do nothing otherwise + + when removing also removes all incoming and outgoing edges from the coordinate system to which the + element belongs, then the coordinate system itself as well as the element-coordinate system mapping + + Parameters + ---------- + element_name + The name of the element to remove. + """ + try: + element_coordinate_system = self._graph.get_element_coordinate_system(element_name) + except ElementNotRegisteredToAnyCoordinateSystemError as _err: + return # element does not exist, do nothing + + number_of_elements_in_element_cs = len( + self._graph.get_elements_in_coordinate_system(element_coordinate_system.name) + ) + + if number_of_elements_in_element_cs == 1: + for incoming_edge_key, incoming_edge in self._graph.get_incoming_edges(element_name).items(): + self._graph.remove_specific_transformation(incoming_edge_key, incoming_edge) + + for outgoing_edge_key, outgoing_edge in self._graph.get_outgoing_edges(element_name).items(): + self._graph.remove_specific_transformation(outgoing_edge_key, outgoing_edge) + + self._graph._unset_element(element_name) + self._graph.remove_coordinate_system(element_coordinate_system) + else: + self._graph._unset_element(element_name) + + def set_transformations_from_elements_to_common_cs( + self, + target_cs_name: str, + element_name_transformation_config: dict[str, BaseTransformation], + make_target_cs_virtual_when_creating: bool = False, + ) -> None: + """ + Infer and add the transformation edges needed to align a group of spatial elements to a common coordinate system + + In addition to the transformations specified, adds project axis edges as needed. + + Parameters + ---------- + target_cs_name + Name of the target coordinate system. Will be created if it doesn't exist. + element_name_transformation_config + Dictionary mapping names of spatial elements to transformations that map the corresponding element to + the target coordinate system. + make_target_cs_virtual_when_creating + whether to make target coordinate system virtual, if and when it gets created. + """ + intermediate_edge_list: list[BaseTransformationEdge] = [] + for element_name, transformation in element_name_transformation_config.items(): + element_cs = self._graph.get_element_coordinate_system(element_name) + intermediate_edge = transformation._to_ngff_transformation_edge_creating_output_cs( + input_coordinate_system=element_cs, + output_coordinate_system_name=f"{target_cs_name}_{element_name}", + output_coordinate_system_virtual=True, + ) + # if a project axis edge is needed for this element, this edge can be used below + # hence it has a virtual output cs with a specially encoded name. + intermediate_edge_list.append(intermediate_edge) + + expected_output_axes_flat = tuple(chain.from_iterable(x.output.axes for x in intermediate_edge_list)) + # basically flattens a tuple[tuple[T, ...], ...] into tuple[T] + expected_output_axes_union = get_union_of_axes(expected_output_axes_flat) + + target_cs = self._graph.get_coordinate_system(target_cs_name) + if target_cs is None: + new_target_cs = CoordSystem( + name=target_cs_name, axes=expected_output_axes_union, virtual=make_target_cs_virtual_when_creating + ) + self._graph.add_coordinate_system(new_target_cs) + else: + new_target_cs_axes = get_union_of_axes(target_cs.axes + expected_output_axes_union) + created_outputs = set(new_target_cs_axes) - set(target_cs.axes) + if created_outputs: + new_target_cs = CoordSystem( + name=target_cs_name, + axes=new_target_cs_axes, + virtual=target_cs.virtual, + # same "bag" of coordinate systems, so inherits this attribute + ) + self._graph.add_coordinate_system(new_target_cs) + self._graph.add_transformation_edge( + ProjectAxisEdge( + input=target_cs, output=new_target_cs, created_outputs=created_outputs, dropped_inputs=set() + ) + ) + else: + new_target_cs = target_cs + + for ind, (_element_name, transformation) in enumerate(element_name_transformation_config.items()): + created_outputs = set(new_target_cs.axes) - set(intermediate_edge_list[ind].output.axes) + + if created_outputs: + self._graph.add_coordinate_system(intermediate_edge_list[ind].output) + self._graph.add_transformation_edge(intermediate_edge_list[ind]) + # create intermediate virtual cs, create and add edge between element and this virtual cs + + project_axis_edge = ProjectAxisEdge( + input=intermediate_edge_list[ind].output, + output=new_target_cs, + created_outputs=created_outputs, + dropped_inputs=set(), + ) + self._graph.add_transformation_edge(project_axis_edge) + # add an extra projection edge from virtual cs to new_target_cs + else: + transformation_edge = transformation._to_ngff_transformation_edge_existing_output_cs( + input_coordinate_system=intermediate_edge_list[ind].input, + output_coordinate_system=new_target_cs, + ) + self._graph.add_transformation_edge(transformation_edge) + # create and add edge based on element_name and transformation to new_target_cs + + def remove_transformations(self, transformation_edges: dict[str, BaseTransformationEdge]) -> None: + """ + Remove specific transformation edges from the graph. + + Parameters + ---------- + transformation_edges + dictionary with edge keys as keys and edges as values; those to be removed + """ + for edge_key, transformation_edge in transformation_edges.items(): + self._graph.remove_specific_transformation(edge_key, transformation_edge) + + def remove_all_outgoing_transformations_of_element(self, element_name: str) -> None: + """ + Remove all outgoing transformation edges of an element. + + Parameters + ---------- + element_name + The name of the element whose transformations should be removed. + """ + all_outgoing_edges_from_element = self._graph.get_outgoing_edges(element_name) + self.remove_transformations(all_outgoing_edges_from_element) + + def remove_all_incoming_transformations_of_element(self, element_name: str) -> None: + """ + Remove all incoming transformation edges of an element. + + Parameters + ---------- + element_name + The name of the element whose transformations should be removed. + """ + all_incoming_edges_from_element = self._graph.get_incoming_edges(element_name) + self.remove_transformations(all_incoming_edges_from_element) + + def remove_all_transformations_to_coordinate_system(self, coordinate_system_name: str) -> None: + """ + Remove all transformation edges that target a specific coordinate system. + + Parameters + ---------- + coordinate_system_name + The name of the target coordinate system. + + Raises + ------ + CoordinateSystemNotFoundError + If the coordinate system is not found. + AssertionError + If multiple coordinate systems with the same name exist (should never happen). + """ + cs_probables = [x for x in self._graph.list_coordinate_systems() if x.name == coordinate_system_name] + if len(cs_probables) == 1: + for predecessor in self._graph._graph.predecessors(cs_probables[0]): + self._graph.remove_all_transformations_between_coordinate_systems(predecessor, cs_probables[0]) + elif len(cs_probables) == 0: + CoordinateSystemNotFoundError(name=coordinate_system_name) + else: + raise AssertionError( + "This case should never happen, please raise an issue at https://github.com/scverse/spatialdata/issues" + ) diff --git a/src/spatialdata/_io/_utils.py b/src/spatialdata/_io/_utils.py index ba10405bd..57baeb9b4 100644 --- a/src/spatialdata/_io/_utils.py +++ b/src/spatialdata/_io/_utils.py @@ -30,14 +30,11 @@ from spatialdata._io.format import RasterFormatType, RasterFormatV01, RasterFormatV02, RasterFormatV03 from spatialdata._logging import logger from spatialdata._utils import get_pyramid_levels -from spatialdata.models._utils import ( - MappingToCoordinateSystem_t, - SpatialElement, - ValidAxis_t, - _validate_mapping_to_coordinate_system_type, -) +from spatialdata.models import SpatialElement +from spatialdata.transformations._graph.types import ValidAxis_t +from spatialdata.transformations._utils import MappingToCoordinateSystem_t, _validate_mapping_to_coordinate_system_type from spatialdata.transformations.ngff.ngff_transformations import NgffBaseTransformation -from spatialdata.transformations.transformations import BaseTransformation, _get_current_output_axes +from spatialdata.transformations.transformations import BaseTransformation def _get_transformations_from_ngff_dict( @@ -70,7 +67,7 @@ def overwrite_coordinate_transformations_non_raster( _validate_mapping_to_coordinate_system_type(transformations) ngff_transformations = [] for target_coordinate_system, t in transformations.items(): - output_axes = _get_current_output_axes(transformation=t, input_axes=tuple(axes)) + output_axes = t._get_resulting_output_axes(input_axes=tuple(axes)) ngff_transformations.append( t.to_ngff( input_axes=tuple(axes), @@ -112,7 +109,7 @@ def overwrite_coordinate_transformations_raster( # prepare the transformations in the dict representation ngff_transformations = [] for target_coordinate_system, t in transformations.items(): - output_axes = _get_current_output_axes(transformation=t, input_axes=tuple(axes)) + output_axes = t._get_resulting_output_axes(input_axes=tuple(axes)) ngff_transformations.append( t.to_ngff( input_axes=tuple(axes), diff --git a/src/spatialdata/_io/format.py b/src/spatialdata/_io/format.py index a9aa21ad3..ba522bc71 100644 --- a/src/spatialdata/_io/format.py +++ b/src/spatialdata/_io/format.py @@ -15,7 +15,7 @@ ) from shapely import GeometryType -from spatialdata.models.models import ATTRS_KEY, PointsModel, ShapesModel +from spatialdata.models import ATTRS_KEY, PointsModel, ShapesModel CoordinateTransform_t = list[dict[str, Any]] diff --git a/src/spatialdata/_io/io_raster.py b/src/spatialdata/_io/io_raster.py index 3806fabf7..b04ab2c2b 100644 --- a/src/spatialdata/_io/io_raster.py +++ b/src/spatialdata/_io/io_raster.py @@ -3,10 +3,13 @@ import warnings from collections.abc import Mapping, Sequence from pathlib import Path -from typing import Any, Literal, TypeGuard +from typing import Any, Literal, TypeGuard, get_args import dask.array as da import numpy as np +import ome_zarr as oz +import ome_zarr_models.v06.coordinate_transforms as ozm06trans +import xarray as xr import zarr from ome_zarr.format import Format from ome_zarr.io import ZarrLocation @@ -18,6 +21,7 @@ from ome_zarr.writer import write_multiscale as write_multiscale_ngff from ome_zarr.writer import write_multiscale_labels as write_multiscale_labels_ngff from xarray import DataArray, DataTree +from xarray.indexes import RangeIndex from spatialdata._io._utils import ( _get_transformations_from_ngff_dict, @@ -30,9 +34,11 @@ RasterFormatType, get_ome_zarr_format, ) +from spatialdata._types import ELEMENT_TYPE, ELEMENT_TYPE_RASTER from spatialdata._utils import get_pyramid_levels -from spatialdata.models.models import ATTRS_KEY +from spatialdata.models.literals import ATTRS_KEY from spatialdata.models.pyramids_utils import dask_arrays_to_datatree +from spatialdata.transformations._graph import Axis, BaseTransformationEdge, CoordSystem, parse_ngff_transf from spatialdata.transformations._utils import ( _get_transformations, _get_transformations_xarray, @@ -161,11 +167,104 @@ def _prepare_storage_options( return prepared_options +def try_read_ngff06_multiscale(store: Path) -> tuple[DataTree, Sequence[BaseTransformationEdge]]: + multiscale = oz.OMEZarrMultiscale.from_ome_zarr(str(store)) + assert isinstance(multiscale, oz.OMEZarrMultiscale) # disambiguate from OMEZarrLabel + return try_parse_ngff06_multiscale(multiscale) + + +def try_parse_ngff06_multiscale(multiscale: oz.OMEZarrMultiscale) -> tuple[DataTree, Sequence[BaseTransformationEdge]]: + """Parse an OMEZarMultiscale into a DataTree and collects Multiscale-level transforms.""" + name_to_cs: dict[str, CoordSystem] = {} + for cs in multiscale.metadata.coordinateSystems or (): + parsed_cs = CoordSystem.try_from_model(cs) + name_to_cs[cs.name] = parsed_cs + + parsed_transfs: list[BaseTransformationEdge] = [] + for transf in multiscale.metadata.coordinateTransformations or (): + in_cs_id = transf.input + out_cs_ref = transf.output + # these should not be None as per the spec + assert in_cs_id is not None + assert out_cs_ref is not None + + # FIXME: not handling references into labels yet, which use name and path + in_cs_name = in_cs_id.name + out_cs_name = out_cs_ref.name + assert in_cs_name is not None + assert out_cs_name is not None + + # assume CS references are valid via ome-zarr(-models)-py + input = name_to_cs[in_cs_name] + output = name_to_cs[out_cs_name] + parsed = parse_ngff_transf(input=input, output=output, model=transf) + parsed_transfs.append(parsed) + + omero = multiscale.omero + channel_names = None if omero is None else [d.color for d in omero.channels] + + data_tree = xr.DataTree() + for scale_idx, (ds_md, ds) in enumerate(zip(multiscale.metadata.datasets, multiscale.images, strict=True)): + transf = ds_md.coordinateTransformations[0] + + intrinsic_cs = name_to_cs[multiscale.metadata.intrinsic_coordinate_system.name] + assert transf.input is not None + assert transf.input.path is not None + + # This coord system doesn't exist explicitly in the NGFF file, nor will + # it exist in our graph of transformations; It is only created here + # for the sake of creating the transformations that will be expressed + # in levels of a xr.DataTree + pixel_cs = CoordSystem( + name=transf.input.path, + axes=tuple(Axis(name=ax.name, type=ax.type) for ax in intrinsic_cs.axes), + virtual=True, + ) + + pixel_cs_to_intrinsic_ngff = ozm06trans.Sequence(transformations=ds_md.coordinateTransformations) + seq = parse_ngff_transf(input=pixel_cs, output=intrinsic_cs, model=pixel_cs_to_intrinsic_ngff) + ds_shape = np.asarray(ds.data.shape) + transformed_start = seq.transform_points(np.zeros_like(ds_shape)[np.newaxis, :])[0] + transformed_stop = seq.transform_points((ds_shape - 1)[np.newaxis, :])[0] + + coords: xr.Coordinates = xr.Coordinates() + for low, high, ax, extent in zip(transformed_start, transformed_stop, intrinsic_cs.axes, ds_shape, strict=True): + if ax.type == "channel" and channel_names is not None: + coords = coords.merge({ax.name: channel_names}).coords + else: + axis_index = xr.Coordinates.from_xindex( + RangeIndex.linspace( + start=low, + stop=high, + num=extent, + endpoint=True, + dim=ax.name, + ) + ) + coords = coords.merge(axis_index).coords + + # Note: the magic "image" and "scale " strings mimic the current + # behavior from `dask_arrays_to_datatree` + data_tree[f"scale{scale_idx}"] = xr.Dataset( + { + "image": xr.DataArray( + ds.data, + name="image", + dims=intrinsic_cs.axes_names, + coords=coords, + ) + }, + ) + return data_tree, parsed_transfs + + def _read_multiscale( - store: str | Path, raster_type: Literal["image", "labels"], reader_format: Format + store: str | Path, raster_type: ELEMENT_TYPE_RASTER, reader_format: Format ) -> DataArray | DataTree: assert isinstance(store, str | Path) - assert raster_type in ["image", "labels"] + assert raster_type in get_args(ELEMENT_TYPE_RASTER), ValueError( + f"{raster_type} is not a valid raster type. Must be one of {get_args(ELEMENT_TYPE_RASTER)}." + ) nodes: list[Node] = [] image_loc = ZarrLocation(store, fmt=reader_format) @@ -214,7 +313,7 @@ def _read_multiscale( transformations = _get_transformations_from_ngff_dict(encoded_ngff_transformations) # if image, read channels metadata channels: list[Any] | None = None - if raster_type == "image": + if raster_type == ELEMENT_TYPE.IMAGE: if legacy_channels_metadata is not None: channels = [d["label"] for d in legacy_channels_metadata["channels"]] if omero_metadata is not None: @@ -229,7 +328,7 @@ def _read_multiscale( data = loaded_node.array(resolution=datasets[0]) si = DataArray( data, - name="image", + name=ELEMENT_TYPE.IMAGE, dims=axes, coords={"c": channels} if channels is not None else {}, ) @@ -264,9 +363,7 @@ def _get_multiscale_nodes(image_nodes: list[Node], nodes: list[Node]) -> list[No return nodes -def _get_raster_element_group( - raster_type: Literal["image", "labels"], group: zarr.Group, element_name: str -) -> zarr.Group: +def _get_raster_element_group(raster_type: ELEMENT_TYPE_RASTER, group: zarr.Group, element_name: str) -> zarr.Group: """Get the Zarr group holding a raster element that has just been written. Labels are nested one level deeper than images: ome-zarr writes them inside a "labels" group, so for them the @@ -285,7 +382,7 @@ def _get_raster_element_group( ------- The Zarr group of the raster element. """ - if raster_type != "labels": + if raster_type != ELEMENT_TYPE.LABELS: return group labels_group = group["labels"] if not isinstance(labels_group, zarr.Group): @@ -299,7 +396,7 @@ def _get_raster_element_group( def _write_raster( - raster_type: Literal["image", "labels"], + raster_type: ELEMENT_TYPE_RASTER, raster_data: DataArray | DataTree, group: zarr.Group, name: str, @@ -332,12 +429,13 @@ def _write_raster( metadata Additional metadata for the raster element """ - if raster_type not in ["image", "labels"]: - raise ValueError(f"{raster_type} is not a valid raster type. Must be 'image' or 'labels'.") + assert raster_type in get_args(ELEMENT_TYPE_RASTER), ValueError( + f"{raster_type} is not a valid raster type. Must be one of {get_args(ELEMENT_TYPE_RASTER)}." + ) # "name" and "label_metadata" are only used for labels. "name" is written in write_multiscale_ngff() but ignored in # write_image_ngff() (possibly an ome-zarr-py bug). We only use "name" to ensure correct group access in the # ome-zarr API. - if raster_type == "labels": + if raster_type == ELEMENT_TYPE.LABELS: metadata["name"] = name metadata["label_metadata"] = label_metadata @@ -367,7 +465,7 @@ def _write_raster( raise ValueError("Not a valid labels object") group = _get_raster_element_group(raster_type, group, name) - if raster_type == "image": + if raster_type == ELEMENT_TYPE.IMAGE: # ome-zarr-py >= 0.18 no longer writes the omero channel metadata, so we write it ourselves. overwrite_channel_names(group, raster_data) if ATTRS_KEY not in group.attrs: @@ -463,7 +561,7 @@ def _update_dict_v3(d: dict[str, Any]) -> None: def _write_raster_dataarray( - raster_type: Literal["image", "labels"], + raster_type: ELEMENT_TYPE_RASTER, group: zarr.Group, element_name: str, raster_data: DataArray, @@ -493,7 +591,7 @@ def _write_raster_dataarray( metadata Additional metadata for the raster element """ - write_single_scale_ngff = write_image_ngff if raster_type == "image" else write_labels_ngff + write_single_scale_ngff = write_image_ngff if raster_type == ELEMENT_TYPE.IMAGE else write_labels_ngff data = raster_data.data transformations = _get_transformations(raster_data) @@ -532,7 +630,7 @@ def _write_raster_dataarray( def _write_raster_datatree( - raster_type: Literal["image", "labels"], + raster_type: ELEMENT_TYPE_RASTER, group: zarr.Group, element_name: str, raster_data: DataTree, @@ -562,7 +660,9 @@ def _write_raster_datatree( metadata Additional metadata for the raster element """ - write_multi_scale_ngff = write_multiscale_ngff if raster_type == "image" else write_multiscale_labels_ngff + write_multi_scale_ngff = ( + write_multiscale_ngff if raster_type == ELEMENT_TYPE.IMAGE else write_multiscale_labels_ngff + ) data = get_pyramid_levels(raster_data, attr="data") list_of_input_axes: list[Any] = get_pyramid_levels(raster_data, attr="dims") assert len(set(list_of_input_axes)) == 1 @@ -632,7 +732,7 @@ def write_image( ) _write_raster( - raster_type="image", + raster_type=ELEMENT_TYPE.IMAGE, raster_data=image, group=group, name=name, @@ -659,7 +759,7 @@ def write_labels( ) _write_raster( - raster_type="labels", + raster_type=ELEMENT_TYPE.LABELS, raster_data=labels, group=group, name=name, diff --git a/src/spatialdata/_io/io_shapes.py b/src/spatialdata/_io/io_shapes.py index 833406d39..787b72973 100644 --- a/src/spatialdata/_io/io_shapes.py +++ b/src/spatialdata/_io/io_shapes.py @@ -26,6 +26,7 @@ ShapesFormatV03, _parse_version, ) +from spatialdata.literals import TRANSFORM_KEY from spatialdata.models import ShapesModel, get_axes_names from spatialdata.transformations._utils import ( _get_transformations, @@ -187,8 +188,6 @@ def _write_shapes_v02_v03( Whether to use the WKB or geoarrow encoding for GeoParquet. See :meth:`geopandas.GeoDataFrame.to_parquet` for details. """ - from spatialdata.models._utils import TRANSFORM_KEY - element_store = group.store_path.store if not isinstance(element_store, LocalStore): raise TypeError(f"Writing a shapes element requires a local zarr store, got {type(element_store).__name__}.") diff --git a/src/spatialdata/_io/io_zarr.py b/src/spatialdata/_io/io_zarr.py index 186d7c974..cee3aba7d 100644 --- a/src/spatialdata/_io/io_zarr.py +++ b/src/spatialdata/_io/io_zarr.py @@ -5,7 +5,7 @@ from collections.abc import Callable from json import JSONDecodeError from pathlib import Path -from typing import Any, Literal, cast +from typing import Any, Literal, cast, get_args import zarr.storage from anndata import AnnData @@ -28,7 +28,14 @@ from spatialdata._io.io_shapes import _read_shapes from spatialdata._io.io_table import _read_table from spatialdata._logging import logger -from spatialdata._types import JSONValue, Raster_T +from spatialdata._types import ( + ELEMENT_TYPE, + ELEMENT_TYPE_RASTER, + ELEMENT_TYPE_VECTOR, + GROUP_NAME, + JSONValue, + Raster_T, +) def _read_zarr_group_spatialdata_element( @@ -37,8 +44,8 @@ def _read_zarr_group_spatialdata_element( sdata_version: Literal["0.1", "0.2"], selector: set[str], read_func: Callable[..., Any], - group_name: Literal["images", "labels", "shapes", "points", "tables"], - element_type: Literal["image", "labels", "shapes", "points", "tables"], + group_name: GROUP_NAME, + element_type: ELEMENT_TYPE, element_container: (dict[str, Raster_T] | dict[str, DaskDataFrame] | dict[str, GeoDataFrame] | dict[str, AnnData]), on_bad_files: Literal[BadFileHandleMethod.ERROR, BadFileHandleMethod.WARN], ) -> None: @@ -76,14 +83,14 @@ def _read_zarr_group_spatialdata_element( ValueError, ), ): - if element_type in ["image", "labels"]: + if element_type in get_args(ELEMENT_TYPE_RASTER): reader_format = get_raster_format_for_read(elem_group, sdata_version) element = read_func( elem_group_path, element_type, reader_format, ) - elif element_type in ["shapes", "points", "tables"]: + elif element_type in get_args(ELEMENT_TYPE_VECTOR) + (ELEMENT_TYPE.TABLES,): element = read_func(elem_group_path) else: raise ValueError(f"Unknown element type {element_type}") @@ -204,19 +211,19 @@ def read_zarr( # we could make this more readable. One can get lost when looking at this dict and iteration over the items group_readers: dict[ - Literal["images", "labels", "shapes", "points", "tables"], + GROUP_NAME, tuple[ Callable[..., Any], - Literal["image", "labels", "shapes", "points", "tables"], + ELEMENT_TYPE, dict[str, Raster_T] | dict[str, DaskDataFrame] | dict[str, GeoDataFrame] | dict[str, AnnData], ], ] = { - # ome-zarr-py needs a kwargs that has "image" has key. So here we have "image" and not "images" - "images": (_read_multiscale, "image", images), - "labels": (_read_multiscale, "labels", labels), - "points": (_read_points, "points", points), - "shapes": (_read_shapes, "shapes", shapes), - "tables": (_read_table, "tables", tables), + # ome-zarr-py needs a kwargs that has "image" as key. So here we have "image" and not "images" + GROUP_NAME.IMAGES: (_read_multiscale, ELEMENT_TYPE.IMAGE, images), + GROUP_NAME.LABELS: (_read_multiscale, ELEMENT_TYPE.LABELS, labels), + GROUP_NAME.POINTS: (_read_points, ELEMENT_TYPE.POINTS, points), + GROUP_NAME.SHAPES: (_read_shapes, ELEMENT_TYPE.SHAPES, shapes), + GROUP_NAME.TABLES: (_read_table, ELEMENT_TYPE.TABLES, tables), } for group_name, ( read_func, @@ -295,7 +302,7 @@ def _get_groups_for_element( The Zarr groups for the root, element_type and element for a specific element. """ if not isinstance(zarr_path, Path): - raise ValueError("zarr_path should be a Path object") + raise TypeError("zarr_path should be a Path object") if element_type not in [ "images", diff --git a/src/spatialdata/_types.py b/src/spatialdata/_types.py index e41de9fe7..46c3f19c5 100644 --- a/src/spatialdata/_types.py +++ b/src/spatialdata/_types.py @@ -1,11 +1,22 @@ from __future__ import annotations -from typing import Any +from enum import StrEnum +from typing import Any, Literal import numpy as np from xarray import DataArray, DataTree -__all__ = ["ArrayLike", "ColorLike", "DTypeLike", "JSONValue", "Raster_T"] +__all__ = [ + "ELEMENT_TYPE", + "ELEMENT_TYPE_RASTER", + "ELEMENT_TYPE_VECTOR", + "GROUP_NAME", + "ArrayLike", + "ColorLike", + "DTypeLike", + "JSONValue", + "Raster_T", +] from numpy.typing import DTypeLike, NDArray @@ -21,8 +32,29 @@ type Raster_T = DataArray | DataTree ColorLike = tuple[float, ...] | str + # A value that survives a round-trip through JSON, which is the invariant that `SpatialData.attrs` must satisfy: the # attrs are persisted with `zarr.Group.attrs.put()`, which rejects anything that is not JSON-serializable (e.g. numpy # arrays, sets, DataFrames). Note that JSON has no tuples and only string keys, so a tuple is read back as a list and # a non-string key as a string. type JSONValue = dict[str, JSONValue] | list[JSONValue] | str | int | float | bool | None + + +class ELEMENT_TYPE(StrEnum): + IMAGE = "image" + LABELS = "labels" + SHAPES = "shapes" + POINTS = "points" + TABLES = "tables" + + +class GROUP_NAME(StrEnum): + IMAGES = "images" + LABELS = "labels" + SHAPES = "shapes" + POINTS = "points" + TABLES = "tables" + + +ELEMENT_TYPE_RASTER = Literal[ELEMENT_TYPE.IMAGE, ELEMENT_TYPE.LABELS] +ELEMENT_TYPE_VECTOR = Literal[ELEMENT_TYPE.POINTS, ELEMENT_TYPE.SHAPES] diff --git a/src/spatialdata/_utils.py b/src/spatialdata/_utils.py index a0086b7b2..b0ff0b7fd 100644 --- a/src/spatialdata/_utils.py +++ b/src/spatialdata/_utils.py @@ -13,11 +13,17 @@ from anndata import AnnData from dask import array as da from dask import config -from dask.array import Array as DaskArray from xarray import DataArray, Dataset, DataTree from spatialdata._types import ArrayLike, ListOrNDArrayFloating -from spatialdata.transformations import Sequence, Translation, get_transformation, set_transformation +from spatialdata.models import get_axes_names +from spatialdata.transformations import ( + Sequence, + Translation, + compute_coordinates, + get_transformation, + set_transformation, +) RT = TypeVar("RT") @@ -79,8 +85,6 @@ def unpad_raster(raster: DataArray | DataTree) -> DataArray | DataTree: ------- The unpadded raster. """ - from spatialdata.models import get_axes_names - from spatialdata.transformations._utils import compute_coordinates def _compute_paddings(data: DataArray, axis: str) -> tuple[int, int]: others = list(data.dims) @@ -330,37 +334,3 @@ def _error_message_add_element() -> None: "write_labels(), write_points(), write_shapes() and write_table(). We are going to make these calls more " "ergonomic in a follow up PR." ) - - -def _check_match_length_channels_c_dim( - data: DaskArray | DataArray | DataTree, c_coords: str | list[str], dims: tuple[str, ...] -) -> list[str]: - """ - Check whether channel names `c_coords` are of equal length to the `c` dimension of the data. - - Parameters - ---------- - data - The image array - c_coords - The channel names - dims - The axes names in the order that is the same as the `ImageModel` from which it is derived. - - Returns - ------- - c_coords - The channel names as list - """ - c_index = dims.index("c") - c_length = ( - data.shape[c_index] if isinstance(data, DataArray | DaskArray) else data["scale0"]["image"].shape[c_index] - ) - if isinstance(c_coords, str): - c_coords = [c_coords] - if c_coords is not None and len(c_coords) != c_length: - raise ValueError( - f"The number of channel names `{len(c_coords)}` does not match the length of dimension 'c'" - f" with length {c_length}." - ) - return c_coords diff --git a/src/spatialdata/literals/__init__.py b/src/spatialdata/literals/__init__.py new file mode 100644 index 000000000..78041a78a --- /dev/null +++ b/src/spatialdata/literals/__init__.py @@ -0,0 +1,14 @@ +from __future__ import annotations + +from spatialdata.literals.graph_literals import ( + DEFAULT_COORDINATE_SYSTEM, + TRANSFORM_KEY, + C, + NgffAxisType, + X, + Y, + Z, + axis_type_mapping_ngff, +) + +__all__ = ["DEFAULT_COORDINATE_SYSTEM", "C", "X", "Y", "Z", "NgffAxisType", "axis_type_mapping_ngff", "TRANSFORM_KEY"] diff --git a/src/spatialdata/literals/graph_literals.py b/src/spatialdata/literals/graph_literals.py new file mode 100644 index 000000000..05bacdbb6 --- /dev/null +++ b/src/spatialdata/literals/graph_literals.py @@ -0,0 +1,26 @@ +from __future__ import annotations + +from enum import StrEnum + +"""Key used to store transformation data in graph edge attributes.""" +DEFAULT_COORDINATE_SYSTEM = "global" +C = "c" +Z = "z" +Y = "y" +X = "x" + + +class NgffAxisType(StrEnum): + """Enum defining NGFF axis types""" + + CHANNEL = "channel" + SPACE = "space" + + +axis_type_mapping_ngff = { + C: NgffAxisType.CHANNEL, + X: NgffAxisType.SPACE, + Y: NgffAxisType.SPACE, + Z: NgffAxisType.SPACE, +} +TRANSFORM_KEY = "transform" diff --git a/src/spatialdata/models/__init__.py b/src/spatialdata/models/__init__.py index 01db9dc2d..5023e3671 100644 --- a/src/spatialdata/models/__init__.py +++ b/src/spatialdata/models/__init__.py @@ -2,20 +2,17 @@ from spatialdata._core.validation import check_target_region_column_symmetry from spatialdata.models._utils import ( - C, - SpatialElement, - X, - Y, - Z, + convert_region_column_to_categorical, force_2d, get_axes_names, get_channel_names, - get_spatial_axes, + get_raster_model_from_data_dims, points_dask_dataframe_to_geopandas, points_geopandas_to_dask_dataframe, set_channel_names, - validate_axes, - validate_axis_name, +) +from spatialdata.models.literals import ( + ATTRS_KEY, ) from spatialdata.models.models import ( Image2DModel, @@ -23,12 +20,14 @@ Labels2DModel, Labels3DModel, PointsModel, + RasterSchema, ShapesModel, TableModel, get_model, get_table_keys, validate_element, ) +from spatialdata.models.types import SpatialElement __all__ = [ "Labels2DModel", @@ -37,17 +36,10 @@ "Image3DModel", "ShapesModel", "PointsModel", + "RasterSchema", "TableModel", "get_model", "SpatialElement", - "get_spatial_axes", - "validate_axes", - "validate_axis_name", - "X", - "Y", - "Z", - "C", - "get_axes_names", "points_geopandas_to_dask_dataframe", "points_dask_dataframe_to_geopandas", "check_target_region_column_symmetry", @@ -56,4 +48,8 @@ "get_channel_names", "set_channel_names", "force_2d", + "ATTRS_KEY", + "get_axes_names", + "convert_region_column_to_categorical", + "get_raster_model_from_data_dims", ] diff --git a/src/spatialdata/models/_accessor.py b/src/spatialdata/models/_accessor.py index 77acf437c..7dec5acf4 100644 --- a/src/spatialdata/models/_accessor.py +++ b/src/spatialdata/models/_accessor.py @@ -138,7 +138,9 @@ def __repr__(self) -> str: ]: wrap_method_with_attrs(method_name=method_name, dask_class=DaskSeries) -wrap_indexer_with_attrs(indexer_name="loc", dask_class=DaskDataFrame) -wrap_indexer_with_attrs(indexer_name="iloc", dask_class=DaskDataFrame) -wrap_indexer_with_attrs(indexer_name="loc", dask_class=DaskSeries) -# DaskSeries do not have iloc + +def register_attrs_accessor_dask_df_series() -> None: + wrap_indexer_with_attrs(indexer_name="loc", dask_class=DaskDataFrame) + wrap_indexer_with_attrs(indexer_name="iloc", dask_class=DaskDataFrame) + wrap_indexer_with_attrs(indexer_name="loc", dask_class=DaskSeries) + # DaskSeries do not have iloc diff --git a/src/spatialdata/models/_utils.py b/src/spatialdata/models/_utils.py index 972e2d9be..188c96676 100644 --- a/src/spatialdata/models/_utils.py +++ b/src/spatialdata/models/_utils.py @@ -8,6 +8,7 @@ import geopandas import numpy as np import pandas as pd +from dask.array import Array as DaskArray from dask.dataframe import DataFrame as DaskDataFrame from geopandas import GeoDataFrame from shapely import get_coordinate_dimension @@ -16,23 +17,13 @@ from xarray import DataArray, DataTree from spatialdata._logging import logger -from spatialdata._utils import _check_match_length_channels_c_dim -from spatialdata.transformations.transformations import BaseTransformation - -type SpatialElement = DataArray | DataTree | GeoDataFrame | DaskDataFrame -TRANSFORM_KEY = "transform" -DEFAULT_COORDINATE_SYSTEM = "global" -ValidAxis_t = str -MappingToCoordinateSystem_t = dict[str, BaseTransformation] -C = "c" -Z = "z" -Y = "y" -X = "x" +from spatialdata.literals import C, X, Y, Z if TYPE_CHECKING: from anndata import AnnData from spatialdata.models.models import RasterSchema + from spatialdata.models.types import SpatialElement # mypy says that we can't do isinstance(something, SpatialElement), @@ -56,143 +47,6 @@ def has_type_spatial_element(e: Any) -> bool: return isinstance(e, DataArray | DataTree | GeoDataFrame | DaskDataFrame) -# added this code as part of a refactoring to catch errors earlier -def _validate_mapping_to_coordinate_system_type(transformations: MappingToCoordinateSystem_t | None) -> None: - if not ( - transformations is None - or isinstance(transformations, dict) - and all(isinstance(k, str) and isinstance(v, BaseTransformation) for k, v in transformations.items()) - ): - raise TypeError( - f"Transform must be of type {MappingToCoordinateSystem_t} or None, but is of type {type(transformations)}." - ) - - -def validate_axis_name(axis: ValidAxis_t) -> None: - """ - Check if the axis name is valid. - - Parameters - ---------- - axis - The axis name - - Raises - ------ - TypeError - If the axis name not in ["c", "x", "y", "z"] - """ - if axis not in ["c", "x", "y", "z"]: - raise TypeError(f"Invalid axis: {axis}") - - -def validate_axes(axes: tuple[ValidAxis_t, ...]) -> None: - """ - Check if the names of the axes are valid. - - Parameters - ---------- - axes - The names of the axes - - Raises - ------ - TypeError - If not in ["c", "x", "y", "z"] - """ - for ax in axes: - validate_axis_name(ax) - if len(axes) != len(set(axes)): - raise ValueError("Axes must be unique.") - - -def get_spatial_axes(axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: - """ - Get the spatial axes of interest. - - Parameters - ---------- - axes - Should be a subset of ['x', 'y', 'z', 'c'] - - Returns - ------- - The spatial axes, i.e. the input axes but without 'c' - """ - validate_axes(axes) - return tuple(ax for ax in axes if ax in [X, Y, Z]) - - -@singledispatch -def get_axes_names(e: SpatialElement | pd.DataFrame) -> tuple[str, ...]: - """ - Get the dimensions of a SpatialElement. - - Parameters - ---------- - e - SpatialElement - - Returns - ------- - Dimensions of the SpatialElement (e.g. ("z", "y", "x")) - """ - raise TypeError(f"Unsupported type: {type(e)}") - - -@get_axes_names.register(DataArray) -def _(e: DataArray) -> tuple[str, ...]: - dims = tuple(str(dim) for dim in e.dims) - _validate_dims(dims) - return dims - - -@get_axes_names.register(DataTree) -def _(e: DataTree) -> tuple[str, ...]: - if "scale0" in e: - # dims_coordinates = tuple(i for i in e["scale0"].dims.keys()) - - scale0 = e["scale0"] - if not isinstance(scale0, DataTree): - raise TypeError(f"Expected `scale0` to be a DataTree node, got {type(scale0).__name__}.") - variables = list(scale0.values()) - assert len(variables) == 1 - xdata = variables[0] - dims_data = tuple(str(dim) for dim in xdata.dims) - - _validate_dims(dims_data) - return dims_data - raise ValueError("Spatialdata DataTree does not contain the scale0 key") - - -@get_axes_names.register(GeoDataFrame) -def _(e: GeoDataFrame) -> tuple[str, ...]: - all_dims = (X, Y, Z) - n = get_coordinate_dimension(e.geometry.iloc[0]) - dims = all_dims[:n] - if Z not in dims and Z in e.columns: - dims += (Z,) - _validate_dims(dims) - return dims - - -@get_axes_names.register(DaskDataFrame) -@get_axes_names.register(pd.DataFrame) -def _(e: DaskDataFrame) -> tuple[str, ...]: - valid_dims = (X, Y, Z) - dims = tuple([c for c in valid_dims if c in e.columns]) - _validate_dims(dims) - return dims - - -def _validate_dims(dims: tuple[str, ...]) -> None: - for c in dims: - if c not in (X, Y, Z, C): - raise ValueError(f"Invalid dimension: {c}") - if dims not in [(X,), (Y,), (Z,), (C,), (X, Y), (X, Y, Z), (Y, X), (Z, Y, X), (C, Y, X), (C, Z, Y, X)]: - raise ValueError(f"Invalid dimensions: {dims}") - - def points_dask_dataframe_to_geopandas(points: DaskDataFrame, suppress_z_warning: bool = False) -> GeoDataFrame: """ Convert a Dask DataFrame to a GeoDataFrame. @@ -364,7 +218,7 @@ def get_raster_model_from_data_dims(dims: tuple[str, ...]) -> type[RasterSchema] ------- The raster model corresponding to the dimensions of the data. """ - from spatialdata.models.models import Image2DModel, Image3DModel, Labels2DModel, Labels3DModel + from spatialdata.models import Image2DModel, Image3DModel, Labels2DModel, Labels3DModel if not set(dims).issubset({C, Z, Y, X}): raise ValueError(f"Invalid dimensions: {dims}") @@ -375,7 +229,7 @@ def get_raster_model_from_data_dims(dims: tuple[str, ...]) -> type[RasterSchema] def convert_region_column_to_categorical(table: AnnData) -> None: - from spatialdata.models.models import TableModel + from spatialdata.models import TableModel if TableModel.ATTRS_KEY in table.uns: region_key = table.uns[TableModel.ATTRS_KEY][TableModel.REGION_KEY_KEY] @@ -434,3 +288,107 @@ def _get_uint_dtype(value: int) -> str: else: raise ValueError(f"Maximum cell number is {value}. Values higher than {max_uint64} are not supported.") return dtype + + +@singledispatch +def get_axes_names(e: SpatialElement | pd.DataFrame) -> tuple[str, ...]: + """ + Get the dimensions of a SpatialElement. + + Parameters + ---------- + e + SpatialElement + + Returns + ------- + Dimensions of the SpatialElement (e.g. ("z", "y", "x")) + """ + raise TypeError(f"Unsupported type: {type(e)}") + + +@get_axes_names.register(DataTree) +def _(e: DataTree) -> tuple[str, ...]: + if "scale0" in e: + # dims_coordinates = tuple(i for i in e["scale0"].dims.keys()) + + scale0 = e["scale0"] + if not isinstance(scale0, DataTree): + raise TypeError(f"Expected `scale0` to be a DataTree node, got {type(scale0).__name__}.") + variables = list(scale0.values()) + assert len(variables) == 1 + xdata = variables[0] + dims_data = tuple(str(dim) for dim in xdata.dims) + + _validate_dims(dims_data) + return dims_data + raise ValueError("Spatialdata DataTree does not contain the scale0 key") + + +@get_axes_names.register(GeoDataFrame) +def _(e: GeoDataFrame) -> tuple[str, ...]: + all_dims = (X, Y, Z) + n = get_coordinate_dimension(e.geometry.iloc[0]) + dims = all_dims[:n] + if Z not in dims and Z in e.columns: + dims += (Z,) + _validate_dims(dims) + return dims + + +@get_axes_names.register(DaskDataFrame) +@get_axes_names.register(pd.DataFrame) +def _(e: DaskDataFrame | pd.DataFrame) -> tuple[str, ...]: + valid_dims = (X, Y, Z) + dims = tuple([c for c in valid_dims if c in e.columns]) + _validate_dims(dims) + return dims + + +@get_axes_names.register(DataArray) +def _(e: DataArray) -> tuple[str, ...]: + dims = tuple(str(dim) for dim in e.dims) + _validate_dims(dims) + return dims + + +def _validate_dims(dims: tuple[str, ...]) -> None: + for c in dims: + if c not in (X, Y, Z, C): + raise ValueError(f"Invalid dimension: {c}") + if dims not in [(X,), (Y,), (Z,), (C,), (X, Y), (X, Y, Z), (Y, X), (Z, Y, X), (C, Y, X), (C, Z, Y, X)]: + raise ValueError(f"Invalid dimensions: {dims}") + + +def _check_match_length_channels_c_dim( + data: DaskArray | DataArray | DataTree, c_coords: str | list[str], dims: tuple[str, ...] +) -> list[str]: + """ + Check whether channel names `c_coords` are of equal length to the `c` dimension of the data. + + Parameters + ---------- + data + The image array + c_coords + The channel names + dims + The axes names in the order that is the same as the `ImageModel` from which it is derived. + + Returns + ------- + c_coords + The channel names as list + """ + c_index = dims.index("c") + c_length = ( + data.shape[c_index] if isinstance(data, DataArray | DaskArray) else data["scale0"]["image"].shape[c_index] + ) + if isinstance(c_coords, str): + c_coords = [c_coords] + if c_coords is not None and len(c_coords) != c_length: + raise ValueError( + f"The number of channel names `{len(c_coords)}` does not match the length of dimension 'c'" + f" with length {c_length}." + ) + return c_coords diff --git a/src/spatialdata/models/literals.py b/src/spatialdata/models/literals.py new file mode 100644 index 000000000..742cfcab5 --- /dev/null +++ b/src/spatialdata/models/literals.py @@ -0,0 +1,3 @@ +from __future__ import annotations + +ATTRS_KEY = "spatialdata_attrs" diff --git a/src/spatialdata/models/models.py b/src/spatialdata/models/models.py index b0a6c16f6..543f953a0 100644 --- a/src/spatialdata/models/models.py +++ b/src/spatialdata/models/models.py @@ -29,24 +29,21 @@ from spatialdata._core.validation import validate_table_attr_keys from spatialdata._logging import logger -from spatialdata._types import ArrayLike -from spatialdata._utils import _check_match_length_channels_c_dim from spatialdata.config import settings -from spatialdata.models import C, X, Y, Z, get_axes_names +from spatialdata.literals import DEFAULT_COORDINATE_SYSTEM, TRANSFORM_KEY, C, X, Y, Z from spatialdata.models._utils import ( - DEFAULT_COORDINATE_SYSTEM, - TRANSFORM_KEY, - MappingToCoordinateSystem_t, - SpatialElement, - _validate_mapping_to_coordinate_system_type, + _check_match_length_channels_c_dim, convert_region_column_to_categorical, + get_axes_names, ) from spatialdata.models.chunks_utils import Chunks_t +from spatialdata.models.literals import ATTRS_KEY from spatialdata.models.pyramids_utils import ScaleFactors_t # ozp -> ome-zarr-py from spatialdata.models.pyramids_utils import to_multiscale as to_multiscale_ozp from spatialdata.transformations._utils import ( _get_transformations, _set_transformations, + _validate_mapping_to_coordinate_system_type, compute_coordinates, ) from spatialdata.transformations.transformations import Identity @@ -54,9 +51,11 @@ if TYPE_CHECKING: from pandas._typing import DtypeObj -__all__ = ["Chunks_t", "ScaleFactors_t"] + from spatialdata._types import ArrayLike + from spatialdata.models.types import SpatialElement + from spatialdata.transformations._utils import MappingToCoordinateSystem_t -ATTRS_KEY = "spatialdata_attrs" +__all__ = ["Chunks_t", "ScaleFactors_t"] def _parse_transformations(element: SpatialElement, transformations: MappingToCoordinateSystem_t | None = None) -> None: diff --git a/src/spatialdata/models/types.py b/src/spatialdata/models/types.py new file mode 100644 index 000000000..d228c49a8 --- /dev/null +++ b/src/spatialdata/models/types.py @@ -0,0 +1,7 @@ +from __future__ import annotations + +from dask.dataframe import DataFrame as DaskDataFrame +from geopandas import GeoDataFrame +from xarray import DataArray, DataTree + +SpatialElement = DataArray | DataTree | GeoDataFrame | DaskDataFrame diff --git a/src/spatialdata/testing.py b/src/spatialdata/testing.py index f699d2e1c..6a5187e79 100644 --- a/src/spatialdata/testing.py +++ b/src/spatialdata/testing.py @@ -11,8 +11,7 @@ from spatialdata import SpatialData from spatialdata._core._elements import Images, Labels, Points, Shapes, Tables -from spatialdata.models import PointsModel -from spatialdata.models._utils import SpatialElement +from spatialdata.models import PointsModel, SpatialElement from spatialdata.transformations.operations import get_transformation _Elements = Images | Labels | Shapes | Points | Tables diff --git a/src/spatialdata/transformations/__init__.py b/src/spatialdata/transformations/__init__.py index e95a92766..c02c43f0a 100644 --- a/src/spatialdata/transformations/__init__.py +++ b/src/spatialdata/transformations/__init__.py @@ -1,5 +1,7 @@ from __future__ import annotations +from spatialdata.transformations._graph import get_spatial_axes +from spatialdata.transformations._utils import compute_coordinates from spatialdata.transformations.operations import ( align_elements_using_landmarks, get_transformation, @@ -34,4 +36,7 @@ "get_transformation_between_landmarks", "align_elements_using_landmarks", "remove_transformations_to_coordinate_system", + "compute_coordinates", + # utils + "get_spatial_axes", ] diff --git a/src/spatialdata/transformations/_graph/__init__.py b/src/spatialdata/transformations/_graph/__init__.py new file mode 100644 index 000000000..24ad0d46a --- /dev/null +++ b/src/spatialdata/transformations/_graph/__init__.py @@ -0,0 +1,118 @@ +from __future__ import annotations + +from spatialdata.transformations._graph.edges.exceptions import ( + AxisRedefinitionError, + DeterminantDifferentFromOne, + EmptyTransformSequenceError, + IncompatibleCoordSystemsError, + MissingAxisError, + NonZeroDroppedAxisError, + NotOrthonormalError, + UnexpectedShapeError, + UnmappedAxisError, +) +from spatialdata.transformations._graph.edges.ngff_aligned_edges import ( + AffineEdge, + BaseTransformationEdge, + ByDimensionEdge, + CsGen, + IdentityEdge, + MapAxisEdge, + ProjectAxisEdge, + RotationEdge, + ScaleEdge, + SequenceEdge, + TranslationEdge, + parse_affine, + parse_by_dimension, + parse_identity, + parse_map_axis, + parse_ngff_transf, + parse_project_axis, + parse_rotation, + parse_scale, + parse_sequence, + parse_translation, +) +from spatialdata.transformations._graph.exceptions import ( + AxisParsingException, + CannotRemoveCoordinateSystemError, + CoordinateSystemAlreadyExistsError, + CoordinateSystemNotFoundError, + ElementAlreadyExistsError, + ElementNotRegisteredToAnyCoordinateSystemError, + InvalidPathError, + NGFFCompatibilityError, + TransformationNotFoundError, + TransformationPathAmbiguousError, + TransformationPathAmbiguousMultipleEdgeExpectedError, + TransformationPathAmbiguousNoEdgeExpectedError, + TransformationPathNotFoundError, + TransformationPathNotSimple, +) +from spatialdata.transformations._graph.graph import TransformationGraph +from spatialdata.transformations._graph.utils import ( + get_instrinsic_coordinate_system, + get_spatial_axes, + get_union_of_axes, +) +from spatialdata.transformations._graph.vert import Axis, CoordSystem + +__all__ = [ + # from ngff aligned edge + "AffineEdge", + "BaseTransformationEdge", + "ByDimensionEdge", + "CsGen", + "IdentityEdge", + "MapAxisEdge", + "ProjectAxisEdge", + "RotationEdge", + "ScaleEdge", + "SequenceEdge", + "TranslationEdge", + "parse_affine", + "parse_by_dimension", + "parse_identity", + "parse_map_axis", + "parse_ngff_transf", + "parse_project_axis", + "parse_rotation", + "parse_scale", + "parse_sequence", + "parse_translation", + # from edge.exceptions + "NonZeroDroppedAxisError", + # from exceptions + "AxisParsingException", + "AxisRedefinitionError", + "CannotRemoveCoordinateSystemError", + "CoordinateSystemAlreadyExistsError", + "CoordinateSystemNotFoundError", + "DeterminantDifferentFromOne", + "ElementAlreadyExistsError", + "ElementNotRegisteredToAnyCoordinateSystemError", + "EmptyTransformSequenceError", + "IncompatibleCoordSystemsError", + "InvalidPathError", + "MissingAxisError", + "NGFFCompatibilityError", + "NotOrthonormalError", + "TransformationNotFoundError", + "TransformationPathAmbiguousError", + "TransformationPathAmbiguousMultipleEdgeExpectedError", + "TransformationPathAmbiguousNoEdgeExpectedError", + "TransformationPathNotFoundError", + "TransformationPathNotSimple", + "UnexpectedShapeError", + "UnmappedAxisError", + # from graph + "TransformationGraph", + # from utils + "get_instrinsic_coordinate_system", + "get_union_of_axes", + "get_spatial_axes", + # from vert + "Axis", + "CoordSystem", +] diff --git a/src/spatialdata/transformations/_graph/edges/__init__.py b/src/spatialdata/transformations/_graph/edges/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/spatialdata/transformations/_graph/edges/exceptions.py b/src/spatialdata/transformations/_graph/edges/exceptions.py new file mode 100644 index 000000000..88dbeb1b4 --- /dev/null +++ b/src/spatialdata/transformations/_graph/edges/exceptions.py @@ -0,0 +1,102 @@ +from __future__ import annotations + +import numpy as np + +from spatialdata._types import ArrayLike + + +class IncompatibleCoordSystemsError(Exception): + def __init__( + self, + transformation_name: str | None = None, + axes_error_message: str | None = None, + ) -> None: + self.transformation_name = transformation_name + axes_message = axes_error_message or "Output axes can't be mapped to input axes" + axes_transformation_message = ( + transformation_name and f"In the context of a `{transformation_name}` transformation, {axes_message}" + ) + super().__init__(axes_transformation_message or axes_message) + + +class MissingAxisError(Exception): + def __init__(self, axis_name: str, cs_name: str) -> None: + self.axis_name = axis_name + self.cs_name = cs_name + super().__init__(f"Axis {axis_name} is not in coordinate system {cs_name}") + + +class UnexpectedShapeError(Exception): + def __init__( + self, + *, + array_shape: tuple[int, ...], + expected_shape: tuple[int, ...] | str | None = None, + array_name: str | None = None, + ) -> None: + self.array_shape = array_shape + self.expected_shape = expected_shape + message = "Unexpected array shape" + if array_name is not None: + message += f"for '{array_name}'" + message += f": {array_shape}" + if expected_shape is not None: + message += f" instead of {expected_shape}" + super().__init__(message) + + +class DeterminantDifferentFromOne(Exception): + def __init__(self, matrix: ArrayLike) -> None: + self.matrix = matrix + super().__init__("Matrix does not have det(M) == 1") + + +class NotOrthonormalError(Exception): + def __init__(self, matrix: ArrayLike) -> None: + self.matrix = matrix + super().__init__("Matrix is not orthonormal") + + +class EmptyTransformSequenceError(Exception): + def __init__(self) -> None: + super().__init__("Empty sequence of transformations") + + +class AxisRedefinitionError(Exception): + def __init__(self, axis_name: str) -> None: + self.axis_name = axis_name + super().__init__(f"Axis {axis_name} is defined multiple times") + + +class UnmappedAxisError(Exception): + def __init__(self, axis_name: str, cs_name: str) -> None: + self.axis_name = axis_name + self.cs_name = cs_name + super().__init__(f"Axis {axis_name} from coordinate system {cs_name} is not mapped to anything") + + +class NonZeroDroppedAxisError(ValueError): + """ + Raised by :class:`DynamicDropAxisEdge` when a column that is supposed to be + dropped does not contain only zeros. + + Parameters + ---------- + axis_name : str + Name of the axis whose column violated the zero‑only requirement. + expected : float + The constant value that was expected (always ``0`` in the current + implementation). + observed : numpy.ndarray + The full column extracted from the input points that contains the + offending values. + """ + + def __init__(self, axis_name: str, expected: float, observed: np.ndarray) -> None: + self.axis_name = axis_name + self.expected = expected + self.observed = observed + super().__init__( + f"Axis '{axis_name}' contains non‑zero values; " + f"expected all entries to be {expected}, but observed {observed}" + ) diff --git a/src/spatialdata/transformations/_graph/edges/ngff_aligned_edges.py b/src/spatialdata/transformations/_graph/edges/ngff_aligned_edges.py new file mode 100644 index 000000000..aa1665ec4 --- /dev/null +++ b/src/spatialdata/transformations/_graph/edges/ngff_aligned_edges.py @@ -0,0 +1,1267 @@ +# pyright: strict + +from __future__ import annotations + +from abc import ABC, abstractmethod +from collections.abc import Mapping, Sequence +from typing import Final +from warnings import warn + +import numpy as np +import ome_zarr_models.v06.coordinate_transforms as ozm06trans + +from spatialdata._types import ArrayLike +from spatialdata.literals import NgffAxisType +from spatialdata.transformations._graph.edges.exceptions import ( + AxisRedefinitionError, + DeterminantDifferentFromOne, + EmptyTransformSequenceError, + IncompatibleCoordSystemsError, + MissingAxisError, + NotOrthonormalError, + UnexpectedShapeError, + UnmappedAxisError, +) +from spatialdata.transformations._graph.vert import Axis, CoordSystem + + +class BaseTransformationEdge(ABC): + """Base class for all the transformations defined by the NGFF specification.""" + + input: Final[CoordSystem] + output: Final[CoordSystem] + name: str | None + + def __init__( + self, + *, + name: str | None = None, + input: CoordSystem, + output: CoordSystem, + ) -> None: + self.input = input + self.output = output + self.name = name + super().__init__() + + def __repr__(self) -> str: + domain = ", ".join(self.input.axes_names) + codomain = ", ".join(self.output.axes_names) + return f"{type(self).__name__} ({domain} -> {codomain})" + + @abstractmethod + def inverse(self, name: str | None = None) -> BaseTransformationEdge | None: + """Return the inverse of the transformation if it exists""" + + @abstractmethod + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + + @abstractmethod + def to_affine(self, name: str | None = None) -> AffineEdge: + """Convert the transformation to an affine transformation, whenever the conversion can be made.""" + + def _validate_transform_points_shapes(self, points: ArrayLike) -> None: + """ + Validate if the shape of the points (coordinates to be transformed) are consistent with the input size of the + transformation. + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + input_size = len(self.input.axes) + if len(points.shape) != 2 or points.shape[1] != input_size: + raise UnexpectedShapeError( + array_name="points", + expected_shape=f"(, {self.input.num_axes})", + array_shape=points.shape, + ) + + # order of the composition: self is applied first, then the transformation passed as argument + def compose_with(self, transformation: BaseTransformationEdge, name: str | None) -> BaseTransformationEdge: + """ + Compose the transfomation object with another transformation + + Parameters + ---------- + transformation + The transformation to compose with. + + Returns + ------- + The compoesed transformation. + + Notes + ------- + Self is applied first, then the transformation passed as argument. + + Raises + ------ + IncompatibleCoordSystemsError + if this transformation's output coordinate system doesn't match `transformation`'s input + """ + return SequenceEdge(transformations=[self, transformation], name=name) + + +class AffineEdge(BaseTransformationEdge): + """The Affine transformation from the NGFF specification.""" + + linear: Final[ArrayLike] + translation: Final[ArrayLike] + affine: Final[ArrayLike] + + def __init__( + self, + *, + name: str | None = None, + linear: ArrayLike, + translation: ArrayLike | None = None, + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + name + A human readable name for this transformation + linear + The linear part of this transformation, i.e., the one that keeps + the origin in the same place. Shape must be (output.num_axes, input.num_axes) + translation y + The translation part of this transformation, of shape (output.num_axes,) + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + + Raises + ------ + UnexpectedShapeError + if `linear`'s shape isn't (output.num_axes, input.num_axes) or + if `translation`'s shape isn't (output.num_axes,) + """ + num_inputs = input.num_axes + num_outputs = output.num_axes + translation = np.zeros(num_outputs) if translation is None else translation + + expected_linear_shape = (num_outputs, num_inputs) + if linear.shape != expected_linear_shape: + raise UnexpectedShapeError( + array_name="linear", array_shape=linear.shape, expected_shape=expected_linear_shape + ) + expected_translation_shape = (num_outputs,) + if translation.shape != expected_translation_shape: + raise UnexpectedShapeError( + array_name="translation", array_shape=translation.shape, expected_shape=expected_translation_shape + ) + + self.linear = linear + self.translation = translation + + self.affine = np.zeros((num_outputs + 1, num_inputs + 1)) + self.affine[:-1, :-1] = self.linear + self.affine[:-1, -1] = self.translation + self.affine[-1, -1] = 1 + + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + s = super().__repr__() + "\n" + s += "\n".join(str(row) for row in self.affine) + return s + + @classmethod + def from_affine_matrix( + cls, + *, + name: str | None, + affine_matrix: ArrayLike, + input: CoordSystem, + output: CoordSystem, + ) -> AffineEdge: + """Creates an AffineEdge from a raw affine matrix in homogenous coordinates + + Parameters + ---------- + name + A human readable name for this transformation + affine_matrix + row-major, (output.num_axes + 1, input.num_axes + 1) matrix with: + - linear part at the top left + - translation as the rightmost column + - last row is [0, 0, ..., 0, 1] + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + + Raises + ------ + UnexpectedShapeError + if `affine_matrix`'s shape isn't (output.num_axes + 1, input.num_axes + 1) + """ + return AffineEdge( + linear=affine_matrix[:-1, :-1], + translation=affine_matrix[:-1, -1], + input=input, + output=output, + name=name, + ) + + @classmethod + def mapping(cls, input: CoordSystem, output: CoordSystem, name: str | None = None) -> AffineEdge: + """Create an AffineEdge that maps input axes to output axes of the same name.""" + linear: ArrayLike = np.zeros((output.num_axes, input.num_axes), dtype=float) + for i, des_axis in enumerate(output.axes): + for j, src_axis in enumerate(input.axes): + if src_axis.name == des_axis.name: # FIXME: compare the entire axis? + linear[i, j] = 1 + return AffineEdge( + linear=linear, + input=input, + output=output, + name=name, + ) + + def inverse(self, name: str | None = None) -> BaseTransformationEdge | None: + cond = np.linalg.cond(self.affine) + if cond > 1e10: + warn( + f"Inverting an affine matrix with large condition number ({cond:.2e})", + RuntimeWarning, + stacklevel=2, + ) + try: + # FIXME: I think there are more efficient/precise ways to invert a matrix + inv = np.linalg.inv(self.affine) + except np.linalg.LinAlgError: + return None + return AffineEdge( + linear=inv[:-1, :-1], + translation=inv[:-1, -1], + input=self.output, + output=self.input, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + return (points @ self.linear.T) + self.translation + + def to_affine(self, name: str | None = None) -> AffineEdge: + return AffineEdge( + input=self.input, output=self.output, linear=self.linear, translation=self.translation, name=name + ) + + +class IdentityEdge(BaseTransformationEdge): + """The Identity transformation from the NGFF specification.""" + + def __init__( + self, + name: str | None, + *, + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + name + A human readable name for this transformation + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + + Raises + ------ + IncompatibleCoordSystemsError + if `input` and `output` don't have the same number of axes + """ + if input.num_axes != output.num_axes: + raise IncompatibleCoordSystemsError( + axes_error_message="Axes must have the same number of dimensions", + transformation_name=self.__class__.__name__, + ) + super().__init__(input=input, output=output, name=name) + + def inverse(self, name: str | None = None) -> BaseTransformationEdge: + return IdentityEdge(input=self.output, output=self.input, name=name) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + return points + + def to_affine(self, name: str | None = None) -> AffineEdge: + return AffineEdge( + linear=np.identity(self.input.num_axes), + input=self.input, + output=self.output, + name=name, + ) + + +class MapAxisEdge(BaseTransformationEdge): + """The MapAxis transformation from the NGFF specification.""" + + def __init__( + self, + *, + name: str | None = None, + input: CoordSystem, + output: CoordSystem, + input_to_output: Mapping[Axis, Axis], + ) -> None: + """ + Parameters + ---------- + name + A human readable name for this transformation + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. Must + have the same number of axes as `input` + input_to_output + A mapping from + + Raises + ------ + IncompatibleCoordSystemsError + if `input` and `output` don't have the number of axes + MissingAxisError + if an axis specified among keys of `input_to_output` is not in `input` + or + if an axis specified among values of `input_to_output` is not in `output` + UnmappedAxisError + if an axis specified in `input` is not a key of `input_to_output` + or + if an axis specified in `output` is not a value of `input_to_output` + + """ + + if input.num_axes != output.num_axes: + raise IncompatibleCoordSystemsError(axes_error_message="Input and output must have the same number of axes") + + unmapped_inputs = set(input.axes) + unmapped_outputs = set(output.axes) + for inp_ax, out_ax in input_to_output.items(): + try: + unmapped_inputs.remove(inp_ax) + except KeyError as e: + raise MissingAxisError(axis_name=inp_ax.name, cs_name=input.name) from e + try: + unmapped_outputs.remove(out_ax) + except KeyError as e: + raise MissingAxisError(axis_name=out_ax.name, cs_name=output.name) from e + + for ax in unmapped_inputs: + raise UnmappedAxisError(axis_name=ax.name, cs_name=input.name) + for ax in unmapped_outputs: + raise UnmappedAxisError(axis_name=ax.name, cs_name=output.name) + + self.input_to_output = input_to_output + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + s = super().__repr__() + "\n" + s += "\n".join(f" {out.name} <- {inp.name}\n" for out, inp in self.input_to_output.items()) + return s + + def inverse(self, name: str | None = None) -> BaseTransformationEdge: + return MapAxisEdge( + input=self.output, + output=self.input, + input_to_output={out: inp for inp, out in self.input_to_output.items()}, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + output_to_input = {out_ax: in_ax for in_ax, out_ax in self.input_to_output.items()} + new_indices = [self.input.axes.index(output_to_input[out_ax]) for out_ax in self.output.axes] + assert len(new_indices) == self.output.num_axes + mapped = points[:, new_indices] + assert isinstance(mapped, np.ndarray) + return mapped + + def to_affine(self, name: str | None = None) -> AffineEdge: + linear = np.zeros((self.output.num_axes, self.input.num_axes)) + for inp_ax, out_ax in self.input_to_output.items(): + inp_idx = self.input.axes.index(inp_ax) + out_idx = self.output.axes.index(out_ax) + linear[out_idx, inp_idx] = 1 + return AffineEdge(linear=linear, input=self.input, output=self.output) + + +class ProjectAxisEdge(BaseTransformationEdge): + dropped_inputs: Final[set[Axis]] + created_outputs: Final[set[Axis]] + + def __init__( + self, + *, + name: str | None = None, + input: CoordSystem, + output: CoordSystem, + dropped_inputs: set[Axis], + created_outputs: set[Axis], + ) -> None: + """ + Parameters + ---------- + dropped_inputs + axes in `input` that will be dropped by this transformation + created_outputs + axes in `output` that will be set to 0 + + Raises + ------ + MissingAxisError + axis in `dropped_inputs` not in `input` + axis in `created_outputs` not in `output` + IncompatibleCoordSystemsError + when input can't be mapped to output given dropped_inputs and created_outputs + """ + + for axis in dropped_inputs: + if axis not in input.axes: + raise MissingAxisError(axis_name=axis.name, cs_name=input.name) + for axis in created_outputs: + if axis not in output.axes: + raise MissingAxisError(axis_name=axis.name, cs_name=output.name) + if input.num_axes - len(dropped_inputs) + len(created_outputs) != output.num_axes: + message = f"Can't map from {input} to {output}" + if dropped_inputs: + message += f" dropping {dropped_inputs}" + if created_outputs: + message += f" creating {created_outputs}" + raise IncompatibleCoordSystemsError( + axes_error_message=message, + ) + + self.dropped_inputs = set(dropped_inputs) + self.created_outputs = set(created_outputs) + super().__init__(name=name, input=input, output=output) + + def to_affine(self, name: str | None = None) -> AffineEdge: + + linear = np.zeros((self.output.num_axes, self.input.num_axes), dtype=float) + input_indices = iter(range(self.input.num_axes)) + for out_idx, out_ax in enumerate(self.output.axes): + if out_ax in self.created_outputs: + continue + in_idx = next(input_indices) + if in_idx in self.dropped_inputs: + continue + linear[out_idx, in_idx] = 1 + + return AffineEdge(name=name, input=self.input, output=self.output, linear=linear) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + return self.to_affine().transform_points(points) + + def inverse(self, name: str | None = None) -> BaseTransformationEdge | None: + # FIXME: there may be other cases where this is invertible + if self.input.num_axes != self.output.num_axes: + return None + if len(self.dropped_inputs) > 0: + return None + if len(self.created_outputs) > 0: + return None + return ProjectAxisEdge( + input=self.output, + output=self.input, + dropped_inputs=set(), + created_outputs=set(), + name=name, + ) + + +class TranslationEdge(BaseTransformationEdge): + """The Translation transformation from the NGFF specification.""" + + def __init__( + self, + *, + name: str | None = None, + translation: ArrayLike, + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + name + A human readable name for this transformation + translation + A vector of shape (input.num_axes,) specifying the translation along each axis. + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + + Raises + ------ + IncompatibleCoordSystemsError + If the input and output have different number of dimensions + UnexpectedShapeError + if `translation`'s shape doesn't match `input` + """ + if input.num_axes != output.num_axes: + raise IncompatibleCoordSystemsError(axes_error_message="Number of input and output axes must be the same") + expected_translation_shape = (input.num_axes,) + if translation.shape != expected_translation_shape: + raise UnexpectedShapeError( + array_name="translation", + array_shape=translation.shape, + expected_shape=expected_translation_shape, + ) + self.translation = translation + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + return super().__repr__() + str(self.translation) + + def inverse(self, name: str | None = None) -> BaseTransformationEdge: + return TranslationEdge( + translation=-self.translation, + input=self.output, + output=self.input, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + return points + self.translation + + def to_affine(self, name: str | None = None) -> AffineEdge: + return AffineEdge( + linear=np.identity(self.input.num_axes), + translation=self.translation, + input=self.input, + output=self.output, + name=name, + ) + + +class ScaleEdge(BaseTransformationEdge): + """The Scale transformation from the NGFF specification.""" + + def __init__( + self, + *, + name: str | None = None, + scale: ArrayLike, + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + scale + A vector specifying the scale along each axis of `input`. + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + + Raises + ------ + UnexpectedShapeError + If scale doesn't have the same number of elements as input has axes + IncompatibleCoordSystemsError + If input and output have different number of axes + """ + expected_scale_shape = (input.num_axes,) + if scale.shape != expected_scale_shape: + raise UnexpectedShapeError(array_name="scale", array_shape=scale.shape, expected_shape=expected_scale_shape) + if input.num_axes != output.num_axes: + raise IncompatibleCoordSystemsError( + axes_error_message="input and output must have same number of dimensions" + ) + self.scale = scale + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + return super().__repr__() + str(self.scale) + + def inverse(self, name: str | None = None) -> ScaleEdge | None: + if any(s == 0 for s in self.scale): + return None + new_scale = 1 / self.scale + return ScaleEdge( + scale=new_scale, + input=self.output, + output=self.input, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + return points * self.scale + + def to_affine(self, name: str | None = None) -> AffineEdge: + return AffineEdge( + linear=np.diag(self.scale), + input=self.input, + output=self.output, + name=name, + ) + + +class RotationEdge(BaseTransformationEdge): + """The Rotation transformation from the NGFF specification.""" + + rotation: Final[ArrayLike] + + def __init__( + self, + *, + name: str | None = None, + linear_matrix: ArrayLike, + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + linear_matrix + And orthonormal matrix of shape (output.num_axes, input.num_axes) + input + Input coordinate system of the transformation. + output + Output coordinate system of the transformation. + Raises + ------ + UnexpectedShapeError + if linear_matrix's shape isn't (output.num_axes, input.num_axes) + IncompatibleCoordSystemsError + if input and output don't have the same number of axes + DeterminantDifferentFromOne + if linear_matrix doesn't have determinant ~= 1 + NotOrthonormal + if linear_matrix is not orthonormal + """ + if input.num_axes != output.num_axes: + raise IncompatibleCoordSystemsError( + axes_error_message="input and output should have the same number of axes" + ) + expected_shape = (output.num_axes, input.num_axes) + if linear_matrix.shape != expected_shape: + raise UnexpectedShapeError( + array_name="linear_matrix", array_shape=linear_matrix.shape, expected_shape=expected_shape + ) + if not np.allclose(linear_matrix.T @ linear_matrix, np.identity(input.num_axes)): + raise NotOrthonormalError(matrix=linear_matrix) + if not np.isclose(np.linalg.det(linear_matrix), 1.0): + raise DeterminantDifferentFromOne(matrix=linear_matrix) + linear_matrix.flags.writeable = False + self.rotation = linear_matrix + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + s = super().__repr__() + "\n" + s += "\n".join(str(row) for row in self.rotation) + return s + + def inverse(self, name: str | None = None) -> BaseTransformationEdge: + return RotationEdge( + linear_matrix=self.rotation.T, + input=self.output, + output=self.input, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + self._validate_transform_points_shapes(points) + return points @ self.rotation.T + + def to_affine(self, name: str | None = None) -> AffineEdge: + return AffineEdge( + linear=self.rotation, + input=self.input, + output=self.output, + name=name, + ) + + +class SequenceEdge(BaseTransformationEdge): + """The Sequence transformation from the NGFF specification.""" + + def __init__( + self, + *, + name: str | None = None, + transformations: Sequence[BaseTransformationEdge], + ) -> None: + """ + Init the NgffSequence object. + + Parameters + ---------- + transformations + The transformations which compose the sequence. + Raises + ------ + EmptyTransformSequence + if `transformations` is empty + IncompatibleCoordSystemsError + if any item in `transformations` is incompatible with its neighbors + """ + if len(transformations) == 0: + raise EmptyTransformSequenceError() + previous_transf = transformations[0] + for transf_idx, current_transf in enumerate(transformations[1:], start=1): + if previous_transf.output != current_transf.input: + raise IncompatibleCoordSystemsError( + axes_error_message=( + f"Output of transformation #{transf_idx - 1} is different " + f"from input of transformation #{transf_idx}" + ), + ) + previous_transf = current_transf + self.transformations = transformations + super().__init__( + input=transformations[0].input, + output=transformations[-1].output, + name=name, + ) + + def __repr__(self) -> str: + from textwrap import indent + + out = super().__repr__() + " [\n" + for t in self.transformations: + out += indent(repr(t), prefix=" ") + "\n" + out += "]" + return out + + def inverse(self, name: str | None = None) -> SequenceEdge | None: + inverted: list[BaseTransformationEdge] = [] + for t in reversed(self.transformations): + inv = t.inverse() + if inv is None: + return None + inverted.append(inv) + return SequenceEdge(transformations=inverted, name=name) + + def to_affine(self, name: str | None = None) -> AffineEdge: + composed = self.transformations[0].to_affine().affine + for t in self.transformations[1:]: + a = t.to_affine() + composed = a.affine @ composed + return AffineEdge.from_affine_matrix( + affine_matrix=composed, + input=self.input, + output=self.output, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + return self.to_affine().transform_points(points) # FIXME + + +class ByDimensionEdge(BaseTransformationEdge): + """The ByDimension transformation from the NGFF specification.""" + + transformations: Final[Sequence[BaseTransformationEdge]] + + def __init__( + self, + *, + name: str | None = None, + transformations: Sequence[BaseTransformationEdge], + input: CoordSystem, + output: CoordSystem, + ) -> None: + """ + Parameters + ---------- + transformations + A list of transformations, whose set of output coordinate systems partition the output coordinate system of + the ByDimension transformation. + input + The input coordinate system of the transformation. + output + The output coordinate system of the transformation. + Raises + ------ + MissingAxisError + if any input axis from `transformations` is not present in `input` or + if any output axis from `transformations` is not present in `output`. + AxisRedefinitionError + if an output axis is specified by more than one item of `transformations` + UnmappedAxisError + if axis of `output` is not covered by any item of `transformations` + """ + # we check that: + # 1. each input from each transformation in self.transformation must appear in the set of input axes + # 2. each output from each transformation in self.transformation must appear at most once in the set of output + # axes + defined_output_axes: set[str] = set() + for t in transformations: + for ax in t.input.axes: + if ax not in input.axes: + raise MissingAxisError(axis_name=ax.name, cs_name=input.name) + for ax in t.output.axes: + if ax not in output.axes: + raise MissingAxisError(axis_name=ax.name, cs_name=output.name) + if ax.name in defined_output_axes: + raise AxisRedefinitionError(axis_name=ax.name) + defined_output_axes.add(ax.name) + for ax in output.axes: + if ax.name not in defined_output_axes: + raise UnmappedAxisError(axis_name=ax.name, cs_name=output.name) + + self.transformations = tuple(transformations) + super().__init__(input=input, output=output, name=name) + + def __repr__(self) -> str: + from textwrap import indent + + out = super().__repr__() + " [\n" + for t in self.transformations: + out += indent(repr(t), prefix=" ") + "\n" + out += "]" + return out + + def inverse(self, name: str | None = None) -> BaseTransformationEdge | None: + inverse_transformations: list[BaseTransformationEdge] = [] + for t in self.transformations: + inv = t.inverse() + if inv is None: + return None + inverse_transformations.append(inv) + return ByDimensionEdge( + transformations=inverse_transformations, + input=self.output, + output=self.input, + name=name, + ) + + def transform_points(self, points: ArrayLike) -> ArrayLike: + """ + Transform points (coordinates). + + Raises + ------ + UnexpectedShapeError + if `points`'s shape is incompatible with this transformation's input shape + """ + input_axes = self.input.axes_names + output_axes = self.output.axes_names + self._validate_transform_points_shapes(points) + output_columns: dict[str, ArrayLike] = {} + for t in self.transformations: + input_columns = [points[:, input_axes.index(ax)] for ax in t.input.axes_names] + input_columns_stacked: ArrayLike = np.stack(input_columns, axis=1) + output_columns_t = t.transform_points(input_columns_stacked) + for ax, col in zip(t.output.axes_names, output_columns_t.T, strict=True): + output_columns[ax] = col # type: ignore[assignment] + output: ArrayLike = np.stack([output_columns[ax] for ax in output_axes], axis=1) + return output + + def to_affine(self, name: str | None = None) -> AffineEdge: + input_axes = self.input.axes_names + output_axes = self.output.axes_names + m = np.zeros((len(output_axes) + 1, len(input_axes) + 1)) + m[-1, -1] = 1 + for t in self.transformations: + t_affine = t.to_affine() + target_output_indices = [output_axes.index(ax) for ax in t.output.axes_names if ax in output_axes] + source_output_indices = [t.output.axes_names.index(ax) for ax in t.output.axes_names] + target_input_indices = [input_axes.index(ax) for ax in t.input.axes_names] + [-1] + m[np.ix_(target_output_indices, target_input_indices)] = t_affine.affine[source_output_indices, :] + return AffineEdge.from_affine_matrix( + affine_matrix=m, + input=self.input, + output=self.output, + name=name, + ) + + +class CsGen: + """A coordinate system generator + + Use it to create coordinate systems on the fly while avoiding + repeating names + """ + + def __init__(self, base_name: str): + self._base_name = base_name + self._cs_count: int = 0 + super().__init__() + + def get_unique_name(self) -> str: + """ + Generate a unique name for a coordinate system, based on self._base_name and self._cs_count. + + Returns + ------- + A unique name string. + """ + + unique_name = f"{self._base_name}{self._cs_count}" + self._cs_count += 1 + return unique_name + + def generate(self, *, num_axes: int) -> CoordSystem: + """ + Generate a coordinate system with the specified number of axes. + + Parameters + ---------- + num_axes + The number of axes for the coordinate system. + + Returns + ------- + A CoordSystem with the specified number of axes. + """ + return CoordSystem( + name=self.get_unique_name(), + axes=tuple( + Axis( + name=f"axis_{ax_idx}", + type=NgffAxisType.SPACE, # FIXME + ) + for ax_idx in range(num_axes) + ), + virtual=True, + ) + + def generate_with_axes(self, *, axes: tuple[Axis]) -> CoordSystem: + """ + Generate a coordinate system with the specified axes. + + Parameters + ---------- + axes + The axes to use for the coordinate system. + + Returns + ------- + A CoordSystem with the specified axes. + """ + + return CoordSystem( + name=self.get_unique_name(), + axes=axes, + virtual=True, + ) + + def generate_like(self, other: CoordSystem) -> CoordSystem: + """ + Generate a coordinate system similar to another coordinate system. + + Parameters + ---------- + other + The coordinate system to use as a template. + + Returns + ------- + A CoordSystem with the same axes as the template but a unique name. + """ + return CoordSystem( + name=self.get_unique_name(), + axes=tuple( + Axis( + name=axis.name, + type=axis.type, + unit=axis.unit, + long_name=axis.long_name, + ) + for axis in other.axes + ), + virtual=True, + ) + + +def parse_identity( + model: ozm06trans.Identity, + *, + input: CoordSystem, + out: CoordSystem | CsGen, +) -> IdentityEdge: + """Parse an `Identity` NGFF transformation model into an `IdentityEdge`.""" + output = out.generate_like(input) if isinstance(out, CsGen) else out + return IdentityEdge(name=model.name, input=input, output=output) + + +def parse_translation( + model: ozm06trans.Translation, + *, + input: CoordSystem, + out: CoordSystem | CsGen, +) -> TranslationEdge: + """Parse a `Translation` NGFF transformation model into a `TranslationEdge`.""" + output = out.generate_like(input) if isinstance(out, CsGen) else out + return TranslationEdge( + translation=np.asarray(model.translation, dtype=float), + input=input, + output=output, + name=model.name, + ) + + +def parse_project_axis( + model: ozm06trans.ProjectAxis, + *, + input: CoordSystem, + output: CoordSystem | CsGen, +) -> ProjectAxisEdge: + """Parse a `ProjectAxis` NGFF transformation model into a `ProjectAxisEdge`.""" + if isinstance(output, CsGen): + num_dropped_inputs = len(model.droppedInputs or ()) + num_created_outputs = len(model.createdOutputs or ()) + num_output_axes = input.num_axes - num_dropped_inputs + num_created_outputs + output = output.generate(num_axes=num_output_axes) + + return ProjectAxisEdge( + created_outputs={output.axes[i] for i in model.createdOutputs or ()}, + dropped_inputs={input.axes[i] for i in model.droppedInputs or ()}, + input=input, + output=output, + name=model.name, + ) + + +def parse_scale( + model: ozm06trans.Scale, + *, + input: CoordSystem, + out: CoordSystem | CsGen, +) -> ScaleEdge: + """Parse a `Scale` NGFF transformation model into a `ScaleEdge`.""" + output = out.generate_like(input) if isinstance(out, CsGen) else out + return ScaleEdge( + scale=np.asarray(model.scale, dtype=float), + input=input, + output=output, + name=model.name, + ) + + +def parse_map_axis( + model: ozm06trans.MapAxis, + *, + input: CoordSystem, + out: CoordSystem | CsGen, +) -> MapAxisEdge: + """Parse a `MapAxis` NGFF transformation model into a `MapAxisEdge`.""" + if isinstance(out, CoordSystem): + output = out + else: + dummy_cs = out.generate(num_axes=len(model.mapAxis)) + output = CoordSystem( + name=dummy_cs.name, + axes=tuple(input.axes[i] for i in model.mapAxis), + virtual=True, + ) + input_to_output = { + input.axes[inp_idx]: output_axis for inp_idx, output_axis in zip(model.mapAxis, output.axes, strict=True) + } + return MapAxisEdge(input=input, output=output, name=model.name, input_to_output=input_to_output) + + +def parse_affine( + model: ozm06trans.Affine, + *, + input: CoordSystem, + output: CoordSystem | CsGen, +) -> AffineEdge: + """Parse an `Affine` NGFF transformation model into an `AffineEdge`.""" + num_output_axes = len(model.affine_matrix) # spec doesn't save last row + output = output.generate(num_axes=num_output_axes) if isinstance(output, CsGen) else output + affine_array = np.asarray(model.affine_matrix, dtype=float) + return AffineEdge( + name=model.name, linear=affine_array[:, :-1], translation=affine_array[:, -1], input=input, output=output + ) + + +def parse_rotation( + model: ozm06trans.Rotation, + *, + input: CoordSystem, + out: CoordSystem | CsGen, +) -> RotationEdge: + """Parse a `Rotation` NGFF transformation model into a `RotationEdge`.""" + num_output_axes = len(model.rotation_matrix) + output = out.generate(num_axes=num_output_axes) if isinstance(out, CsGen) else out + return RotationEdge( + name=model.name, + linear_matrix=np.asarray(model.rotation_matrix, dtype=float), + input=input, + output=output, + ) + + +def parse_sequence( + model: ozm06trans.Sequence, + *, + input: CoordSystem, + output: CoordSystem | CsGen, +) -> SequenceEdge: + """Parse a `Sequence` NGFF transformation model into a `SequenceEdge`.""" + parsed_inners: list[BaseTransformationEdge] = [] + + base_name = "intermediate" + ("" if not model.name else f"_for_{model.name}") + cs_gen: CsGen = output if isinstance(output, CsGen) else CsGen(base_name=base_name) + parsed = parse_ngff_transf( + input=input, + model=model.transformations[0], + output=cs_gen if len(model.transformations) > 1 else output, + ) + parsed_inners.append(parsed) + + for t in model.transformations[1:-1]: + parsed = parse_ngff_transf(input=parsed.output, output=cs_gen, model=t) + parsed_inners.append(parsed) + + if len(model.transformations) > 1: + parsed = parse_ngff_transf(input=parsed.output, output=output, model=model.transformations[-1]) + parsed_inners.append(parsed) + + return SequenceEdge(name=model.name, transformations=parsed_inners) + + +def parse_by_dimension( + model: ozm06trans.ByDimension, + *, + input: CoordSystem, + output: CoordSystem | CsGen, +) -> ByDimensionEdge: + """Parse a `ByDimension` NGFF transformation model into a `ByDimensionEdge`.""" + if not isinstance(output, CoordSystem): + max_out_idx = max(ax_idx for t in model.transformations for ax_idx in t.outputAxes) + output = output.generate(num_axes=max_out_idx + 1) + + piecewise_transforms: list[BaseTransformationEdge] = [] + for t in model.transformations: + inp_axes = tuple(input.axes[i] for i in t.inputAxes) + partial_input = CoordSystem( + axes=inp_axes, + name=f"{input.name}_{','.join(ax.name for ax in inp_axes)}", + virtual=True, + ) + + out_axes = tuple(output.axes[i] for i in t.outputAxes) + partial_out = CoordSystem( + axes=out_axes, + name=f"{output.name}_{','.join(ax.name for ax in inp_axes)}", + virtual=True, + ) + + parsed_t = parse_ngff_transf(model=t.transformation, input=partial_input, output=partial_out) + piecewise_transforms.append(parsed_t) + + return ByDimensionEdge( + input=input, + output=output, + name=model.name, + transformations=piecewise_transforms, + ) + + +def parse_ngff_transf( + input: CoordSystem, + model: ozm06trans.AnyTransform, + output: CoordSystem | CsGen, +) -> BaseTransformationEdge: + """Parse an NGFF coordinate transformation model into a `BaseTransformationEdge`""" + if isinstance(model, ozm06trans.Identity): + return parse_identity(model, input=input, out=output) + elif isinstance(model, ozm06trans.Translation): + return parse_translation(model, input=input, out=output) + elif isinstance(model, ozm06trans.Scale): + return parse_scale(model, input=input, out=output) + elif isinstance(model, ozm06trans.MapAxis): + return parse_map_axis(model, input=input, out=output) + elif isinstance(model, ozm06trans.Affine): + return parse_affine(model, input=input, output=output) + elif isinstance(model, ozm06trans.Rotation): + return parse_rotation(model, input=input, out=output) + elif isinstance(model, ozm06trans.Sequence): + return parse_sequence(model, input=input, output=output) + elif isinstance(model, ozm06trans.ByDimension): + return parse_by_dimension(model, input=input, output=output) + elif isinstance(model, ozm06trans.ProjectAxis): + return parse_project_axis(model, input=input, output=output) + else: + raise NotImplementedError(f"Unsupported transformation: {model.type}") diff --git a/src/spatialdata/transformations/_graph/exceptions.py b/src/spatialdata/transformations/_graph/exceptions.py new file mode 100644 index 000000000..8761a8f18 --- /dev/null +++ b/src/spatialdata/transformations/_graph/exceptions.py @@ -0,0 +1,277 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, override + +if TYPE_CHECKING: + pass + + +class CoordinateSystemNotFoundError(ValueError): + """ + Exception raised when a coordinate system is not found in the transformation manager. + + Attributes + ---------- + name : str + The name of the coordinate system that was not found. + """ + + def __init__(self, name: str) -> None: + self.name: str = name + super().__init__(f"Coordinate system '{name}' not found in the transformation manager.") + + +class CoordinateSystemAlreadyExistsError(ValueError): + """ + Exception raised when a coordinate system already exists in the transformation manager. + + Attributes + ---------- + name : str + The name of the coordinate system to check. + """ + + def __init__(self, name: str) -> None: + self.name: str = name + super().__init__(f"Coordinate system '{name}' already exists in the transformation manager.") + + +class ElementNotRegisteredToAnyCoordinateSystemError(KeyError): + """ + Exception raised when an element is not found in the transformation manager. + + Attributes + ---------- + element_name : str + The name of the element that was not found. + """ + + def __init__(self, element_name: str) -> None: + self.element_name: str = element_name + super().__init__(f"Element '{element_name}' has not been registered to any coordinate system.") + + +class TransformationNotFoundError(KeyError): + """ + Exception raised when a transformation is not found between coordinate systems. + + Attributes + ---------- + input_cs_name : str + The name of the input coordinate system. + output_cs_name : str + The name of the output coordinate system. + edge_key: str or None + key used when adding transformation + """ + + def __init__(self, source_cs_name: str, target_cs_name: str, edge_key: str | None = None) -> None: + self.input_cs_name: str = source_cs_name + self.output_cs_name: str = target_cs_name + self.edge_key: str | None = edge_key + msg = f"Transformation from '{source_cs_name}' to '{target_cs_name}' not found" + if edge_key is not None: + msg += f" with key '{edge_key}'" + super().__init__(msg) + + +class ElementAlreadyExistsError(ValueError): + """ + Exception raised when trying to add an Element that already exists. + + Attributes + ---------- + name : str + The name of the Element that already exists. + """ + + def __init__(self, name: str) -> None: + self.name: str = name + super().__init__(f"Element '{name}' already exists in the transformation manager") + + +class InvalidPathError(ValueError): + """ + Exception raised when a path is defined with less than 2 nodes. + + Attributes + ---------- + invalid_path : + list of coordinate system names + """ + + def __init__(self, invalid_path: list[str]) -> None: + super().__init__(f"Found an invalid path with less than 2 nodes: {invalid_path}") + + +class TransformationPathNotFoundError(ValueError): + """ + Exception raised when no transformation path exists between coordinate systems. + + Attributes + ---------- + source_cs_name : str + The name of the source coordinate system. + target_cs_name : str + The name of the target coordinate system. + """ + + def __init__(self, source_cs_name: str, target_cs_name: str) -> None: + self.source_cs_name: str = source_cs_name + self.target_cs_name: str = target_cs_name + super().__init__(f"No transformation path found from {source_cs_name} to {target_cs_name}") + + +class TransformationPathAmbiguousError(ValueError): + """ + Exception raised when multiple transformation path exists between coordinate systems. + + Attributes + ---------- + source_cs_name : str + The name of the source coordinate system. + target_cs_name : str + The name of the target coordinate system. + """ + + def __init__(self, source_cs_name: str, target_cs_name: str) -> None: + self.source_cs_name: str = source_cs_name + self.target_cs_name: str = target_cs_name + base_msg = f"Transformation Path ambiguous from {source_cs_name} to {target_cs_name}." + cause_of_confusion = self.cause_of_confusion() + msg = f"{base_msg} {cause_of_confusion}" if cause_of_confusion else base_msg + super().__init__(msg) + + def cause_of_confusion(self) -> str: + return "" + + +class TransformationPathAmbiguousNoEdgeExpectedError(TransformationPathAmbiguousError): + def __init__(self, source_cs_name: str, target_cs_name: str) -> None: + super().__init__(source_cs_name, target_cs_name) + + @override + def cause_of_confusion(self) -> str: + + return "Multiple edges found. None of them were specified to be expected" + + +class TransformationPathAmbiguousMultipleEdgeExpectedError(TransformationPathAmbiguousError): + def __init__(self, source_cs_name: str, target_cs_name: str, number_of_edges_expected: int) -> None: + self.number_of_edges_expected: int = number_of_edges_expected + super().__init__(source_cs_name, target_cs_name) + + @override + def cause_of_confusion(self) -> str: + return f"Multiple ({self.number_of_edges_expected}) edges were specified to be expected" + + +class TransformationPathNotSimple(ValueError): + """ + Exception raised when a path represented by a list of coordinate systems is not simple. + + A simple path is one in which each coordinate system appears only once + + Attributes + ---------- + path: + list of coordinate system names + """ + + def __init__(self, path: list[str]) -> None: + self.path = path + super().__init__( + "Transformation Path not simple, i.e., some coordinate systems appear multiple times:\n" + f"{path}" + ) + + +class CannotRemoveCoordinateSystemError(ValueError): + """Exception raised when trying to remove a coordinate system. + + Attributes + ---------- + name: str + name of the coordinate system + """ + + def __init__(self, name: str) -> None: + self.name: str = name + super().__init__(f"Cannot remove coordinate system with name {name}.") + + +class CoordinateSystemHasTransformationsError(ValueError): + """ + Exception raised when trying to remove a coordinate system that has associated transformations. + + Attributes + ---------- + name : str + The name of the coordinate system that has transformations. + """ + + def __init__(self, name: str) -> None: + self.name: str = name + super().__init__(f"Coordinate System ('{name}') has transformations.") + + +class CoordinateSystemHasElementsError(ValueError): + """ + Exception raised when trying to remove a coordinate system that has associated elements. + + Attributes + ---------- + name : str + The name of the coordinate system that has elements. + associated_elements : list[str] + List of element names associated with the coordinate system. + """ + + def __init__(self, name: str, associated_elements: list[str]) -> None: + self.name: str = name + self.associated_elements: list[str] = associated_elements + super().__init__(f"Coordinate system '{name}' has elements belonging to it: {associated_elements}") + + +class TransformationManagerWarning(UserWarning): + """Base warning category for TransformationManager.""" + + pass + + +class AxisParsingException(Exception): + pass + + +class NGFFCompatibilityError(Exception): + """ + Exception raised when a transformation cannot be mapped to an NGFF transformation edge. + + Attributes + ---------- + msg : str + The error message describing the compatibility issue. + """ + + def __init__(self, msg: str): + self.msg = msg + super().__init__(msg) + + +class IncompatibleTransformationOutputCoordinateSystemError(Exception): + """ + Exception raised the expected output coordinate system is incompatible with an input coordinate system and a transformation edge. + """ + + def __init__( + self, input_coordinate_system_name: str, transformation_class_name: str, output_coordinate_system_name: str + ) -> None: + + self.input_coordinate_system_name = input_coordinate_system_name + self.output_coordinate_system_name = output_coordinate_system_name + self.transformation_class_name = transformation_class_name + super().__init__( + f"Transformation incompatible with input/output coordinate systems: " + f"\nInput coordinate system name: {input_coordinate_system_name}" + f"\nOutput coordinate system name: {output_coordinate_system_name}" + f"\nTransformation type: {transformation_class_name}" + ) diff --git a/src/spatialdata/transformations/_graph/graph.py b/src/spatialdata/transformations/_graph/graph.py new file mode 100644 index 000000000..0ab3421c4 --- /dev/null +++ b/src/spatialdata/transformations/_graph/graph.py @@ -0,0 +1,742 @@ +from __future__ import annotations + +import hashlib +import uuid +from collections.abc import Sequence +from typing import Any, override + +import networkx as nx + +from spatialdata.literals import TRANSFORM_KEY +from spatialdata.transformations._graph.edges.ngff_aligned_edges import BaseTransformationEdge, SequenceEdge +from spatialdata.transformations._graph.exceptions import ( + CannotRemoveCoordinateSystemError, + CoordinateSystemAlreadyExistsError, + CoordinateSystemHasElementsError, + CoordinateSystemHasTransformationsError, + CoordinateSystemNotFoundError, + ElementAlreadyExistsError, + ElementNotRegisteredToAnyCoordinateSystemError, + InvalidPathError, + TransformationNotFoundError, + TransformationPathAmbiguousError, + TransformationPathAmbiguousMultipleEdgeExpectedError, + TransformationPathAmbiguousNoEdgeExpectedError, + TransformationPathNotFoundError, + TransformationPathNotSimple, +) +from spatialdata.transformations._graph.vert import CoordSystem + + +class TransformationGraph: + def __init__(self) -> None: + """ + Initialize a TransformationManager with empty graph and mappings. + + The TransformationManager maintains a graph of coordinate systems and transformations between them, + as well as a mapping of element names to their respective coordinate systems. + """ + self._graph: nx.MultiDiGraph[CoordSystem, dict[str, Any], dict[str, BaseTransformationEdge]] = nx.MultiDiGraph() + # MultiDiGraph with CoordSystem objects as nodes and BaseTransformationEdge as edge attributes + self._element_to_cs_mapping: dict[str, CoordSystem] = {} + # mapping element_name to the coordinate system to which the element belongs + + def assert_element_exists(self, element_name: str) -> None: + """ + Assert that an element exists in the transformation manager. + + Parameters + ---------- + element_name + The name of the element to check. + + Raises + ------ + ElementNotRegisteredToAnyCoordinateSystemError + If the element does not exist. + """ + if element_name not in self._element_to_cs_mapping: + raise ElementNotRegisteredToAnyCoordinateSystemError(element_name) + + def assert_element_does_not_exist(self, element_name: str) -> None: + """ + Assert that an element doesn't exist in the transformation manager. + + Parameters + ---------- + element_name + The name of the element to check. + + Raises + ------ + ElementAlreadyExistsError + If the element already exists. + """ + if element_name in self._element_to_cs_mapping: + raise ElementAlreadyExistsError(element_name) + + def assert_coordinate_system_exists(self, cs: CoordSystem) -> None: + """ + Assert that coordinate system exists in the graph. + + Parameters + ---------- + cs + The coordinate system to check. + + Raises + ------ + CoordinateSystemNotFoundError + If the coordinate system does not exist. + """ + if cs not in self._graph: + raise CoordinateSystemNotFoundError(cs.name) + + def assert_coordinate_system_does_not_exist(self, cs: CoordSystem) -> None: + """ + Assert that coordinate system doesn't exist in the graph. + + Parameters + ---------- + cs + The coordinate system to check. + + Raises + ------ + CoordinateSystemAlreadyExistsError + If the coordinate system already exists. + """ + if cs in self._graph: + raise CoordinateSystemAlreadyExistsError(cs.name) + + def assert_an_edge_exists_between_coordinate_systems( + self, source_cs: CoordSystem, target_cs: CoordSystem, edge_key: str | None = None + ) -> None: + """ + Assert that an edge exists between coordinate systems. + + Parameters + ---------- + source_cs + The input coordinate system. + target_cs + The output coordinate system. + edge_key + If specified, asserts that the specific edge exists between coordinate systems. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + TransformationNotFoundError + If the edge does not exist. + """ + self.assert_coordinate_system_exists(source_cs) + self.assert_coordinate_system_exists(target_cs) + if not self._graph.has_edge(source_cs, target_cs): + raise TransformationNotFoundError(source_cs.name, target_cs.name, edge_key=edge_key) + + def assert_coordinate_system_has_no_transformations(self, cs: CoordSystem) -> None: + """ + Assert that coordinate system has no associated transformations. + + Parameters + ---------- + cs + The coordinate system to check. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + CoordinateSystemHasTransformationsError + If the coordinate system has associated transformations. + + """ + self.assert_coordinate_system_exists(cs) + has_successors = next(self._graph.successors(cs), None) is not None + has_predecessors = next(self._graph.predecessors(cs), None) is not None + + if has_successors or has_predecessors: + raise CoordinateSystemHasTransformationsError(cs.name) + + def assert_coordinate_system_has_no_elements(self, cs: CoordSystem) -> None: + """ + Assert that a coordinate system doesn't have elements belonging to it. + + Parameters + ---------- + cs + The coordinate system to check. + + Raises + ------ + CoordinateSystemHasElementsError + If the coordinate system has elements belonging to it. + """ + elements = self._get_elements_belonging_to_cs(cs) + # also checks if cs exists + if elements: + raise CoordinateSystemHasElementsError(cs.name, elements) + + def add_coordinate_system(self, cs: CoordSystem) -> None: + """ + Add a coordinate system to the transformation manager. + + Parameters + ---------- + cs + The coordinate system to add. + + Raises + ------ + CoordinateSystemAlreadyExistsError + If the coordinate system already exists. + """ + self.assert_coordinate_system_does_not_exist(cs) + self._graph.add_node(cs) + + def get_coordinate_system(self, cs_name: str) -> CoordSystem | None: + """ + Search for a coordinate system by name, return it if found, else return None + + Parameters + ---------- + cs_name + Name of a coordinate system + """ + + for cs in self.list_coordinate_systems(): + if cs.name == cs_name: + return cs + + return None + + def get_elements_in_coordinate_system(self, cs_name: str) -> list[str]: + """ + Returns a list of element names that belong to this coordinate system + + Parameters + ---------- + cs_name + Name of a coordinate system + + Returns + ---- + list of element names + """ + + return [element_name for element_name, cs in self._element_to_cs_mapping.items() if cs.name == cs_name] + + def remove_coordinate_system(self, cs: CoordSystem) -> None: + """ + Remove a coordinate system from the transformation manager. + + Parameters + ---------- + cs + The coordinate system to remove. + + Raises + ------ + CoordinateSystemNotFoundError + If the coordinate system is not found. + CannotRemoveCoordinateSystemError + If the coordinate system cannot be removed. + Caused by CoordinateSystemHasTransformationsError or + CoordinateSystemHasElementsError when the system has dependencies. + """ + try: + self.assert_coordinate_system_has_no_transformations(cs) + # also asserts that cs exists + self.assert_coordinate_system_has_no_elements(cs) + except (CoordinateSystemHasTransformationsError, CoordinateSystemHasElementsError) as err: + raise CannotRemoveCoordinateSystemError(cs.name) from err + + self._graph.remove_node(cs) + + def list_coordinate_systems(self) -> list[CoordSystem]: + """ + List all registered coordinate systems. + + Returns + ------- + A list of coordinate system objects. + """ + return list(self._graph.nodes()) + + def add_element(self, element_name: str, coordinate_system: CoordSystem) -> None: + """ + Register an element and associate it with a coordinate system. + + Parameters + ---------- + element_name + The name of the element. + coordinate_system + The coordinate system to which the element belongs. + + Raises + ------ + CoordinateSystemNotFoundError + If the coordinate system is not found. + """ + self.assert_element_does_not_exist(element_name) + self.assert_coordinate_system_exists(coordinate_system) + self._element_to_cs_mapping[element_name] = coordinate_system + + def get_element_coordinate_system(self, element_name: str) -> CoordSystem: + """ + Get the coordinate system to which an element belongs. + + Parameters + ---------- + element_name + The name of the element. + + Returns + ------- + The coordinate system. + + Raises + ------ + ElementNotFoundError + If the element does not exist. + """ + self.assert_element_exists(element_name) + return self._element_to_cs_mapping[element_name] + + def get_outgoing_edges( + self, element_name: str, target_coordinate_system_names: Sequence[str] | None = None + ) -> dict[str, BaseTransformationEdge]: + """ + Get outgoing transformation edges from an element's coordinate system. + + Parameters + ---------- + element_name + The name of the element. + target_coordinate_system_names + Optional list of target coordinate system names to filter by. + + Returns + ------- + dict, with edge keys as keys and outgoing transformation edges from the element's coordinate system as values. + """ + element_coordinate_system = self.get_element_coordinate_system(element_name) + outgoing_edges: dict[str, BaseTransformationEdge] = {} + + successors_to_consider = list(self._graph.successors(element_coordinate_system)) + if target_coordinate_system_names is not None: + successors_to_consider = [x for x in successors_to_consider if x.name in target_coordinate_system_names] + + for successor_coordinate_system in successors_to_consider: + outgoing_edge_this_successor = self.get_direct_transformations( + element_coordinate_system, successor_coordinate_system + ) + outgoing_edges |= outgoing_edge_this_successor + + return outgoing_edges + + def get_incoming_edges( + self, element_name: str, source_coordinate_system_names: Sequence[str] | None = None + ) -> dict[str, BaseTransformationEdge]: + """ + Get incoming transformation edges to an element's coordinate system. + + Parameters + ---------- + element_name + The name of the element. + source_coordinate_system_names + Optional list of source coordinate system names to filter by. + + Returns + ------- + dict, with edge keys as keys and transformations from the element's coordinate system as values. + """ + element_coordinate_system = self.get_element_coordinate_system(element_name) + incoming_edges: dict[str, BaseTransformationEdge] = {} + + predecessors_to_consider = list(self._graph.successors(element_coordinate_system)) + if source_coordinate_system_names is not None: + predecessors_to_consider = [x for x in predecessors_to_consider if x.name in source_coordinate_system_names] + + for predecessor_coordinate_system in predecessors_to_consider: + incoming_edge_this_successor = self.get_direct_transformations( + predecessor_coordinate_system, element_coordinate_system + ) + incoming_edges |= incoming_edge_this_successor + + return incoming_edges + + def _unset_element(self, element_name: str) -> None: + """ + Unregister an element from the coordinate system to which it belongs. + + Parameters + ---------- + element_name + The name of the element. + + Raises + ------ + ElementNotFoundError + If the element has not been registered to any coordinate system. + """ + self.assert_element_exists(element_name) + del self._element_to_cs_mapping[element_name] + + @staticmethod + def _calc_edge_key(edge: BaseTransformationEdge) -> str: + """ + Encode attributes of an edge into a string to be used as graph-wide edge identifier. + + The string consists of a deterministic part -- {input cs name}->{output cs name}:{edge class name}, followed + by a random part. + + It is better to not make edge key a deterministic function of the transformation edge (e.g., using id()), + because, when the transformation graph gets deep copied, the transformation edge objects will be recreated, + but the edge keys will just be copied, leading to inconsistencies. + + Parameters + ---------- + edge + The transformation edge to encode. + + Returns + ------- + edge_key + str, the encoded edge identifier. + """ + base = f"{edge.input.name}->{edge.output.name}:{type(edge).__name__}" + salt = uuid.uuid4().hex + suffix = hashlib.sha256((base + salt).encode()).hexdigest()[:12] # short hash + return f"{base}-{suffix}" + + def add_transformation_edge( + self, + transformation_edge: BaseTransformationEdge, + ) -> str: + """ + Add an edge representing a transformation between coordinate systems. + + Parameters + ---------- + transformation_edge + The transformation edge to add. + + Returns + ------- + edge_key + str, the encoded edge identifier. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + """ + source_cs, target_cs = transformation_edge.input, transformation_edge.output + self.assert_coordinate_system_exists(source_cs) + self.assert_coordinate_system_exists(target_cs) + + edge_key = self._calc_edge_key(transformation_edge) + edge_attributes = {TRANSFORM_KEY: transformation_edge} + self._graph.add_edge(source_cs, target_cs, key=edge_key, **edge_attributes) + + return edge_key + + def get_direct_transformations( + self, source_cs: CoordSystem, target_cs: CoordSystem + ) -> dict[str, BaseTransformationEdge]: + """ + Retrieve transformations directly defined between coordinate systems. + + Parameters + ---------- + source_cs + The input coordinate system. + target_cs + The output coordinate system. + + Returns + ------- + dict, with edge keys as keys and transformations as values (empty if there are no direct transformations). + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + """ + self.assert_coordinate_system_exists(source_cs) + self.assert_coordinate_system_exists(target_cs) + + transforms = {} + if self._graph.has_edge(source_cs, target_cs): + for edge_key, edge_data in self._graph[source_cs][target_cs].items(): + transform: BaseTransformationEdge = edge_data[TRANSFORM_KEY] + transforms[edge_key] = transform + return transforms + + def remove_specific_transformation( + self, + edge_key: str, + transformation_edge: BaseTransformationEdge, + ) -> None: + """ + Remove a specific transformation between coordinate systems. + + Parameters + ---------- + edge_key + unique id of the edge + transformation_edge + The transformation edge to remove. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + TransformationNotFoundError + If the transformation does not exist. + """ + source_cs, target_cs = transformation_edge.input, transformation_edge.output + self.assert_an_edge_exists_between_coordinate_systems(source_cs, target_cs, edge_key=edge_key) + # also checks if source_cs and target_cs exist + self._graph.remove_edge(source_cs, target_cs, key=edge_key) + + def remove_all_transformations_between_coordinate_systems( + self, + source_cs: CoordSystem, + target_cs: CoordSystem, + ) -> None: + """ + Remove all transformation between coordinate systems. + + Parameters + ---------- + source_cs + The input coordinate system. + target_cs + The output coordinate system. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + TransformationNotFoundError + If no transformation exists between the coordinate systems. + """ + self.assert_coordinate_system_exists(source_cs) + self.assert_coordinate_system_exists(target_cs) + if self._graph.has_edge(source_cs, target_cs): + transform_edges_to_remove = {} + for edge_key in self._graph[source_cs][target_cs]: + transform_edges_to_remove[edge_key] = self._graph[source_cs][target_cs][edge_key][TRANSFORM_KEY] + + for edge_key, transform_edge in transform_edges_to_remove.items(): + self._graph.remove_edge(transform_edge.input, transform_edge.output, key=edge_key) + + def _get_transformation_sequences_from_simple_paths_after_disambiguation( + self, + paths: list[list[CoordSystem]], + expected_intermediate_edge_keys: Sequence[str] | None, + ) -> list[SequenceEdge]: + """ + Traverses paths to form sequence of Transformations. + + In case of ambiguity looks into `expected_intermediate_transformations` to disambiguate. + + Parameters + ---------- + paths + Sequence of list of nodes. + expected_intermediate_edge_keys + List of edges, for use when multiple edges are found between two coordinate systems in a path. + + Returns + ------- + List of transformation edges, each an instance of the transformation class Sequence. + + Raises + ------ + TransformationPathNotSimple + If any path in `paths` is not simple, i.e., has recurring coordinate systems. + """ + for path in paths: + if len(set(path)) != len(path): + raise TransformationPathNotSimple([x.name for x in path]) + + expected_intermediate_transformation_edge_keys = set[str]() + if expected_intermediate_edge_keys is not None: + for edge_key in expected_intermediate_edge_keys: + expected_intermediate_transformation_edge_keys |= {edge_key} + + all_sequences = [] + deduplicated_paths = list({repr(x): x for x in paths}.values()) + for path in deduplicated_paths: + if len(path) <= 1: + raise InvalidPathError([x.name for x in path]) + + transformation_list: list[BaseTransformationEdge] = [] + for i in range(len(path) - 1): + source_cs_here, target_cs_here = path[i], path[i + 1] + edge_data = self._graph[source_cs_here][target_cs_here] + if len(edge_data) > 1: + # when there are multiple edges between a pair of coordinate systems in the path + # find the expected edge based on key + expected_intermediate_transformation_key_here = ( + expected_intermediate_transformation_edge_keys & set(edge_data.keys()) + ) + if len(expected_intermediate_transformation_key_here) == 0: + # transformation was not specified in `expected_intermediate_transformations` for disambiguation + raise TransformationPathAmbiguousNoEdgeExpectedError(source_cs_here.name, target_cs_here.name) + if len(expected_intermediate_transformation_key_here) > 1: + # multiple transformations were specified in `expected_intermediate_transformations` + # for disambiguation + raise TransformationPathAmbiguousMultipleEdgeExpectedError( + source_cs_here.name, target_cs_here.name, len(expected_intermediate_transformation_key_here) + ) + edge_key_to_use = list(expected_intermediate_transformation_key_here)[0] + transformation_list.append(edge_data[edge_key_to_use][TRANSFORM_KEY]) + else: + # Only one edge, no ambiguity + edge_key = next(iter(edge_data.keys())) + transformation_list.append(edge_data[edge_key][TRANSFORM_KEY]) + all_sequences.append(SequenceEdge(transformations=transformation_list)) + return all_sequences + + def get_all_shortest_transformation_sequences( + self, + source_cs: CoordSystem, + target_cs: CoordSystem, + expected_intermediate_edge_keys: Sequence[str] = (), + ) -> list[SequenceEdge]: + """ + Get all shortest sequences of transformations between two coordinate systems. + + Parameters + ---------- + source_cs + The source coordinate system. + target_cs + The target coordinate system. + expected_intermediate_edge_keys + List of keys of intermediate edges. + Used to choose an edge when multiple edges are found between the same coordinate systems. + + Returns + ------- + list[Sequence] + All shortest sequences of transformations from source_cs to target_cs. + When multiple transformations are defined between the same coordinate systems, only those containing + a transformation among intermediate transformations are included. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + TransformationPathNotFoundError + If no path exists between the source and target coordinate systems. + TransformationPathAmbiguousError + When multiple transformations are defined between the same coordinate systems and transformations + are not specified in `expeceted_intermediate_transformations` for disambiguation. + """ + try: + paths = list(nx.all_shortest_paths(self._graph, source=source_cs, target=target_cs)) + except nx.NetworkXNoPath as nxe: + raise TransformationPathNotFoundError(source_cs.name, target_cs.name) from nxe + + try: + return self._get_transformation_sequences_from_simple_paths_after_disambiguation( + paths, expected_intermediate_edge_keys + ) + except TransformationPathAmbiguousError as tpae: + raise TransformationPathAmbiguousError(source_cs.name, target_cs.name) from tpae + + def get_all_transformation_sequences( + self, + source_cs: CoordSystem, + target_cs: CoordSystem, + expected_intermediate_edge_keys: Sequence[str] = (), + ) -> list[SequenceEdge]: + """ + Get all existing sequences of transformations between two coordinate systems. + + Parameters + ---------- + source_cs + The source coordinate system. + target_cs + The target coordinate system. + expected_intermediate_edge_keys + List of keys of intermediate edges expected. + Used to choose an edge when multiple edges are found between the same coordinate systems. + + Returns + ------- + list[Sequence] + All existing sequences of transformations from source_cs to target_cs. + When multiple transformations are defined between the same coordinate systems, only those containing + a transformation among intermediate transformations are included. + + Raises + ------ + CoordinateSystemNotFoundError + If either coordinate system does not exist. + TransformationPathNotFoundError + If no path exists between the source and target coordinate systems. + TransformationPathAmbiguousError + When multiple transformations are defined between the same coordinate systems and transformations + are not specified in `expected_intermediate_transformations` for disambiguation. + """ + try: + paths = list(nx.all_simple_paths(self._graph, source=source_cs, target=target_cs)) + except nx.NetworkXNoPath as nxe: + raise TransformationPathNotFoundError(source_cs.name, target_cs.name) from nxe + + try: + return self._get_transformation_sequences_from_simple_paths_after_disambiguation( + paths, expected_intermediate_edge_keys + ) + except TransformationPathAmbiguousError as tpae: + raise tpae from TransformationPathAmbiguousError(source_cs.name, target_cs.name) + + def _get_elements_belonging_to_cs(self, cs: CoordSystem) -> list[str]: + """ + Get all elements belonging to a coordinate system. + + Parameters + ---------- + cs + The coordinate system to check. + + Returns + ------- + List of element names belonging to the coordinate system. + + Raises + ------ + CoordinateSystemNotFoundError + If the coordinate system does not exist. + + """ + self.assert_coordinate_system_exists(cs) + + elements = [] + for element_name, element_cs in self._element_to_cs_mapping.items(): + if element_cs == cs: + elements.append(element_name) + + return elements + + @override + def __repr__(self) -> str: + """Return a string representation of the TransformationManager.""" + return ( + f"TransformationManager(" + f" coordinate_systems={list(self._graph.nodes())}, " + f" coordinate_transforms={[x[TRANSFORM_KEY] for *_, x in self._graph.edges(data=True)]}, " + f" elements={list(self._element_to_cs_mapping.keys())}" + f")" + ) diff --git a/src/spatialdata/transformations/_graph/types.py b/src/spatialdata/transformations/_graph/types.py new file mode 100644 index 000000000..03f20636e --- /dev/null +++ b/src/spatialdata/transformations/_graph/types.py @@ -0,0 +1,3 @@ +from __future__ import annotations + +ValidAxis_t = str diff --git a/src/spatialdata/transformations/_graph/utils.py b/src/spatialdata/transformations/_graph/utils.py new file mode 100644 index 000000000..c997313fd --- /dev/null +++ b/src/spatialdata/transformations/_graph/utils.py @@ -0,0 +1,150 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +import pandas as pd + +from spatialdata.literals import X, Y, Z +from spatialdata.models import get_axes_names +from spatialdata.transformations._graph.vert import Axis, CoordSystem + +if TYPE_CHECKING: + from spatialdata.models import SpatialElement + from spatialdata.transformations._graph.types import ValidAxis_t + + +def get_intrinsic_coordinate_system_name(element_key: str) -> str: + + return f"{element_key}_intrinsic" + + +def get_instrinsic_coordinate_system( + element_to_be_added_to_graph: SpatialElement, coordinate_system_name: str +) -> CoordSystem: + """ + Create a coordinate system object from a spatial element, which can be added to a graph. + + Parameters + ---------- + element_to_be_added_to_graph + The spatial element for which to create a coordinate system. + coordinate_system_name + The name to use for the coordinate system. + + Returns + ------- + A CoordSystem object representing the default coordinate system for the element. + """ + ngff_axes = [] + for axis_name in get_axes_names(element_to_be_added_to_graph): + ngff_axis = Axis.from_spatialdata_axis_name(axis_name) + ngff_axes.append(ngff_axis) + return CoordSystem(name=coordinate_system_name, axes=tuple(ngff_axes), virtual=False) + + +def validate_axis_name(axis: ValidAxis_t) -> None: + """ + Check if the axis name is valid. + + Parameters + ---------- + axis + The axis name + + Raises + ------ + TypeError + If the axis name not in ["c", "x", "y", "z"] + """ + if axis not in ["c", "x", "y", "z"]: + raise TypeError(f"Invalid axis: {axis}") + + +def validate_axes(axes: tuple[ValidAxis_t, ...]) -> None: + """ + Check if the names of the axes are valid. + + Parameters + ---------- + axis + The names of the axes + + Raises + ------ + TypeError + If not in ["c", "x", "y", "z"] + """ + for ax in axes: + validate_axis_name(ax) + if len(axes) != len(set(axes)): + raise ValueError("Axes must be unique.") + + +def get_spatial_axes(axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + """ + Get the spatial axes of interest. + + Parameters + ---------- + axes + Should be a subset of ['x', 'y', 'z', 'c'] + + Returns + ------- + The spatial axes, i.e. the input axes but without 'c' + """ + validate_axes(axes) + return tuple(ax for ax in axes if ax in [X, Y, Z]) + + +def get_union_of_axes(axes: tuple[Axis, ...]) -> tuple[Axis, ...]: + """ + Produce a union of ``Axis`` objects, dropping duplicates and ordering reverse alphabetically. + + The ordering favors non‑None values for ``unit`` and ``long_name``. + Axes that share the same ``name`` **and** ``type`` are considered duplicates; + only the first occurrence (according to the sorted order) is retained. + """ + # ------------------------------------------------------------------ + # Build a DataFrame that holds the original Axis objects together with + # their attributes and auxiliary “non‑null” flags. + # ------------------------------------------------------------------ + df = pd.DataFrame( + { + "axis": pd.Series(axes, dtype="object"), # original Axis objects + "name": pd.Series([ax.name for ax in axes], dtype="string"), + "type": pd.Series([ax.type for ax in axes], dtype="object"), + "unit": pd.Series([ax.unit for ax in axes], dtype="string"), + "long_name": pd.Series([ax.long_name for ax in axes], dtype="string"), + } + ) + + # Remove duplicates based on the logical key (name, type, unit, long_name). + df = df.drop_duplicates(subset=["name", "type", "unit", "long_name"], keep="first") + + # Boolean columns that indicate the presence of a value. + df["unit_not_none"] = df["unit"].notna() + df["long_name_not_none"] = df["long_name"].notna() + + # Sort axes in descending (alphabetical) order + df = df.sort_values( + by=[ + "type", + "name", + "unit_not_none", + "unit", + "long_name_not_none", + "long_name", + ], + ascending=[ + True, # NgffAxisType has CHANNEL first, SPACE next, so that we get channel axes before space axes + False, + False, # so that those with unit are placed before those without + False, + False, # so that those with long_name are placed before those without + False, + ], + kind="mergesort", # stable sort to preserve order of equal keys + ) + + return tuple(df["axis"]) diff --git a/src/spatialdata/transformations/_graph/vert.py b/src/spatialdata/transformations/_graph/vert.py new file mode 100644 index 000000000..c7688d720 --- /dev/null +++ b/src/spatialdata/transformations/_graph/vert.py @@ -0,0 +1,184 @@ +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass +from typing import Final + +import ome_zarr_models.v06.coordinate_transforms as ozm06ct + +from spatialdata.literals import NgffAxisType, axis_type_mapping_ngff +from spatialdata.transformations._graph.exceptions import AxisParsingException + + +@dataclass(frozen=True) +class Axis: + """Representation of a coordinate system axis""" + + name: Final[str] + type: Final[NgffAxisType] + unit: Final[str | None] = None + "unit of the axis. For a set of valid options see https://ngff.openmicroscopy.org/" + long_name: Final[str | None] = None + "a longer, human-friendly name for this axis" + + def cloned_with(self, *, unit: str | None) -> Axis: + return Axis(name=self.name, type=self.type, unit=unit or self.unit, long_name=self.long_name) + + def __hash__(self) -> int: + return hash((self.name, self.type, self.unit, self.long_name)) + + def __repr__(self) -> str: + return f"Axis(name={self.name}, type={self.type})" + + def __eq__(self, value: object, /) -> bool: + if not isinstance(value, Axis): + return False + return ( + self.name == value.name + and self.type == value.type + and self.unit == value.unit + and self.long_name == value.long_name + ) + + @classmethod + def try_from_model(cls, model: ozm06ct.Axis) -> Axis: + """ + Parse an `Axis` from an ome-zarr-models axis model. + + Raises + ------ + AxisParsingException + if `model` doesn't have a name, has a type other than "channel" or "space", or has a unit + that isn't a string or None + """ + + name = model.name + if name is None: + raise AxisParsingException("Axis doesn't have a name") + if model.type is None: + raise AxisParsingException("Axis doesn't have a type") + try: + axis_type = NgffAxisType(model.type) + except ValueError as ve: + raise AxisParsingException(f"Can't handle axis of type {model.type}") from ve + if not isinstance(model.unit, (str, type(None))): + raise AxisParsingException("Can't handle axis unit") + return Axis( + name=name, + type=axis_type, + unit=model.unit, + long_name=model.longName, + ) + + @classmethod + def from_spatialdata_axis_name(cls, axis_name: str, unit: None = None) -> Axis: + """ + Create an Axis from a spatialdata axis name. + + Parameters + ---------- + axis_name + The name of the axis (e.g., 'x', 'y', 'z', 'c'). + unit + Optional unit for the axis. + + Returns + ------- + An Axis object with the specified name and corresponding type. + """ + return cls(name=axis_name, type=axis_type_mapping_ngff[axis_name], unit=unit) + + +class CoordSystemParsingException(Exception): + pass + + +class DuplicateAxisNameError(Exception): + def __init__(self, *, axis_name: str) -> None: + self.axis_name = axis_name + super().__init__(f"Axis name '{axis_name}' is used more than once") + + +@dataclass(frozen=True) +class CoordSystem: + """ + Representation of a coordinate system. + """ + + name: Final[str] + axes: Final[tuple[Axis, ...]] + virtual: Final[bool] + """A virtual coordinate system exists as an intermediate step between + non-virtual coordinate systems and is usually ignored during serialization""" + + def __post_init__(self) -> None: + """ + Raises + ------ + DuplicateAxisNameError + if `axes` contains axes with duplicate names + """ + seen_names: set[str] = set() + for axis in self.axes: + if axis.name in seen_names: + raise DuplicateAxisNameError(axis_name=axis.name) + seen_names.add(axis.name) + + def __repr__(self) -> str: + return f"{self.__class__.__name__}({self.name!r}, {self.axes})" + + def __hash__(self) -> int: + return hash((self.name, self.axes, self.virtual)) + + @classmethod + def try_from_model(cls, model: ozm06ct.CoordinateSystem) -> CoordSystem: + """ + Parse a `CoordSystem` from an ome-zarr-models coordinate system model. + + Raises + ------ + AxisParsingException + if any axis in `model.axes` fails to parse + DuplicateAxisNameError + if `model.axes` contains axes with duplicate names + """ + return CoordSystem( + name=model.name, + axes=tuple(Axis.try_from_model(axis) for axis in model.axes), + virtual=False, + ) + + @property + def num_axes(self) -> int: + return len(self.axes) + + @property + def axes_names(self) -> tuple[str, ...]: + """Get axes' names""" + return tuple([ax.name for ax in self.axes]) + + @property + def axes_types(self) -> tuple[str, ...]: + """Get axes' types""" + return tuple([ax.type for ax in self.axes]) + + def has_axis(self, name: str) -> bool: + """ + Check the coordinate system has an axis of the given name. + + Parameters + ---------- + name + name of the axis. + """ + return any(axis.name == name for axis in self.axes) + + def get_axis(self, name: str) -> Axis | None: + """Get the axis by name""" + for axis in self.axes: + if axis.name == name: + return axis + return None + + def get_spatial_axes(self) -> Sequence[Axis]: + return [axis for axis in self.axes if axis.type == "space"] diff --git a/src/spatialdata/transformations/_utils.py b/src/spatialdata/transformations/_utils.py index 032c89c13..e1bca3b05 100644 --- a/src/spatialdata/transformations/_utils.py +++ b/src/spatialdata/transformations/_utils.py @@ -10,16 +10,18 @@ from xarray import DataArray, Dataset, DataTree from spatialdata._types import ArrayLike +from spatialdata.literals import TRANSFORM_KEY +from spatialdata.transformations.transformations import Affine, BaseTransformation, Scale, Sequence if TYPE_CHECKING: from spatialdata._core.spatialdata import SpatialData from spatialdata.models import SpatialElement - from spatialdata.models._utils import MappingToCoordinateSystem_t - from spatialdata.transformations.transformations import Affine, BaseTransformation, Scale + + +MappingToCoordinateSystem_t = dict[str, BaseTransformation] def _get_transformations_from_dict_container(dict_container: Any) -> MappingToCoordinateSystem_t | None: - from spatialdata.models._utils import TRANSFORM_KEY if TRANSFORM_KEY in dict_container: d = dict_container[TRANSFORM_KEY] @@ -38,7 +40,6 @@ def _get_transformations(e: SpatialElement) -> MappingToCoordinateSystem_t | Non def _set_transformations_to_dict_container(dict_container: Any, transformations: MappingToCoordinateSystem_t) -> None: - from spatialdata.models._utils import TRANSFORM_KEY if TRANSFORM_KEY not in dict_container: dict_container[TRANSFORM_KEY] = {} @@ -48,7 +49,6 @@ def _set_transformations_to_dict_container(dict_container: Any, transformations: def _set_transformations_to_element(element: Any, transformations: MappingToCoordinateSystem_t) -> None: - from spatialdata.models._utils import TRANSFORM_KEY attrs = element.attrs if TRANSFORM_KEY not in attrs: @@ -88,7 +88,6 @@ def _(e: DataTree, transformations: MappingToCoordinateSystem_t) -> None: # set the transformation at the highest level and concatenate with the appropriate scale at each level dims = get_axes_names(e) - from spatialdata.transformations.transformations import Scale, Sequence old_shape: ArrayLike | None = None for i, (scale, node) in enumerate(dict(e).items()): @@ -122,7 +121,7 @@ def _(e: DataTree, transformations: MappingToCoordinateSystem_t) -> None: @_set_transformations.register(GeoDataFrame) @_set_transformations.register(DaskDataFrame) -def _(e: GeoDataFrame | GeoDataFrame, transformations: MappingToCoordinateSystem_t) -> None: +def _(e: GeoDataFrame | DaskDataFrame, transformations: MappingToCoordinateSystem_t) -> None: _set_transformations_to_element(e, transformations) # _set_transformations_to_dict_container(e.attrs, transformations) @@ -134,7 +133,6 @@ def _(e: DataArray) -> MappingToCoordinateSystem_t | None: @_get_transformations.register(DataTree) def _(e: DataTree) -> MappingToCoordinateSystem_t | None: - from spatialdata.models._utils import TRANSFORM_KEY if TRANSFORM_KEY in e.attrs: raise ValueError( @@ -212,7 +210,7 @@ def _(data: DataArray) -> DataArray: @compute_coordinates.register(DataTree) def _(data: DataTree) -> DataTree: - from spatialdata.models import get_axes_names + from spatialdata.models._utils import get_axes_names spatial_coords = [ax for ax in get_axes_names(data) if ax in ["x", "y", "z"]] img_name = list(data["scale0"].data_vars.keys())[0] @@ -295,3 +293,15 @@ def convert_transformations_to_affine(sdata: SpatialData, coordinate_system: str if not isinstance(t, Affine): affine = _get_affine_for_element(element, t) set_transformation(element, transformation=affine, to_coordinate_system=coordinate_system) + + +# added this code as part of a refactoring to catch errors earlier +def _validate_mapping_to_coordinate_system_type(transformations: MappingToCoordinateSystem_t | None) -> None: + if not ( + transformations is None + or isinstance(transformations, dict) + and all(isinstance(k, str) and isinstance(v, BaseTransformation) for k, v in transformations.items()) + ): + raise TypeError( + f"Transform must be of type {MappingToCoordinateSystem_t} or None, but is of type {type(transformations)}." + ) diff --git a/src/spatialdata/transformations/ngff/_utils.py b/src/spatialdata/transformations/ngff/_utils.py index 07dff152c..8652252c4 100644 --- a/src/spatialdata/transformations/ngff/_utils.py +++ b/src/spatialdata/transformations/ngff/_utils.py @@ -2,7 +2,7 @@ import copy -from spatialdata.models import C, X, Y, Z +from spatialdata.literals import C, X, Y, Z from spatialdata.transformations.ngff.ngff_coordinate_system import NgffAxis, NgffCoordinateSystem __all__ = "get_default_coordinate_system" diff --git a/src/spatialdata/transformations/ngff/ngff_transformations.py b/src/spatialdata/transformations/ngff/ngff_transformations.py index 1758159bf..9caf84584 100644 --- a/src/spatialdata/transformations/ngff/ngff_transformations.py +++ b/src/spatialdata/transformations/ngff/ngff_transformations.py @@ -8,6 +8,7 @@ import numpy as np from spatialdata._types import ArrayLike +from spatialdata.literals import C, X, Y, Z from spatialdata.transformations.ngff.ngff_coordinate_system import NgffCoordinateSystem __all__ = [ @@ -366,8 +367,6 @@ def _affine_matrix_from_input_and_output_axes( output_axes the output axes. """ - from spatialdata.models import C, X, Y, Z - assert all(ax in (X, Y, Z, C) for ax in input_axes) assert all(ax in (X, Y, Z, C) for ax in output_axes) m = np.zeros((len(output_axes) + 1, len(input_axes) + 1)) diff --git a/src/spatialdata/transformations/operations.py b/src/spatialdata/transformations/operations.py index 9091cec38..9b9abfcb3 100644 --- a/src/spatialdata/transformations/operations.py +++ b/src/spatialdata/transformations/operations.py @@ -1,21 +1,23 @@ from __future__ import annotations import contextlib -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, Literal, overload import numpy as np from dask.dataframe import DataFrame as DaskDataFrame from geopandas import GeoDataFrame +from skimage.transform import estimate_transform from spatialdata._logging import logger +from spatialdata.literals import DEFAULT_COORDINATE_SYSTEM from spatialdata.transformations._utils import _get_transformations, _set_transformations +from spatialdata.transformations.transformations import Affine, BaseTransformation, Sequence if TYPE_CHECKING: import networkx as nx from spatialdata._core.spatialdata import SpatialData - from spatialdata.models._utils import SpatialElement - from spatialdata.transformations.transformations import Affine, BaseTransformation + from spatialdata.models import SpatialElement def set_transformation( @@ -49,7 +51,6 @@ def set_transformation( to the SpatialData object, and the SpatialData object needs to be backed. """ - from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM from spatialdata.transformations import BaseTransformation if write_to_sdata is None: @@ -88,6 +89,19 @@ def set_transformation( write_to_sdata.write_transformations(element_name=element_name) +# these signatures help the static type checker understand the correct return type based on the inputs +@overload +def get_transformation( + element: SpatialElement, to_coordinate_system: str | None = None, get_all: Literal[False] = False +) -> BaseTransformation: ... + + +@overload +def get_transformation( + element: SpatialElement, to_coordinate_system: None = None, get_all: Literal[True] = True +) -> dict[str, BaseTransformation]: ... + + def get_transformation( element: SpatialElement, to_coordinate_system: str | None = None, get_all: bool = False ) -> BaseTransformation | dict[str, BaseTransformation]: @@ -112,7 +126,6 @@ def get_transformation( The transformation, if `to_coordinate_system` is not None, otherwise a dictionary of transformations to all the coordinate systems. """ - from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM transformations = _get_transformations(element) assert isinstance(transformations, dict) @@ -159,7 +172,6 @@ def remove_transformation( If not None, the element needs to belong to the SpatialData object, and the SpatialData object needs to be backed. """ - from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM if write_to_sdata is None: if remove_all is False: @@ -378,11 +390,8 @@ def get_transformation_between_landmarks( >>> references_coords = PointsModel(points_reference) >>> transformation = get_transformation_between_landmarks(references_coords, moving_coords) """ - from skimage.transform import estimate_transform - from spatialdata import transform from spatialdata.models import get_axes_names - from spatialdata.transformations.transformations import Affine, Sequence assert get_axes_names(references_coords) == ("x", "y") assert get_axes_names(moving_coords) == ("x", "y") diff --git a/src/spatialdata/transformations/transformations.py b/src/spatialdata/transformations/transformations.py index 32b81b1be..400b69db9 100644 --- a/src/spatialdata/transformations/transformations.py +++ b/src/spatialdata/transformations/transformations.py @@ -1,7 +1,7 @@ from __future__ import annotations from abc import ABC, abstractmethod -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, override from warnings import warn import numpy as np @@ -9,6 +9,26 @@ import xarray as xr from xarray import DataArray +from spatialdata.literals import DEFAULT_COORDINATE_SYSTEM +from spatialdata.models import get_axes_names +from spatialdata.transformations._graph import ( + AffineEdge, + Axis, + BaseTransformationEdge, + CoordSystem, + CsGen, + IdentityEdge, + IncompatibleCoordSystemsError, + MapAxisEdge, + MissingAxisError, + NGFFCompatibilityError, + ScaleEdge, + SequenceEdge, + TranslationEdge, + UnmappedAxisError, +) +from spatialdata.transformations._graph.exceptions import IncompatibleTransformationOutputCoordinateSystemError +from spatialdata.transformations.ngff._utils import get_default_coordinate_system from spatialdata.transformations.ngff.ngff_coordinate_system import NgffCoordinateSystem, _get_spatial_axes from spatialdata.transformations.ngff.ngff_transformations import ( NgffAffine, @@ -23,7 +43,7 @@ if TYPE_CHECKING: from spatialdata._types import ArrayLike, ListOrNDArrayFloating from spatialdata.models import SpatialElement - from spatialdata.models._utils import ValidAxis_t + from spatialdata.transformations._graph.types import ValidAxis_t TRANSFORMATIONS_MAP: dict[type[NgffBaseTransformation], type[BaseTransformation]] = {} @@ -38,7 +58,7 @@ def validate_axes(axes: tuple[ValidAxis_t, ...]) -> None: This function is to allow to call validate_axes() from this file in multiple places while avoiding circular imports. """ - from spatialdata.models._utils import validate_axes + from spatialdata.transformations._graph.utils import validate_axes validate_axes(axes) @@ -63,6 +83,22 @@ def _repr_indent(self, indent: int = 0) -> str: def __repr__(self) -> str: return self._repr_indent(0) + @abstractmethod + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + """ + Get the resulting output axes for the transformation, given a tuple of input axes. + + Parameters + ---------- + input_axes + The input axes. + + Returns + ------- + The output axes after applying the transformation. + """ + pass + @classmethod @abstractmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: @@ -85,6 +121,126 @@ def to_ngff( ) -> NgffBaseTransformation: pass + def _to_ngff_transformation_edge_creating_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system_name: str, + output_coordinate_system_virtual: bool = False, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + """ + Convert the transformation to an NGFF transformation edge, creating its output coordinate system + + Parameters + ---------- + input_coordinate_system + The input coordinate system. + output_coordinate_system_name + The name for the output coordinate system. + output_coordinate_system_virtual + Whether the output coordinate system is virtual + transformation_edge_name + Optional name for the transformation edge. + + Returns + ------- + A subclass of BaseTransformationEdge representing an equivalent NGFF coordinate transformation. + """ + + output_coordinate_system = self._get_ngff_output_coordinate_system( + input_coordinate_system=input_coordinate_system, + output_coordinate_system_name=output_coordinate_system_name, + output_coordinate_system_virtual=output_coordinate_system_virtual, + ) + + return self._define_ngff_transformation_edge_given_valid_output_cs( + input_coordinate_system=input_coordinate_system, + output_coordinate_system=output_coordinate_system, + transformation_edge_name=transformation_edge_name, + ) + + @abstractmethod + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + pass + + def _to_ngff_transformation_edge_existing_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + """ + Convert the transformation to an NGFF transformation edge, creating its output coordinate system + + Parameters + ---------- + input_coordinate_system + The input coordinate system. + output_coordinate_system + Existing coordinate system to use as output + transformation_edge_name + Optional name for the transformation edge. + + Returns + ------- + A subclass of BaseTransformationEdge representing an equivalent NGFF coordinate transformation. + """ + + self.assert_input_output_cs_compatible(input_cs=input_coordinate_system, output_cs=output_coordinate_system) + + return self._define_ngff_transformation_edge_given_valid_output_cs( + input_coordinate_system=input_coordinate_system, + output_coordinate_system=output_coordinate_system, + ) + + def assert_input_output_cs_compatible(self, input_cs: CoordSystem, output_cs: CoordSystem) -> None: + + expected_cs = self._get_ngff_output_coordinate_system( + input_coordinate_system=input_cs, + output_coordinate_system_name="dummy", + ) + + if set(expected_cs.axes) != set(output_cs.axes): + raise IncompatibleTransformationOutputCoordinateSystemError( + input_coordinate_system_name=input_cs.name, + output_coordinate_system_name=output_cs.name, + transformation_class_name=self.__class__.__name__, + ) + + def _get_ngff_output_coordinate_system( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system_name: str, + output_coordinate_system_virtual: bool = False, + ) -> CoordSystem: + """ + Get the NGFF output coordinate system for a transformation. + + Parameters + ---------- + input_coordinate_system + The input coordinate system. + output_coordinate_system_name + The name for the output coordinate system. + + Returns + ------- + The output coordinate system. + """ + + output_axes_names = self._get_resulting_output_axes(input_coordinate_system.axes_names) + + output_axes = tuple(Axis.from_spatialdata_axis_name(name) for name in output_axes_names) + return CoordSystem( + name=output_coordinate_system_name, axes=output_axes, virtual=output_coordinate_system_virtual + ) + def _get_default_coordinate_system( self, axes: tuple[ValidAxis_t, ...], @@ -92,7 +248,6 @@ def _get_default_coordinate_system( name: str | None = None, default_to_global: bool = False, ) -> NgffCoordinateSystem: - from spatialdata.transformations.ngff._utils import get_default_coordinate_system cs = get_default_coordinate_system(axes) if unit is not None: @@ -102,8 +257,6 @@ def _get_default_coordinate_system( if name is not None: cs.name = name elif default_to_global: - from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM - cs.name = DEFAULT_COORDINATE_SYSTEM return cs @@ -230,6 +383,9 @@ def _repr_transformation_description(self, indent: int = 0) -> str: def _transform_coordinates(self, data: DataArray) -> DataArray: return data + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + return input_axes + @classmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: assert isinstance(t, NgffIdentity) @@ -252,6 +408,17 @@ def to_ngff( ngff_transformation = NgffIdentity(input_coordinate_system=input_cs, output_coordinate_system=output_cs) return ngff_transformation + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + return IdentityEdge( + input=input_coordinate_system, output=output_coordinate_system, name=transformation_edge_name + ) + def __eq__(self, other: Any) -> bool: return isinstance(other, Identity) @@ -272,7 +439,7 @@ class MapAxis(BaseTransformation): def __init__(self, map_axis: dict[ValidAxis_t, ValidAxis_t]) -> None: # to avoid circular imports - from spatialdata.models import validate_axis_name + from spatialdata.transformations._graph.utils import validate_axis_name assert isinstance(map_axis, dict) for des_ax, src_ax in map_axis.items(): @@ -319,7 +486,7 @@ def _repr_transformation_description(self, indent: int = 0) -> str: def _transform_coordinates(self, data: DataArray) -> DataArray: self._xarray_coords_validate_axes(data) data_input_axes = self._xarray_coords_get_coords(data) - data_output_axes = _get_current_output_axes(self, data_input_axes) + data_output_axes = self._get_resulting_output_axes(input_axes=data_input_axes) transformed = [] for ax in data_output_axes: @@ -333,6 +500,20 @@ def _transform_coordinates(self, data: DataArray) -> DataArray: to_return = self._xarray_coords_reorder_axes(to_return) return to_return + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + map_axis_input_axes = set(self.map_axis.values()) + set(self.map_axis.keys()) + to_return = [] + for ax in input_axes: + if ax not in map_axis_input_axes: + assert ax not in to_return + to_return.append(ax) + else: + mapped = [ax_out for ax_out, ax_in in self.map_axis.items() if ax_in == ax] + assert all(ax_out not in to_return for ax_out in mapped) + to_return.extend(mapped) + return tuple(to_return) + @classmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: assert isinstance(t, NgffMapAxis) @@ -357,6 +538,46 @@ def to_ngff( ) return ngff_transformation + def _to_ngff_input_output_mapping( + self, input_coordinate_system: CoordSystem, output_coordinate_system: CoordSystem + ) -> dict[Axis, Axis]: + + input_to_output_axes_mapping = { + Axis.from_spatialdata_axis_name(input_axis_name): Axis.from_spatialdata_axis_name(output_axis_name) + for output_axis_name, input_axis_name in self.map_axis.items() + } + + for input_cs_axis in input_coordinate_system.axes: + if (input_cs_axis not in input_to_output_axes_mapping) and (input_cs_axis in output_coordinate_system.axes): + input_to_output_axes_mapping[input_cs_axis] = input_cs_axis + + return input_to_output_axes_mapping + + @override + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + input_to_output_axes_mapping = self._to_ngff_input_output_mapping( + input_coordinate_system, output_coordinate_system + ) + + try: + map_axis_edge = MapAxisEdge( + input=input_coordinate_system, + output=output_coordinate_system, + name=transformation_edge_name, + input_to_output=input_to_output_axes_mapping, + ) + except (IncompatibleCoordSystemsError, MissingAxisError, UnmappedAxisError) as e: + raise NGFFCompatibilityError( + f"Cannot map the following transformation to an NGFF transformation Edge.\n{repr(self)}" + ) from e + return map_axis_edge + def __eq__(self, other: Any) -> bool: return isinstance(other, MapAxis) and self.map_axis == other.map_axis @@ -410,6 +631,9 @@ def _transform_coordinates(self, data: DataArray) -> DataArray: to_return = self._xarray_coords_reorder_axes(transformed) return to_return + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + return input_axes + @classmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: assert isinstance(t, NgffTranslation) @@ -440,6 +664,23 @@ def to_ngff( ) return ngff_transformation + @override + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + translation_vector = self.to_translation_vector(axes=output_coordinate_system.axes_names) + + return TranslationEdge( + input=input_coordinate_system, + output=output_coordinate_system, + translation=translation_vector, + name=transformation_edge_name, + ) + def __eq__(self, other: Any) -> bool: return ( isinstance(other, Translation) @@ -496,6 +737,9 @@ def _transform_coordinates(self, data: DataArray) -> DataArray: to_return = self._xarray_coords_reorder_axes(transformed) return to_return + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + return input_axes + @classmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: assert isinstance(t, NgffScale) @@ -523,6 +767,23 @@ def to_ngff( ) return ngff_transformation + @override + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + scale_vector = self.to_scale_vector(input_coordinate_system.axes_names) + + return ScaleEdge( + input=input_coordinate_system, + output=output_coordinate_system, + scale=scale_vector, + name=transformation_edge_name, + ) + def __eq__(self, other: Any) -> bool: return isinstance(other, Scale) and np.allclose(self.scale, other.scale) and self.axes == other.axes @@ -782,7 +1043,7 @@ def _repr_transformation_description(self, indent: int = 0) -> str: def _transform_coordinates(self, data: DataArray) -> DataArray: self._xarray_coords_validate_axes(data) data_input_axes = self._xarray_coords_get_coords(data) - data_output_axes = _get_current_output_axes(self, data_input_axes) + data_output_axes = self._get_resulting_output_axes(data_input_axes) matrix = self.to_affine_matrix(data_input_axes, data_output_axes) transformed = (matrix @ np.vstack((data.data.T, np.ones(len(data))))).T[:, :-1] to_return = DataArray(transformed, coords={"points": data.coords["points"], "dim": list(data_output_axes)}) @@ -799,6 +1060,27 @@ def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: output_axes = tuple(t.output_coordinate_system.axes_names) return Affine(matrix=t.affine, input_axes=input_axes, output_axes=output_axes) + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + to_return = [] + add_affine_output_axes = False + for ax in input_axes: + if ax not in self.input_axes: + assert ax not in to_return + to_return.append(ax) + else: + add_affine_output_axes = True + if add_affine_output_axes: + for ax in self.output_axes: + if ax not in to_return: + to_return.append(ax) + else: + raise ValueError( + f"Trying to query an invalid representation of an affine matrix: the ax {ax} is not " + f"an input axis of the affine matrix but it appears both as output as input of the " + f"matrix representation being queried" + ) + return tuple(to_return) + def to_ngff( self, input_axes: tuple[ValidAxis_t, ...], @@ -819,6 +1101,27 @@ def to_ngff( ) return ngff_transformation + @override + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + affine_matrix = self.to_affine_matrix( + input_axes=input_coordinate_system.axes_names, + output_axes=output_coordinate_system.axes_names, + ) + + return AffineEdge( + input=input_coordinate_system, + output=output_coordinate_system, + linear=affine_matrix[:-1, :-1], + translation=affine_matrix[:-1, -1], + name=transformation_edge_name, + ) + def __eq__(self, other: Any) -> bool: if not isinstance(other, Affine): return False @@ -848,7 +1151,7 @@ def _to_affine_matrix_wrapper( raise ValueError("Input axes must be a subset of output axes.") current_input_axes = input_axes - current_output_axes = _get_current_output_axes(self.transformations[0], current_input_axes) + current_output_axes = self.transformations[0]._get_resulting_output_axes(input_axes=current_input_axes) m = self.transformations[0].to_affine_matrix(current_input_axes, current_output_axes) if DEBUG_SEQUENCE: print(f"# 0: current_input_axes = {current_input_axes}, current_output_axes = {current_output_axes}") @@ -856,7 +1159,7 @@ def _to_affine_matrix_wrapper( print() for i, t in enumerate(self.transformations[1:]): current_input_axes = current_output_axes - current_output_axes = _get_current_output_axes(t, current_input_axes) + current_output_axes = t._get_resulting_output_axes(input_axes=current_input_axes) if DEBUG_SEQUENCE: print( f"# {i + 1}: current_input_axes = {current_input_axes}, current_output_axes = {current_output_axes}" @@ -917,6 +1220,13 @@ def _transform_coordinates(self, data: DataArray) -> DataArray: self._xarray_coords_validate_axes(data) return data + def _get_resulting_output_axes(self, input_axes: tuple[ValidAxis_t, ...]) -> tuple[ValidAxis_t, ...]: + + input_axes_this_iteration = input_axes + for t in self.transformations: + input_axes_this_iteration = t._get_resulting_output_axes(input_axes_this_iteration) + return input_axes_this_iteration + @classmethod def _from_ngff(cls, t: NgffBaseTransformation) -> BaseTransformation: assert isinstance(t, NgffSequence) @@ -939,7 +1249,7 @@ def to_ngff( converted_transformations = [] latest_input_axes = input_axes for t in self.transformations: - latest_output_axes = _get_current_output_axes(t, latest_input_axes) + latest_output_axes = t._get_resulting_output_axes(input_axes=latest_input_axes) converted_transformations.append( t.to_ngff( input_axes=latest_input_axes, @@ -956,63 +1266,45 @@ def to_ngff( ) return ngff_transformation + @override + def _define_ngff_transformation_edge_given_valid_output_cs( + self, + input_coordinate_system: CoordSystem, + output_coordinate_system: CoordSystem, + transformation_edge_name: str | None = None, + ) -> BaseTransformationEdge: + + cs_gen = CsGen(base_name=f"Sequence_from_{input_coordinate_system.name}_to_{output_coordinate_system.name}") + + individual_edges = [] + current_input_cs = input_coordinate_system + for t in self.transformations[:-1]: + current_edge = t._to_ngff_transformation_edge_creating_output_cs( + input_coordinate_system=current_input_cs, + output_coordinate_system_name=cs_gen.get_unique_name(), + output_coordinate_system_virtual=True, + ) + individual_edges.append(current_edge) + current_input_cs = current_edge.output + + last_edge = self.transformations[-1]._define_ngff_transformation_edge_given_valid_output_cs( + input_coordinate_system=current_input_cs, + output_coordinate_system=output_coordinate_system, + ) + individual_edges.append(last_edge) + + return SequenceEdge(name=transformation_edge_name, transformations=individual_edges) + def __eq__(self, other: Any) -> bool: if not isinstance(other, Sequence): return False return self.transformations == other.transformations -def _get_current_output_axes( - transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...] -) -> tuple[ValidAxis_t, ...]: - if isinstance(transformation, Identity | Translation | Scale): - return input_axes - elif isinstance(transformation, MapAxis): - map_axis_input_axes = set(transformation.map_axis.values()) - set(transformation.map_axis.keys()) - to_return = [] - for ax in input_axes: - if ax not in map_axis_input_axes: - assert ax not in to_return - to_return.append(ax) - else: - mapped = [ax_out for ax_out, ax_in in transformation.map_axis.items() if ax_in == ax] - assert all(ax_out not in to_return for ax_out in mapped) - to_return.extend(mapped) - return tuple(to_return) - elif isinstance(transformation, Affine): - to_return = [] - add_affine_output_axes = False - for ax in input_axes: - if ax not in transformation.input_axes: - assert ax not in to_return - to_return.append(ax) - else: - add_affine_output_axes = True - if add_affine_output_axes: - for ax in transformation.output_axes: - if ax not in to_return: - to_return.append(ax) - else: - raise ValueError( - f"Trying to query an invalid representation of an affine matrix: the ax {ax} is not " - f"an input axis of the affine matrix but it appears both as output as input of the " - f"matrix representation being queried" - ) - return tuple(to_return) - elif isinstance(transformation, Sequence): - for t in transformation.transformations: - input_axes = _get_current_output_axes(t, input_axes) - return input_axes - else: - raise ValueError("Unknown transformation type.") - - def _get_affine_for_element(element: SpatialElement, transformation: BaseTransformation) -> Affine: - from spatialdata.models import get_axes_names input_axes = get_axes_names(element) - output_axes = _get_current_output_axes(transformation, input_axes) + output_axes = transformation._get_resulting_output_axes(input_axes) matrix = transformation.to_affine_matrix(input_axes=input_axes, output_axes=output_axes) return Affine(matrix, input_axes=input_axes, output_axes=output_axes) @@ -1060,7 +1352,7 @@ def _validate_square_affine_for_decomposition( If the linear part of the affine has a large condition number, in which case the decomposition may be numerically inaccurate. """ - output_axes = _get_current_output_axes(transformation=transformation, input_axes=input_axes) + output_axes = transformation._get_resulting_output_axes(input_axes=input_axes) if set(input_axes) != set(output_axes): raise ValueError("The transformation should leave the set of input axes unmodified.") # the axes may come out in a different order than input_axes; querying in input_axes order makes the matrix diff --git a/tests/core/operations/test_rasterize.py b/tests/core/operations/test_rasterize.py index 53d362d15..48fa4a9bf 100644 --- a/tests/core/operations/test_rasterize.py +++ b/tests/core/operations/test_rasterize.py @@ -16,9 +16,9 @@ from spatialdata._core.operations.rasterize import rasterize from spatialdata._core.query.relational_query import get_element_instances from spatialdata._utils import get_pyramid_levels -from spatialdata.models import PointsModel, ShapesModel, TableModel, get_axes_names -from spatialdata.models._utils import get_spatial_axes -from spatialdata.transformations import MapAxis +from spatialdata.models import PointsModel, ShapesModel, TableModel +from spatialdata.models._utils import get_axes_names +from spatialdata.transformations import MapAxis, get_spatial_axes from tests.conftest import _get_images, _get_labels diff --git a/tests/core/operations/test_rasterize_bins.py b/tests/core/operations/test_rasterize_bins.py index afd807657..ba5479bc0 100644 --- a/tests/core/operations/test_rasterize_bins.py +++ b/tests/core/operations/test_rasterize_bins.py @@ -22,7 +22,7 @@ from spatialdata._core.spatialdata import SpatialData from spatialdata._logging import logger from spatialdata._types import ArrayLike -from spatialdata.models.models import ( +from spatialdata.models import ( Image2DModel, PointsModel, ShapesModel, diff --git a/tests/core/operations/test_transform.py b/tests/core/operations/test_transform.py index ef307ac7b..3286fe312 100644 --- a/tests/core/operations/test_transform.py +++ b/tests/core/operations/test_transform.py @@ -16,7 +16,8 @@ from spatialdata._core.data_extent import are_extents_equal, get_extent from spatialdata._core.spatialdata import SpatialData from spatialdata._utils import disable_dask_tune_optimization, unpad_raster -from spatialdata.models import Image2DModel, PointsModel, ShapesModel, get_axes_names +from spatialdata.models import Image2DModel, PointsModel, ShapesModel +from spatialdata.models._utils import get_axes_names from spatialdata.transformations.operations import ( align_elements_using_landmarks, get_transformation, diff --git a/tests/core/operations/test_vectorize.py b/tests/core/operations/test_vectorize.py index b0368a12a..4a962c655 100644 --- a/tests/core/operations/test_vectorize.py +++ b/tests/core/operations/test_vectorize.py @@ -15,7 +15,7 @@ to_polygons, ) from spatialdata.datasets import blobs -from spatialdata.models.models import ShapesModel +from spatialdata.models import ShapesModel from spatialdata.testing import assert_elements_are_identical diff --git a/tests/core/query/test_relational_query.py b/tests/core/query/test_relational_query.py index 07d232559..2e7aaf9dd 100644 --- a/tests/core/query/test_relational_query.py +++ b/tests/core/query/test_relational_query.py @@ -18,7 +18,7 @@ get_element_annotators, join_spatialelement_table, ) -from spatialdata.models.models import TableModel +from spatialdata.models import TableModel from spatialdata.testing import assert_anndata_equal, assert_geodataframe_equal diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index 3b6ff3581..27ca2093b 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -23,6 +23,7 @@ ) from spatialdata._core.spatialdata import SpatialData from spatialdata.models import ( + ATTRS_KEY, Image2DModel, Image3DModel, Labels2DModel, @@ -31,7 +32,6 @@ ShapesModel, TableModel, ) -from spatialdata.models.models import ATTRS_KEY from spatialdata.testing import assert_spatial_data_objects_are_identical from spatialdata.transformations import Identity, MapAxis, set_transformation from tests.conftest import _make_points, _make_squares diff --git a/tests/core/test_centroids.py b/tests/core/test_centroids.py index f8c8be1da..5c1f05ae6 100644 --- a/tests/core/test_centroids.py +++ b/tests/core/test_centroids.py @@ -10,7 +10,8 @@ from spatialdata._core.centroids import get_centroids from spatialdata._core.query.relational_query import get_element_instances -from spatialdata.models import Labels2DModel, Labels3DModel, PointsModel, TableModel, get_axes_names +from spatialdata.models import Labels2DModel, Labels3DModel, PointsModel, TableModel +from spatialdata.models._utils import get_axes_names from spatialdata.transformations import Affine, Identity, get_transformation, set_transformation RNG = default_rng(42) diff --git a/tests/io/test_ngff_06.py b/tests/io/test_ngff_06.py new file mode 100644 index 000000000..222bae81e --- /dev/null +++ b/tests/io/test_ngff_06.py @@ -0,0 +1,54 @@ +from __future__ import annotations + +import numpy as np +from ome_zarr import OMEZarrImage, OMEZarrMultiscale +from ome_zarr_models.v06.coordinate_transforms import Scale as ModelsScale +from ome_zarr_models.v06.coordinate_transforms import Sequence as ModelsSequence +from ome_zarr_models.v06.coordinate_transforms import Translation as ModelsTranslation + +from spatialdata._io.io_raster import try_parse_ngff06_multiscale +from tests.conftest import SEED + + +def test_parse_multiscale(): + full_shape = (256, 256) + data = np.random.default_rng(seed=SEED).random(full_shape) + omz_img = OMEZarrImage(data=data, axes="yx") + scale_factors = (2, 4, 8, 16) + ms = OMEZarrMultiscale(image=omz_img, scale_factors=scale_factors) + + data_tree, _transforms = try_parse_ngff06_multiscale(ms) + for xr_scale_node, ngff_scale, ngff_meta in zip( + data_tree.children.values(), ms.images, ms.metadata.datasets, strict=True + ): + xr_scale = xr_scale_node["image"] + assert xr_scale.shape == ngff_scale.data.shape + + coords_x = xr_scale.coords["x"] + coords_y = xr_scale.coords["y"] + + seq = ngff_meta.coordinateTransformations[0] + assert isinstance(seq, ModelsSequence) + scale = seq.transformations[0] + assert isinstance(scale, ModelsScale) + translate = seq.transformations[1] + assert isinstance(translate, ModelsTranslation) + + start_indices = (0, 0) + start = (coords_x[0], coords_y[0]) + + end_indices = (xr_scale.shape[0] - 1, xr_scale.shape[1] - 1) + end = (coords_x[-1], coords_y[-1]) + + for indices, point in [(start_indices, start), (end_indices, end)]: + expected = np.asarray(scale.scale) * indices + translate.translation + assert np.allclose(expected, point) + + sliced_data_tree = data_tree.sel(x=slice(0, 256), y=slice(0, 256), method="nearest") + assert sliced_data_tree.equals(data_tree) + + quarter_tree = data_tree.sel(x=slice(0, 128), y=slice(0, 128), method="nearest") + for xr_scale_node, scale_factor in zip(quarter_tree.children.values(), [1, *scale_factors], strict=True): + xr_scale = xr_scale_node["image"] + expected_shape = np.asarray(full_shape) / 2 / scale_factor + np.testing.assert_equal(xr_scale.shape, expected_shape) diff --git a/tests/models/test_models.py b/tests/models/test_models.py index 4ee3210b7..d39c17de4 100644 --- a/tests/models/test_models.py +++ b/tests/models/test_models.py @@ -31,14 +31,7 @@ from spatialdata._core.validation import ValidationError from spatialdata._types import ArrayLike from spatialdata.config import LARGE_CHUNK_THRESHOLD_BYTES -from spatialdata.models import get_table_keys -from spatialdata.models._utils import ( - force_2d, - points_dask_dataframe_to_geopandas, - points_geopandas_to_dask_dataframe, - validate_axis_name, -) -from spatialdata.models.models import ( +from spatialdata.models import ( Image2DModel, Image3DModel, Labels2DModel, @@ -47,10 +40,17 @@ RasterSchema, ShapesModel, TableModel, - get_axes_names, get_model, + get_table_keys, +) +from spatialdata.models._utils import ( + force_2d, + get_axes_names, + points_dask_dataframe_to_geopandas, + points_geopandas_to_dask_dataframe, ) from spatialdata.testing import assert_elements_are_identical +from spatialdata.transformations._graph.utils import validate_axis_name from spatialdata.transformations._utils import _set_transformations from spatialdata.transformations.operations import get_transformation, set_transformation from spatialdata.transformations.transformations import Identity, Scale diff --git a/tests/transformations/graph/__init__.py b/tests/transformations/graph/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/transformations/graph/conftest.py b/tests/transformations/graph/conftest.py new file mode 100644 index 000000000..9a6f21752 --- /dev/null +++ b/tests/transformations/graph/conftest.py @@ -0,0 +1,255 @@ +""" +Global pytest fixtures for graph transformations. +""" + +from __future__ import annotations + +import pytest + +from spatialdata.literals import NgffAxisType +from spatialdata.transformations._graph import Axis, CoordSystem + + +# ---------------------------------------------------------------------- +# Simple axes – one per common name / type combination +# ---------------------------------------------------------------------- +@pytest.fixture +def axis_x_space() -> Axis: + """Axis(name='x', type=SPACE)""" + return Axis(name="x", type=NgffAxisType.SPACE) + + +@pytest.fixture +def axis_y_space() -> Axis: + """Axis(name='y', type=SPACE)""" + return Axis(name="y", type=NgffAxisType.SPACE) + + +@pytest.fixture +def axis_z_space() -> Axis: + """Axis(name='z', type=SPACE)""" + return Axis(name="z", type=NgffAxisType.SPACE) + + +@pytest.fixture +def axis_c_channel() -> Axis: + """Axis(name='c', type=CHANNEL)""" + return Axis(name="c", type=NgffAxisType.CHANNEL) + + +# ---------------------------------------------------------------------- +# Axes with optional metadata +# ---------------------------------------------------------------------- +@pytest.fixture +def axis_x_space_with_unit() -> Axis: + """Axis(name='x', type=SPACE, unit='micrometer')""" + return Axis(name="x", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_y_space_with_unit() -> Axis: + """Axis(name='y', type=SPACE, unit='micrometer')""" + return Axis(name="y", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_z_space_with_unit() -> Axis: + """Axis(name='z', type=SPACE, unit='micrometer')""" + return Axis(name="z", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_a_space_with_unit() -> Axis: + """Axis(name='a', type=SPACE, unit='micrometer')""" + return Axis(name="a", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_b_space_with_unit() -> Axis: + """Axis(name='b', type=SPACE, unit='micrometer')""" + return Axis(name="b", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_c_space_with_unit() -> Axis: + """Axis(name='x', type=SPACE, unit='micrometer')""" + return Axis(name="c", type=NgffAxisType.SPACE, unit="micrometer") + + +@pytest.fixture +def axis_x_space_with_long_name() -> Axis: + """Axis(name='x', type=SPACE, long_name='X axis')""" + return Axis(name="x", type=NgffAxisType.SPACE, long_name="X axis") + + +@pytest.fixture +def axis_x_space_full() -> Axis: + """Axis(name='x', type=SPACE, unit='micrometer', long_name='X axis')""" + return Axis( + name="x", + type=NgffAxisType.SPACE, + unit="micrometer", + long_name="X axis", + ) + + +# ---------------------------------------------------------------------- +# Composite groups that appear frequently in the code base +# ---------------------------------------------------------------------- +@pytest.fixture +def axes_xy(axis_x_space, axis_y_space) -> tuple[Axis, ...]: + """Typical 2‑D spatial system: ('x', 'y').""" + return axis_x_space, axis_y_space + + +@pytest.fixture +def axes_zyx(axis_z_space, axis_y_space, axis_x_space) -> tuple[Axis, ...]: + """Typical 3‑D spatial system in reverse order: ('z', 'y', 'x').""" + return axis_z_space, axis_y_space, axis_x_space + + +@pytest.fixture +def axes_czyx(axis_c_channel, axis_z_space, axis_y_space, axis_x_space) -> tuple[Axis, ...]: + """Channel plus full 3‑D spatial axes: ('c', 'z', 'y', 'x').""" + return axis_c_channel, axis_z_space, axis_y_space, axis_x_space + + +@pytest.fixture +def axes_cyx(axis_c_channel, axis_y_space, axis_x_space) -> tuple[Axis, ...]: + """Channel plus 2‑D spatial axes: ('c', 'y', 'x').""" + return axis_c_channel, axis_y_space, axis_x_space + + +# ---------------------------------------------------------------------- +# CoordSystem Fixtures +# ---------------------------------------------------------------------- + + +@pytest.fixture +def x_cs(axis_x_space_with_unit) -> CoordSystem: + return CoordSystem(name="x", axes=(axis_x_space_with_unit,), virtual=False) + + +@pytest.fixture +def y_cs(axis_y_space_with_unit) -> CoordSystem: + return CoordSystem(name="y", axes=(axis_y_space_with_unit,), virtual=False) + + +@pytest.fixture +def z_cs(axis_z_space_with_unit) -> CoordSystem: + return CoordSystem(name="z", axes=(axis_z_space_with_unit,), virtual=False) + + +@pytest.fixture +def c_cs(axis_c_channel) -> CoordSystem: + return CoordSystem(name="c", axes=(axis_c_channel,), virtual=False) + + +@pytest.fixture +def xy_cs(axis_x_space_with_unit, axis_y_space_with_unit) -> CoordSystem: + return CoordSystem( + name="xy", + axes=( + axis_x_space_with_unit, + axis_y_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def yx_cs(axis_y_space_with_unit, axis_x_space_with_unit) -> CoordSystem: + return CoordSystem( + name="yx", + axes=( + axis_y_space_with_unit, + axis_x_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def xyz_cs(axis_x_space_with_unit, axis_y_space_with_unit, axis_z_space_with_unit) -> CoordSystem: + return CoordSystem( + name="xyz", + axes=( + axis_x_space_with_unit, + axis_y_space_with_unit, + axis_z_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def zyx_cs(axis_z_space_with_unit, axis_y_space_with_unit, axis_x_space_with_unit) -> CoordSystem: + return CoordSystem( + name="zyx", + axes=( + axis_z_space_with_unit, + axis_y_space_with_unit, + axis_x_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def czyx_cs(axis_c_channel, axis_z_space_with_unit, axis_y_space_with_unit, axis_x_space_with_unit) -> CoordSystem: + return CoordSystem( + name="czyx", + axes=( + axis_c_channel, + axis_z_space_with_unit, + axis_y_space_with_unit, + axis_x_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def xyc_cs(axis_x_space_with_unit, axis_y_space_with_unit, axis_c_channel) -> CoordSystem: + return CoordSystem( + name="xyc", + axes=( + axis_x_space_with_unit, + axis_y_space_with_unit, + axis_c_channel, + ), + virtual=False, + ) + + +@pytest.fixture +def cyx_cs(axis_c_channel, axis_y_space_with_unit, axis_x_space_with_unit) -> CoordSystem: + return CoordSystem( + name="cyx", + axes=( + axis_c_channel, + axis_y_space_with_unit, + axis_x_space_with_unit, + ), + virtual=False, + ) + + +@pytest.fixture +def abc_cs(axis_a_space_with_unit, axis_b_space_with_unit, axis_c_space_with_unit) -> CoordSystem: + return CoordSystem( + name="abc", axes=(axis_a_space_with_unit, axis_b_space_with_unit, axis_c_space_with_unit), virtual=False + ) + + +@pytest.fixture +def cba_cs(axis_c_space_with_unit, axis_b_space_with_unit, axis_a_space_with_unit) -> CoordSystem: + return CoordSystem( + name="cba", + axes=( + axis_c_space_with_unit, + axis_b_space_with_unit, + axis_a_space_with_unit, + ), + virtual=False, + ) diff --git a/tests/transformations/graph/edge/__init__.py b/tests/transformations/graph/edge/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/transformations/graph/edge/test_edge_parsing.py b/tests/transformations/graph/edge/test_edge_parsing.py new file mode 100644 index 000000000..bc26a44d1 --- /dev/null +++ b/tests/transformations/graph/edge/test_edge_parsing.py @@ -0,0 +1,224 @@ +# pyright: strict + +from __future__ import annotations + +import numpy as np +import ome_zarr_models.v06.coordinate_transforms as ozm06trans + +from spatialdata.transformations._graph import ( + AffineEdge, + ByDimensionEdge, + CsGen, + MapAxisEdge, + RotationEdge, + ScaleEdge, + SequenceEdge, + TranslationEdge, + parse_affine, + parse_by_dimension, + parse_identity, + parse_map_axis, + parse_project_axis, + parse_rotation, + parse_scale, + parse_sequence, + parse_translation, +) + +# The ome-zarr-models transform models used below all leave `input`/`output` at their +# default of None: those fields identify coordinate systems by name/path in the NGFF +# metadata and aren't consumed by the parse_* functions, which instead take the already +# resolved `CoordSystem` (or a `CsGen`) as separate arguments. Only the transform-specific +# fields (e.g. `translation`, `scale`, `name`) are exercised here. + +POINTS_2D = np.array([[1.0, 2.0], [3.0, 4.0]]) + + +class TestParseIdentity: + def test_parses_name(self, xy_cs): + model = ozm06trans.Identity(name="my-identity") + edge = parse_identity(model, input=xy_cs, out=xy_cs) + assert edge.name == "my-identity" + assert edge.input is xy_cs + assert edge.output is xy_cs + + def test_csgen_output_mirrors_input(self, xy_cs): + model = ozm06trans.Identity(name="my-identity") + edge = parse_identity(model, input=xy_cs, out=CsGen(base_name="gen")) + assert edge.output.axes_names == xy_cs.axes_names + assert edge.output.virtual is True + + +class TestParseTranslation: + def test_parses_translation_vector(self, xy_cs): + model = ozm06trans.Translation(name="my-translation", translation=(1.0, 2.0)) + edge = parse_translation(model, input=xy_cs, out=xy_cs) + np.testing.assert_allclose(edge.translation, [1.0, 2.0]) + assert edge.name == model.name + assert edge.input is xy_cs + assert edge.output is xy_cs + + def test_csgen_output_mirrors_input(self, xy_cs): + model = ozm06trans.Translation(translation=(1.0, 2.0)) + edge = parse_translation(model, input=xy_cs, out=CsGen(base_name="gen")) + assert edge.output.axes_names == xy_cs.axes_names + assert edge.output.virtual is True + + +class TestParseProjectAxis: + def test_parses_dropped_inputs(self, xyz_cs, xy_cs, axis_z_space_with_unit): + model = ozm06trans.ProjectAxis(name="drop-z", droppedInputs=(2,)) + edge = parse_project_axis(model, input=xyz_cs, output=xy_cs) + assert edge.dropped_inputs == {axis_z_space_with_unit} + assert edge.created_outputs == set() + assert edge.name == "drop-z" + + def test_parses_created_outputs(self, xyz_cs, xy_cs, axis_z_space_with_unit): + model = ozm06trans.ProjectAxis(name="create-z", createdOutputs=(2,)) + edge = parse_project_axis(model, input=xy_cs, output=xyz_cs) + assert edge.created_outputs == { + axis_z_space_with_unit, + } + assert edge.dropped_inputs == set() + + def test_csgen_output_computes_dimensionality(self, xyz_cs, axis_z_space_with_unit): + model = ozm06trans.ProjectAxis(droppedInputs=(2,)) + edge = parse_project_axis(model, input=xyz_cs, output=CsGen(base_name="gen")) + assert edge.output.num_axes == 2 + assert edge.dropped_inputs == { + axis_z_space_with_unit, + } + + +class TestParseScale: + def test_parses_scale_vector_and_name(self, xy_cs): + model = ozm06trans.Scale(name="my-scale", scale=(2.0, 3.0)) + edge = parse_scale(model, input=xy_cs, out=xy_cs) + assert isinstance(edge, ScaleEdge) + np.testing.assert_allclose(edge.scale, [2.0, 3.0]) + assert edge.name == "my-scale" + + def test_csgen_output_mirrors_input(self, xy_cs): + model = ozm06trans.Scale(scale=(2.0, 3.0)) + edge = parse_scale(model, input=xy_cs, out=CsGen(base_name="gen")) + assert edge.output.axes_names == xy_cs.axes_names + assert edge.output.virtual is True + + +class TestParseMapAxis: + def test_parses_with_explicit_output(self, xy_cs, yx_cs): + model = ozm06trans.MapAxis(name="my-map-axis", mapAxis=(1, 0)) + edge = parse_map_axis(model, input=xy_cs, out=yx_cs) + assert isinstance(edge, MapAxisEdge) + assert edge.name == "my-map-axis" + assert edge.output is yx_cs + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[2.0, 1.0], [4.0, 3.0]])) + + def test_csgen_output_is_permuted_input_axes(self, xy_cs, yx_cs): + model = ozm06trans.MapAxis(mapAxis=(1, 0)) + edge = parse_map_axis(model, input=xy_cs, out=CsGen(base_name="gen")) + assert edge.output.axes_names == yx_cs.axes_names + assert edge.output.virtual is True + + +class TestParseAffine: + def test_parses_linear_translation_and_name(self, xy_cs): + model = ozm06trans.Affine(name="my-affine", affine=((2.0, 0.0, 1.0), (0.0, 3.0, 5.0))) + edge = parse_affine(model, input=xy_cs, output=xy_cs) + assert isinstance(edge, AffineEdge) + np.testing.assert_allclose(edge.linear, [[2.0, 0.0], [0.0, 3.0]]) + np.testing.assert_allclose(edge.translation, [1.0, 5.0]) + assert edge.name == "my-affine" + + def test_csgen_output_derived_from_matrix_row_count(self, xy_cs): + model = ozm06trans.Affine(affine=((2.0, 0.0, 1.0), (0.0, 3.0, 5.0))) + edge = parse_affine(model, input=xy_cs, output=CsGen(base_name="gen")) + assert edge.output.num_axes == 2 + assert edge.output.virtual is True + + +class TestParseRotation: + def test_parses_rotation_matrix_and_name(self, xy_cs): + model = ozm06trans.Rotation(name="my-rotation", rotation=((0.0, -1.0), (1.0, 0.0))) + edge = parse_rotation(model, input=xy_cs, out=xy_cs) + assert isinstance(edge, RotationEdge) + np.testing.assert_allclose(edge.rotation, [[0.0, -1.0], [1.0, 0.0]]) + assert edge.name == "my-rotation" + + def test_csgen_output_derived_from_matrix_row_count(self, xy_cs): + model = ozm06trans.Rotation(rotation=((0.0, -1.0), (1.0, 0.0))) + edge = parse_rotation(model, input=xy_cs, out=CsGen(base_name="gen")) + assert edge.output.num_axes == 2 + assert edge.output.virtual is True + + +class TestParseSequence: + def test_parses_name_and_composes_inner_transformations_in_order(self, xy_cs): + inner_translation = ozm06trans.Translation(translation=(1.0, 2.0)) + inner_scale = ozm06trans.Scale(scale=(2.0, 3.0)) + model = ozm06trans.Sequence(name="my-sequence", transformations=(inner_translation, inner_scale)) + edge = parse_sequence(model, input=xy_cs, output=xy_cs) + assert isinstance(edge, SequenceEdge) + assert edge.name == "my-sequence" + assert edge.input is xy_cs + assert edge.output is xy_cs + assert [type(t) for t in edge.transformations] == [TranslationEdge, ScaleEdge] + expected = (POINTS_2D + np.array([1.0, 2.0])) * np.array([2.0, 3.0]) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), expected) + + def test_single_inner_transformation_uses_given_output_directly(self, xy_cs): + inner_translation = ozm06trans.Translation(translation=(1.0, 2.0)) + model = ozm06trans.Sequence(transformations=(inner_translation,)) + edge = parse_sequence(model, input=xy_cs, output=xy_cs) + assert len(edge.transformations) == 1 + assert edge.transformations[0].output is xy_cs + + def test_multiple_inner_transformations_use_virtual_intermediate_coord_systems(self, xy_cs): + inner_translation = ozm06trans.Translation(translation=(1.0, 2.0)) + inner_scale = ozm06trans.Scale(scale=(2.0, 3.0)) + model = ozm06trans.Sequence(transformations=(inner_translation, inner_scale)) + edge = parse_sequence(model, input=xy_cs, output=xy_cs) + intermediate = edge.transformations[0].output + assert intermediate is not xy_cs + assert intermediate.virtual is True + assert intermediate.axes_names == xy_cs.axes_names + assert edge.transformations[1].output is xy_cs + + def test_csgen_output_mirrors_input(self, xy_cs): + inner_translation = ozm06trans.Translation(translation=(1.0, 2.0)) + inner_scale = ozm06trans.Scale(scale=(2.0, 3.0)) + model = ozm06trans.Sequence(transformations=(inner_translation, inner_scale)) + edge = parse_sequence(model, input=xy_cs, output=CsGen(base_name="gen")) + assert edge.output.axes_names == xy_cs.axes_names + assert edge.output.virtual is True + + +class TestParseByDimension: + def test_parses_name_and_splits_axes_between_inner_transformations(self, xy_cs): + scale_x = ozm06trans.Scale(scale=(2.0,)) + translate_y = ozm06trans.Translation(translation=(5.0,)) + model = ozm06trans.ByDimension( + name="my-by-dim", + transformations=( + ozm06trans.ByDimensionTransform(transformation=scale_x, inputAxes=(0,), outputAxes=(0,)), + ozm06trans.ByDimensionTransform(transformation=translate_y, inputAxes=(1,), outputAxes=(1,)), + ), + ) + edge = parse_by_dimension(model, input=xy_cs, output=xy_cs) + assert isinstance(edge, ByDimensionEdge) + assert edge.name == "my-by-dim" + assert [type(t) for t in edge.transformations] == [ScaleEdge, TranslationEdge] + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[2.0, 7.0], [6.0, 9.0]])) + + def test_csgen_output_derived_from_max_output_axis_index(self, xy_cs): + scale_x = ozm06trans.Scale(scale=(2.0,)) + translate_y = ozm06trans.Translation(translation=(5.0,)) + model = ozm06trans.ByDimension( + transformations=( + ozm06trans.ByDimensionTransform(transformation=scale_x, inputAxes=(0,), outputAxes=(0,)), + ozm06trans.ByDimensionTransform(transformation=translate_y, inputAxes=(1,), outputAxes=(1,)), + ), + ) + edge = parse_by_dimension(model, input=xy_cs, output=CsGen(base_name="gen")) + assert edge.output.num_axes == 2 + assert edge.output.virtual is True diff --git a/tests/transformations/graph/edge/test_edge_transformations.py b/tests/transformations/graph/edge/test_edge_transformations.py new file mode 100644 index 000000000..31884b091 --- /dev/null +++ b/tests/transformations/graph/edge/test_edge_transformations.py @@ -0,0 +1,414 @@ +# pyright: strict + +from __future__ import annotations + +import numpy as np +import pytest + +from spatialdata.transformations._graph import ( + AffineEdge, + AxisRedefinitionError, + ByDimensionEdge, + DeterminantDifferentFromOne, + EmptyTransformSequenceError, + IdentityEdge, + IncompatibleCoordSystemsError, + MapAxisEdge, + MissingAxisError, + NotOrthonormalError, + ProjectAxisEdge, + RotationEdge, + ScaleEdge, + SequenceEdge, + TranslationEdge, + UnexpectedShapeError, + UnmappedAxisError, +) + +POINTS_2D = np.array([[1.0, 2.0], [3.0, 4.0]]) +POINTS_3D = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) + + +class TestAffineEdge: + def test_constructor_rejects_wrong_linear_shape(self, xy_cs): + with pytest.raises(UnexpectedShapeError): + AffineEdge(linear=np.eye(3), input=xy_cs, output=xy_cs) + + def test_constructor_rejects_wrong_translation_shape(self, xy_cs): + with pytest.raises(UnexpectedShapeError): + AffineEdge(linear=np.eye(2), translation=np.zeros(3), input=xy_cs, output=xy_cs) + + def test_from_affine_matrix(self, xy_cs): + # fmt: off + affine_matrix = np.array([ + [2.0, 0.0, 1.0], + [0.0, 3.0, 5.0], + [0.0, 0.0, 1.0] + ]) + linear = [ + [2.0, 0.0], + [0.0, 3.0], + ] + translation = [ + 1.0, + 5.0 + ] + # fmt: on + edge = AffineEdge.from_affine_matrix(name=None, affine_matrix=affine_matrix, input=xy_cs, output=xy_cs) + np.testing.assert_equal(edge.linear, linear) + np.testing.assert_equal(edge.translation, translation) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[3.0, 11.0], [7.0, 17.0]])) + + def test_transform_points_scale_and_translate(self, xy_cs): + edge = AffineEdge( + linear=np.array([[2.0, 0.0], [0.0, 3.0]]), translation=np.array([1.0, 1.0]), input=xy_cs, output=xy_cs + ) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[3.0, 7.0], [7.0, 13.0]])) + + def test_transform_points_2d_to_3d(self, xy_cs, xyz_cs): + linear = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + edge = AffineEdge(linear=linear, input=xy_cs, output=xyz_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[1.0, 2.0, 3.0], [3.0, 4.0, 7.0]])) + + def test_transform_points_rejects_wrong_shape(self, xy_cs): + edge = AffineEdge(linear=np.eye(2), input=xy_cs, output=xy_cs) + with pytest.raises(UnexpectedShapeError): + # expecting (n, 2), gets (2,) + edge.transform_points(np.asarray([1.0, 2.0])) + with pytest.raises(UnexpectedShapeError): + # expecting (n, 2), gets (n, 3) + edge.transform_points(POINTS_3D) + + def test_inverse_roundtrips(self, xy_cs): + edge = AffineEdge( + linear=np.array([[2.0, 0.0], [0.0, 4.0]]), translation=np.array([1.0, -1.0]), input=xy_cs, output=xy_cs + ) + inv = edge.inverse() + assert inv is not None + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + def test_inverse_returns_none_for_singular_matrix(self, xy_cs): + edge = AffineEdge(linear=np.zeros((2, 2)), input=xy_cs, output=xy_cs) + assert edge.inverse() is None + + +class TestIdentityEdge: + def test_constructor_rejects_mismatched_number_of_axes(self, xy_cs, xyz_cs): + with pytest.raises(IncompatibleCoordSystemsError): + IdentityEdge(None, input=xy_cs, output=xyz_cs) + + def test_transform_points_is_a_noop(self, xy_cs): + edge = IdentityEdge(None, input=xy_cs, output=xy_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), POINTS_2D) + + def test_to_affine_matches_transform_points(self, xy_cs): + edge = IdentityEdge(None, input=xy_cs, output=xy_cs) + affine = edge.to_affine() + np.testing.assert_allclose(affine.linear, np.eye(2)) + np.testing.assert_allclose(affine.transform_points(POINTS_2D), POINTS_2D) + + def test_inverse_is_still_the_identity(self, xy_cs): + edge = IdentityEdge(None, input=xy_cs, output=xy_cs) + inv = edge.inverse() + np.testing.assert_allclose(inv.transform_points(POINTS_2D), POINTS_2D) + + +class TestMapAxisEdge: + def test_transform_points_swaps_axes( + self, + xyz_cs, + abc_cs, + axis_x_space_with_unit, + axis_y_space_with_unit, + axis_z_space_with_unit, + axis_a_space_with_unit, + axis_b_space_with_unit, + axis_c_space_with_unit, + ): + edge = MapAxisEdge( + input=xyz_cs, + output=abc_cs, + input_to_output={ + axis_x_space_with_unit: axis_b_space_with_unit, + axis_y_space_with_unit: axis_c_space_with_unit, + axis_z_space_with_unit: axis_a_space_with_unit, + }, + ) + np.testing.assert_equal(edge.transform_points(POINTS_3D), np.array([[3.0, 1.0, 2.0], [6.0, 4.0, 5.0]])) + + def test_transform_points_permutation_of_three_axes( + self, + xyz_cs, + abc_cs, + axis_x_space_with_unit, + axis_y_space_with_unit, + axis_z_space_with_unit, + axis_a_space_with_unit, + axis_b_space_with_unit, + axis_c_space_with_unit, + ): + edge = MapAxisEdge( + input=xyz_cs, + output=abc_cs, + input_to_output={ + axis_x_space_with_unit: axis_c_space_with_unit, + axis_y_space_with_unit: axis_b_space_with_unit, + axis_z_space_with_unit: axis_a_space_with_unit, + }, + ) + np.testing.assert_allclose(edge.transform_points(POINTS_3D), np.array([[3.0, 2.0, 1.0], [6.0, 5.0, 4.0]])) + + def test_to_affine_matches_transform_points(self, xy_cs, abc_cs, axis_x_space_with_unit, axis_y_space_with_unit): + edge = MapAxisEdge( + input=xy_cs, + output=xy_cs, + input_to_output={ + axis_x_space_with_unit: axis_y_space_with_unit, + axis_y_space_with_unit: axis_x_space_with_unit, + }, + ) + affine = edge.to_affine() + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips( + self, + xyz_cs, + zyx_cs, + axis_x_space_with_unit, + axis_y_space_with_unit, + axis_z_space_with_unit, + ): + edge = MapAxisEdge( + input=xyz_cs, + output=zyx_cs, + input_to_output={ + axis_x_space_with_unit: axis_z_space_with_unit, + axis_y_space_with_unit: axis_y_space_with_unit, + axis_z_space_with_unit: axis_x_space_with_unit, + }, + ) + inv = edge.inverse() + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_3D)), POINTS_3D) + + +class TestProjectAxisEdge: + def test_constructor_rejects_dropped_axis_missing_from_input(self, xy_cs, axis_z_space_with_unit): + with pytest.raises(MissingAxisError): + ProjectAxisEdge(input=xy_cs, output=xy_cs, dropped_inputs={axis_z_space_with_unit}, created_outputs=set()) + + def test_constructor_rejects_created_axis_missing_from_output(self, xy_cs, axis_z_space_with_unit): + with pytest.raises(MissingAxisError): + ProjectAxisEdge(input=xy_cs, output=xy_cs, dropped_inputs=set(), created_outputs={axis_z_space_with_unit}) + + def test_general_coord_system_incompatibility(self, xy_cs, xyz_cs): + with pytest.raises(IncompatibleCoordSystemsError): + ProjectAxisEdge(input=xy_cs, output=xyz_cs, dropped_inputs=set(), created_outputs=set()) + with pytest.raises(IncompatibleCoordSystemsError): + ProjectAxisEdge(input=xyz_cs, output=xy_cs, dropped_inputs=set(), created_outputs=set()) + + def test_to_affine_for_identity_case(self, xy_cs, xyz_cs): + edge = ProjectAxisEdge(input=xy_cs, output=xy_cs, dropped_inputs=set(), created_outputs=set()) + affine = edge.to_affine() + np.testing.assert_allclose(affine.linear, np.eye(2)) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), POINTS_2D) + + def test_dropping_one_axis(self, xy_cs, xyz_cs, axis_z_space_with_unit): + edge = ProjectAxisEdge( + input=xyz_cs, output=xy_cs, dropped_inputs={axis_z_space_with_unit}, created_outputs=set() + ) + affine = edge.to_affine() + assert affine is not None + np.testing.assert_allclose(edge.transform_points(POINTS_3D), np.array([[1.0, 2.0], [4.0, 5.0]])) + + def test_creating_one_axis(self, xy_cs, xyz_cs, axis_z_space_with_unit): + edge = ProjectAxisEdge( + input=xy_cs, output=xyz_cs, dropped_inputs=set(), created_outputs={axis_z_space_with_unit} + ) + affine = edge.to_affine() + assert affine is not None + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[1.0, 2.0, 0.0], [3.0, 4.0, 0.0]])) + + def test_inverse_none_when_axes_are_dropped_or_created(self, xy_cs, xyz_cs, axis_z_space_with_unit): + dropped = ProjectAxisEdge( + input=xyz_cs, output=xy_cs, dropped_inputs={axis_z_space_with_unit}, created_outputs=set() + ) + assert dropped.inverse() is None + created = ProjectAxisEdge( + input=xy_cs, output=xyz_cs, dropped_inputs=set(), created_outputs={axis_z_space_with_unit} + ) + assert created.inverse() is None + + def test_inverse_roundtrip(self, xy_cs, yx_cs): + edge = ProjectAxisEdge(input=xy_cs, output=yx_cs, dropped_inputs=set(), created_outputs=set()) + inv = edge.inverse() + assert inv is not None + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + +class TestTranslationEdge: + def test_constructor_rejects_mismatched_number_of_axes(self, xy_cs, xyz_cs): + with pytest.raises(IncompatibleCoordSystemsError): + TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xyz_cs) + + def test_transform_points_adds_translation(self, xy_cs): + edge = TranslationEdge(translation=np.array([10.0, 20.0]), input=xy_cs, output=xy_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[11.0, 22.0], [13.0, 24.0]])) + + def test_to_affine_matches_transform_points(self, xy_cs): + edge = TranslationEdge(translation=np.array([10.0, 20.0]), input=xy_cs, output=xy_cs) + affine = edge.to_affine() + assert affine is not None + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips(self, xy_cs): + edge = TranslationEdge(translation=np.array([10.0, 20.0]), input=xy_cs, output=xy_cs) + inv = edge.inverse() + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + +class TestScaleEdge: + def test_constructor_rejects_wrong_scale_shape(self, xy_cs): + with pytest.raises(UnexpectedShapeError): + ScaleEdge(scale=np.array([1.0, 2.0, 3.0]), input=xy_cs, output=xy_cs) + + def test_constructor_rejects_mismatched_number_of_axes(self, xy_cs, xyz_cs): + with pytest.raises(IncompatibleCoordSystemsError): + ScaleEdge(scale=np.array([1.0, 2.0]), input=xy_cs, output=xyz_cs) + + def test_transform_points(self, xy_cs): + edge = ScaleEdge(scale=np.array([2.0, 4.0]), input=xy_cs, output=xy_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[2.0, 8.0], [6.0, 16.0]])) + + def test_to_affine_matches_transform_points(self, xy_cs): + edge = ScaleEdge(scale=np.array([2.0, 4.0]), input=xy_cs, output=xy_cs) + affine = edge.to_affine() + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips(self, xy_cs): + edge = ScaleEdge(scale=np.array([2.0, 4.0]), input=xy_cs, output=xy_cs) + inv = edge.inverse() + assert inv is not None + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + def test_inverse_returns_none_when_scale_is_zero(self, xy_cs): + edge = ScaleEdge(scale=np.array([0.0, 4.0]), input=xy_cs, output=xy_cs) + assert edge.inverse() is None + + +class TestRotationEdge: + def test_constructor_rejects_mismatched_number_of_axes(self, xy_cs, xyz_cs): + with pytest.raises(IncompatibleCoordSystemsError): + RotationEdge(linear_matrix=np.eye(2), input=xy_cs, output=xyz_cs) + + def test_constructor_rejects_wrong_shape(self, xy_cs): + with pytest.raises(UnexpectedShapeError): + RotationEdge(linear_matrix=np.eye(3), input=xy_cs, output=xy_cs) + + def test_constructor_rejects_matrix_with_det_diff_from_1(self, xy_cs): + with pytest.raises(DeterminantDifferentFromOne): + RotationEdge(linear_matrix=np.array([[1.0, 0.0], [0.0, -1.0]]), input=xy_cs, output=xy_cs) + + def test_constructor_rejects_non_orthonormal_matrix(self, xy_cs): + with pytest.raises(NotOrthonormalError): + RotationEdge(linear_matrix=np.array([[1.0, 0.0], [0.0, 2.0]]), input=xy_cs, output=xy_cs) + + def test_transform_points_rotates_90_degrees(self, xy_cs): + edge = RotationEdge(linear_matrix=np.array([[0.0, -1.0], [1.0, 0.0]]), input=xy_cs, output=xy_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[-2.0, 1.0], [-4.0, 3.0]])) + + def test_to_affine_matches_transform_points(self, xy_cs): + edge = RotationEdge(linear_matrix=np.array([[0.0, -1.0], [1.0, 0.0]]), input=xy_cs, output=xy_cs) + affine = edge.to_affine() + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips(self, xy_cs): + edge = RotationEdge(linear_matrix=np.array([[0.0, -1.0], [1.0, 0.0]]), input=xy_cs, output=xy_cs) + inv = edge.inverse() + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + +class TestSequenceEdge: + def test_constructor_rejects_empty_sequence(self): + with pytest.raises(EmptyTransformSequenceError): + SequenceEdge(transformations=[]) + + def test_constructor_rejects_incompatible_neighbors(self, xy_cs, xyz_cs): + first = TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xy_cs) + second = TranslationEdge(translation=np.array([1.0, 2.0, 3.0]), input=xyz_cs, output=xyz_cs) + with pytest.raises(IncompatibleCoordSystemsError): + SequenceEdge(transformations=[first, second]) + + def test_transform_points_composes_in_order(self, xy_cs): + translate = TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xy_cs) + scale = ScaleEdge(scale=np.array([3.0, 4.0]), input=xy_cs, output=xy_cs) + edge = SequenceEdge(transformations=[translate, scale]) + expected = (POINTS_2D + np.array([1.0, 2.0])) * np.array([3.0, 4.0]) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), expected) + + def test_to_affine_matches_transform_points(self, xy_cs): + translate = TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xy_cs) + scale = ScaleEdge(scale=np.array([3.0, 4.0]), input=xy_cs, output=xy_cs) + edge = SequenceEdge(transformations=[translate, scale]) + affine = edge.to_affine() + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips(self, xy_cs): + translate = TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xy_cs) + scale = ScaleEdge(scale=np.array([3.0, 4.0]), input=xy_cs, output=xy_cs) + edge = SequenceEdge(transformations=[translate, scale]) + inv = edge.inverse() + assert inv is not None + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) + + def test_inverse_returns_none_if_any_component_is_not_invertible(self, xy_cs): + scale = ScaleEdge(scale=np.array([0.0, 4.0]), input=xy_cs, output=xy_cs) + translate = TranslationEdge(translation=np.array([1.0, 2.0]), input=xy_cs, output=xy_cs) + edge = SequenceEdge(transformations=[scale, translate]) + assert edge.inverse() is None + + +class TestByDimensionEdge: + def test_constructor_rejects_input_axis_missing_from_overall_input(self, xyc_cs, xy_cs): + sub = IdentityEdge(None, input=xyc_cs, output=xyc_cs) + with pytest.raises(MissingAxisError): + ByDimensionEdge(transformations=[sub], input=xy_cs, output=xyc_cs) + + def test_constructor_rejects_output_axis_missing_from_overall_output(self, xy_cs, x_cs): + sub = IdentityEdge(None, input=xy_cs, output=xy_cs) + with pytest.raises(MissingAxisError): + ByDimensionEdge(transformations=[sub], input=xy_cs, output=x_cs) + + def test_constructor_rejects_output_axis_defined_more_than_once(self, x_cs): + first = IdentityEdge(None, input=x_cs, output=x_cs) + second = IdentityEdge(None, input=x_cs, output=x_cs) + with pytest.raises(AxisRedefinitionError): + ByDimensionEdge(transformations=[first, second], input=x_cs, output=x_cs) + + def test_constructor_rejects_unmapped_output_axis(self, x_cs, xy_cs, axis_y_space_with_unit): + sub = IdentityEdge(None, input=x_cs, output=x_cs) + try: + ByDimensionEdge(transformations=[sub], input=xy_cs, output=xy_cs) + raise AssertionError(f"Expected {UnmappedAxisError.__name__} to be raised") + except UnmappedAxisError as e: + assert e.axis_name == axis_y_space_with_unit.name + + def test_transform_points_applies_each_transformation_per_axis(self, x_cs, y_cs, xy_cs): + scale_x = ScaleEdge(scale=np.array([2.0]), input=x_cs, output=x_cs) + translate_y = TranslationEdge(translation=np.array([5.0]), input=y_cs, output=y_cs) + edge = ByDimensionEdge(transformations=[scale_x, translate_y], input=xy_cs, output=xy_cs) + np.testing.assert_allclose(edge.transform_points(POINTS_2D), np.array([[2.0, 7.0], [6.0, 9.0]])) + + def test_to_affine_matches_transform_points(self, x_cs, y_cs, xy_cs): + scale_x = ScaleEdge(scale=np.array([2.0]), input=x_cs, output=x_cs) + translate_y = TranslationEdge(translation=np.array([5.0]), input=y_cs, output=y_cs) + edge = ByDimensionEdge(transformations=[scale_x, translate_y], input=xy_cs, output=xy_cs) + affine = edge.to_affine() + np.testing.assert_allclose(affine.transform_points(POINTS_2D), edge.transform_points(POINTS_2D)) + + def test_inverse_roundtrips(self, x_cs, y_cs, xy_cs): + scale_x = ScaleEdge(scale=np.array([2.0]), input=x_cs, output=x_cs) + translate_y = TranslationEdge(translation=np.array([5.0]), input=y_cs, output=y_cs) + edge = ByDimensionEdge(transformations=[scale_x, translate_y], input=xy_cs, output=xy_cs) + inv = edge.inverse() + assert inv is not None + np.testing.assert_allclose(inv.transform_points(edge.transform_points(POINTS_2D)), POINTS_2D) diff --git a/tests/transformations/graph/graph/__init__.py b/tests/transformations/graph/graph/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/transformations/graph/graph/blobs_dataset_related.py b/tests/transformations/graph/graph/blobs_dataset_related.py new file mode 100644 index 000000000..976c55e91 --- /dev/null +++ b/tests/transformations/graph/graph/blobs_dataset_related.py @@ -0,0 +1,37 @@ +from __future__ import annotations + +from spatialdata.datasets import blobs +from spatialdata.transformations._graph import IdentityEdge, ProjectAxisEdge +from spatialdata.transformations._graph.utils import get_intrinsic_coordinate_system_name + + +def test_blobs_dataset_visualization_export(): + sdata = blobs() + # has elements of all types mapping to "global" cs with Identity transformation + graph = sdata.transformations._graph + + global_cs = graph.get_coordinate_system("global") + assert global_cs.axes_names == ("c", "y", "x") + + for element_name, _element in (sdata.images | sdata.labels | sdata.shapes | sdata.points).items(): + element_cs = graph.get_element_coordinate_system(element_name) + assert element_cs.name == get_intrinsic_coordinate_system_name(element_name) + + created_axes = set(global_cs.axes) - set(element_cs.axes) + + if created_axes: + intermediate_cs = graph.get_coordinate_system(f"global_{element_name}") + direct_edges_element_to_intermediate = graph.get_direct_transformations(element_cs, intermediate_cs) + assert len(direct_edges_element_to_intermediate) == 1 + assert isinstance(direct_edges_element_to_intermediate[0], IdentityEdge) + + direct_edges_intermediate_to_global = graph.get_direct_transformations(intermediate_cs, global_cs) + assert len(direct_edges_intermediate_to_global) == 1 + assert isinstance(direct_edges_intermediate_to_global[0], ProjectAxisEdge) + assert direct_edges_intermediate_to_global[0].created_outputs == created_axes + assert direct_edges_intermediate_to_global[0].dropped_inputs == set() + + else: + direct_edges_element_to_global = graph.get_direct_transformations(element_cs, global_cs) + assert len(direct_edges_element_to_global) == 1 + assert isinstance(direct_edges_element_to_global[0], IdentityEdge) diff --git a/tests/transformations/graph/graph/conftest.py b/tests/transformations/graph/graph/conftest.py new file mode 100644 index 000000000..592b73df7 --- /dev/null +++ b/tests/transformations/graph/graph/conftest.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 + +""" +Fixtures for transformation manager tests. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from spatialdata.transformations._graph import AffineEdge, Axis, CoordSystem, ScaleEdge, TranslationEdge + + +def get_coord_system(cs_name: str) -> CoordSystem: + return CoordSystem( + name=cs_name, + axes=(Axis(name="x", type="space", unit="micrometer"), Axis(name="y", type="space", unit="micrometer")), + virtual=False, + ) + + +@pytest.fixture +def one_point_graph() -> list[CoordSystem]: + """Fixture providing a single point graph with one coordinate system.""" + return [get_coord_system("cs1")] + + +@pytest.fixture +def fully_connected_two_point_graph() -> tuple[list[CoordSystem], list[TranslationEdge]]: + """Fixture providing a fully connected two-point graph with two coordinate systems and transformations.""" + coordinate_systems = [get_coord_system("cs1"), get_coord_system("cs2")] + cs1, cs2 = coordinate_systems + transformations = [ + TranslationEdge(translation=np.array([-1.0, -2.0]), input=cs1, output=cs2), + ] + return coordinate_systems, transformations + + +@pytest.fixture +def four_point_graph() -> tuple[list[CoordSystem], list[ScaleEdge | TranslationEdge]]: + """Fixture providing a four-point graph with four coordinate systems and transformations.""" + coordinate_systems = [ + get_coord_system("cs1"), + get_coord_system("cs2"), + get_coord_system("cs3"), + get_coord_system("cs4"), + ] + cs1, cs2, cs3, cs4 = coordinate_systems + transformations = [ + TranslationEdge(translation=np.array([1.0, 2.0]), input=cs1, output=cs2), # cs1 -> cs2 + TranslationEdge(translation=np.array([3.0, 4.0]), input=cs2, output=cs3), # cs2 -> cs3 + ScaleEdge(scale=np.array([2, 2]), input=cs3, output=cs4), # cs3 -> cs4 + TranslationEdge(translation=np.array([4.0, 6.0]), input=cs1, output=cs3), + # cs1 -> cs3 (consistent with cs1->cs2 and cs2->cs3) + ] + return coordinate_systems, transformations + + +@pytest.fixture +def five_point_graph() -> tuple[list[CoordSystem], list[ScaleEdge | TranslationEdge | AffineEdge]]: + """Fixture providing a five-point graph with five coordinate systems and five transformations.""" + + coordinate_systems = [ + get_coord_system("cs1"), + get_coord_system("cs2"), + get_coord_system("cs3"), + get_coord_system("cs4"), + get_coord_system("cs5"), + ] + cs1, cs2, cs3, cs4, cs5 = coordinate_systems + transformations = [ + TranslationEdge(translation=np.array([1.0, 2.0]), input=cs1, output=cs2), # cs1 -> cs2, + TranslationEdge(translation=np.array([1.0, 2.0]), input=cs4, output=cs3), # cs4 -> cs3 + TranslationEdge(translation=np.array([3.0, 4.0]), input=cs2, output=cs3), # cs2 -> cs3 + TranslationEdge(translation=np.array([3.0, 4.0]), input=cs1, output=cs4), # cs1 -> cs4 + TranslationEdge(translation=np.array([4.0, 6.0]), input=cs1, output=cs3), + # cs1 -> cs3 (consistent with cs1->cs2 and cs2->cs3) + ScaleEdge(scale=np.array([2, 2]), input=cs3, output=cs5), # cs3 -> cs5; + AffineEdge( + linear=np.identity(2), + translation=np.array([0, 0]), + input=cs3, + output=cs5, + ), # cs3 -> cs5 + ] + return coordinate_systems, transformations diff --git a/tests/transformations/graph/graph/test_graph.py b/tests/transformations/graph/graph/test_graph.py new file mode 100644 index 000000000..c5ca7e7ef --- /dev/null +++ b/tests/transformations/graph/graph/test_graph.py @@ -0,0 +1,667 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import pytest + +from spatialdata.literals import TRANSFORM_KEY +from spatialdata.transformations._graph import ( + BaseTransformationEdge, + CannotRemoveCoordinateSystemError, + CoordinateSystemAlreadyExistsError, + CoordinateSystemNotFoundError, + ElementAlreadyExistsError, + ElementNotRegisteredToAnyCoordinateSystemError, + InvalidPathError, + SequenceEdge, + TransformationGraph, + TransformationNotFoundError, + TransformationPathAmbiguousError, + TransformationPathAmbiguousMultipleEdgeExpectedError, + TransformationPathAmbiguousNoEdgeExpectedError, + TransformationPathNotFoundError, + TransformationPathNotSimple, +) + + +def test_initialization(): + """Test that TransformationManager initializes correctly.""" + + tm = TransformationGraph() + assert len(tm._graph.nodes()) == 0 + assert len(tm._graph.edges()) == 0 + assert tm._element_to_cs_mapping == {} + + +def test_add_coordinate_system(one_point_graph): + """Test adding a coordinate system.""" + + tm = TransformationGraph() + [cs1] = one_point_graph + + tm.add_coordinate_system(cs1) + assert cs1 in tm._graph.nodes() + + +def test_add_coordinate_system_duplicate(one_point_graph): + """Test adding an already existing coordinate system""" + tm = TransformationGraph() + [cs1] = one_point_graph + + tm.add_coordinate_system(cs1) + with pytest.raises(CoordinateSystemAlreadyExistsError, match=f"Coordinate system '{cs1.name}' already exists"): + tm.add_coordinate_system(cs1) + + +def test_remove_coordinate_system(one_point_graph): + """Test removing a coordinate system.""" + tm = TransformationGraph() + [cs1] = one_point_graph + + tm.add_coordinate_system(cs1) + tm.remove_coordinate_system(cs1) + + assert cs1 not in tm._graph.nodes() + + +def test_remove_coordinate_system_nonexistent(one_point_graph): + """Test that removing a non-existent coordinate system raises CoordinateSystemNotFoundError.""" + tm = TransformationGraph() + [cs1] = one_point_graph + + # Add the coordinate system first + tm.add_coordinate_system(cs1) + + # Remove it + tm.remove_coordinate_system(cs1) + + # Try to remove it again - should raise CoordinateSystemNotFoundError + with pytest.raises(CoordinateSystemNotFoundError): + tm.remove_coordinate_system(cs1) + + +def test_remove_coordinate_system_with_associated_transformations(fully_connected_two_point_graph): + """ + Test that removing a coordinate system with associated transformations raises CannotRemoveCoordinateSystemError. + """ + tm = TransformationGraph() + [cs1, cs2], [transformation] = fully_connected_two_point_graph + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + tm.add_transformation_edge(transformation) + + # Should raise CannotRemoveCoordinateSystem when trying to remove a coordinate system with transformations + with pytest.raises(CannotRemoveCoordinateSystemError, match="Cannot remove coordinate system"): + tm.remove_coordinate_system(cs1) + + # Should raise CannotRemoveCoordinateSystem when trying to remove a coordinate system with transformations + with pytest.raises(CannotRemoveCoordinateSystemError, match="Cannot remove coordinate system"): + tm.remove_coordinate_system(cs2) + + +def test_remove_coordinate_system_with_belonging_elements(one_point_graph): + """Test that removing a coordinate system with element associations raises CannotRemoveCoordinateSystemError.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + tm.add_element("image1", cs1) + + with pytest.raises(CannotRemoveCoordinateSystemError, match="Cannot remove coordinate system"): + tm.remove_coordinate_system(cs1) + + +def test_list_coordinate_systems(fully_connected_two_point_graph): + """Test listing all coordinate systems.""" + tm = TransformationGraph() + [cs1, cs2], _ = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + systems = tm.list_coordinate_systems() + assert len(systems) == 2 + assert cs1 in systems + assert cs2 in systems + + +def test_add_element(one_point_graph): + """Test adding an element with an existing coordinate system.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + tm.add_element("image1", cs1) + + mapping = tm._element_to_cs_mapping + assert "image1" in mapping + assert mapping["image1"] == cs1 + + +def test_add_element_duplicate(one_point_graph): + """Test that adding a duplicate element raises ElementAlreadyExistsError.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + element_name = "image1" + tm.add_element(element_name, cs1) + + # Try to add the same element again - should raise ElementAlreadyExistsError + with pytest.raises( + ElementAlreadyExistsError, match=f"Element '{element_name}' already exists in the transformation manager" + ): + tm.add_element(element_name, cs1) + + +def test_unset_element(one_point_graph): + """Test unsetting an element.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + tm.add_element("image1", cs1) + tm._unset_element("image1") + + assert "image1" not in tm._element_to_cs_mapping + + +def test_unset_element_nonexistent(one_point_graph): + """Test that unsetting a non-existent element raises ElementNotFoundError.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + element_name = "image1" + + # Try to unset non-existent element + with pytest.raises(ElementNotRegisteredToAnyCoordinateSystemError): + tm._unset_element(element_name) + + +def test_add_transformation(fully_connected_two_point_graph): + """Test adding a transformation between coordinate systems.""" + + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + edge_key = tm.add_transformation_edge(transform_edge) + + assert tm._graph.has_edge(cs1, cs2) + assert tm._graph[cs1][cs2][edge_key][TRANSFORM_KEY] is transform_edge + + +def test_add_transformation_existing_edge(fully_connected_two_point_graph): + """Test adding a transformation between coordinate systems that already have an edge.""" + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + tm.add_transformation_edge(transform_edge) + assert tm._graph.has_edge(cs1, cs2) + + # Add same transformation again between the same coordinate systems + tm.add_transformation_edge(transform_edge) + + assert tm._graph.has_edge(cs1, cs2) + assert len(tm._graph[cs1][cs2]) == 2 + + +def test_add_transformation_nonexistent_cs(fully_connected_two_point_graph): + """Test that adding a transformation with non-existent coordinate systems raises CoordinateSystemNotFoundError.""" + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + + expected_msg = "Coordinate system .* not found in the transformation manager." + with pytest.raises(CoordinateSystemNotFoundError, match=expected_msg): + tm.add_transformation_edge(transform_edge) + + # Add one coordinate system + tm.add_coordinate_system(cs1) + + with pytest.raises(CoordinateSystemNotFoundError, match=expected_msg): + tm.add_transformation_edge(transform_edge) + + +def test_get_direct_transformations(fully_connected_two_point_graph): + """Test getting an existing transformation.""" + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + tm.add_transformation_edge(transform_edge) + + retrieved = tm.get_direct_transformations(cs1, cs2) + assert len(retrieved) == 1 + assert next(iter(retrieved.values())) is transform_edge + + +def test_get_direct_transformations_nonexistent(fully_connected_two_point_graph): + """Test getting a non-existent transformation.""" + tm = TransformationGraph() + [cs1, cs2], _ = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + retrieved = tm.get_direct_transformations(cs1, cs2) + assert retrieved == {} + + +def test_remove_all_transformations_between_coordinate_systems(fully_connected_two_point_graph, mocker): + """Test removing all transformations between coordinate systems.""" + + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + tm.add_transformation_edge(transform_edge) + + spy = mocker.spy(tm._graph, "remove_edge") + tm.remove_all_transformations_between_coordinate_systems(cs1, cs2) + assert spy.call_count == 1 + assert not tm._graph.has_edge(cs1, cs2) + + +def test_remove_all_transformation_nonexistent(fully_connected_two_point_graph, mocker): + """Test that removing non-existent transformations between coordinate systems raises TransformationNotFoundError.""" + + tm = TransformationGraph() + [cs1, cs2], _ = fully_connected_two_point_graph + tm._graph.add_node(cs1) + tm._graph.add_node(cs2) + + spy = mocker.spy(tm._graph, "remove_edge") + tm.remove_all_transformations_between_coordinate_systems(cs1, cs2) + assert spy.call_count == 0 + + +def test_remove_specific_transformation_between_coordinate_systems(five_point_graph): + """Test removing specific transformation between coordinate systems.""" + + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [_cs1, _cs2, cs3, _cs4, cs5] = coord_systems + [_transform1, _transform2, _transform3, _transform4, _transform5, transform6, transform7] = transform_edges + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform6) + edge_key7 = tm.add_transformation_edge(transform7) + + tm.remove_specific_transformation(edge_key7, transform7) + assert not tm._graph.has_edge(cs3, cs5, key=edge_key7) + + +def test_remove_specific_transformation_between_coordinate_systems_non_existent(five_point_graph): + """ + Test that removing non-existent specific transformation between coordinate systems raises + TransformationNotFoundError. + """ + + tm = TransformationGraph() + coord_systems, transform_edges = five_point_graph + [_cs1, _cs2, cs3, _cs4, cs5] = coord_systems + [_transform1, _transform2, _transform3, _transform4, _transform5, transform6, transform7] = transform_edges + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs5) + + with pytest.raises( + TransformationNotFoundError, match=f"Transformation from '{cs3.name}' to '{cs5.name}' not found" + ): + tm.remove_specific_transformation("irrelevant", transform7) + + +def assert_sequences_equivalent( + result_sequences: list[SequenceEdge], expected_transform_lists: list[list[BaseTransformationEdge]] +): + + assert len(result_sequences) == len(expected_transform_lists) + assert all(isinstance(x, SequenceEdge) for x in result_sequences) + transformation_sequences_found = [x.transformations for x in result_sequences] + assert set(map(tuple, expected_transform_lists)) == set(map(tuple, transformation_sequences_found)) + + +def test_get_all_shortest_transformation_sequences_direct(four_point_graph): + """Test getting all shortest transformation sequences for a direct transformation path.""" + tm = TransformationGraph() + [cs1, cs2, cs3, cs4], [transform1, transform2, _transform3, transform4] = four_point_graph + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform4) + + sequences = tm.get_all_shortest_transformation_sequences(cs1, cs3) + assert_sequences_equivalent(sequences, [[transform4]]) + + +def test_get_all_shortest_transformation_sequences_indirect_one_path(four_point_graph): + """Test getting all shortest transformation sequences for one indirect transformation path.""" + tm = TransformationGraph() + [cs1, cs2, cs3, _cs4], [transform1, transform2, _transform3, _transform4] = four_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform2) + + sequences = tm.get_all_shortest_transformation_sequences(cs1, cs3) + assert_sequences_equivalent(sequences, [[transform1, transform2]]) + + +def test_get_all_shortest_transformation_sequences_indirect_two_paths(five_point_graph): + """Test getting all shortest transformation sequences for two indirect transformation paths.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, _cs5] = coord_systems + [transform1, transform2, transform3, transform4, _transform5, _transform6, _transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + + sequences = tm.get_all_shortest_transformation_sequences(cs1, cs3) + assert_sequences_equivalent(sequences, [[transform4, transform2], [transform1, transform3]]) + + +def test_get_all_shortest_transformation_sequences_no_path(four_point_graph): + """Test that getting all shortest transformation sequences with no path raises TransformationPathNotFoundError.""" + tm = TransformationGraph() + [cs1, cs2, cs3, cs4], [transform1, transform2, _transform3, _transform4] = four_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform2) + + expected_error_msg = f"No transformation path found from {cs1.name} to {cs4.name}" + with pytest.raises(TransformationPathNotFoundError, match=expected_error_msg): + tm.get_all_shortest_transformation_sequences(cs1, cs4) + + +def test_get_all_shortest_transformation_sequences_multiple_paths_multiple_edges_success(five_point_graph): + """Test getting all shortest transformation sequences, with multiple path and multiple edges between nodes.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, _transform5, transform6, transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) # cs1 -> cs2 + tm.add_transformation_edge(transform3) # cs2 -> cs3 + tm.add_transformation_edge(transform4) # cs1 -> cs4 + tm.add_transformation_edge(transform2) # cs4 -> cs3 + # tm.add_transformation(cs1, cs3, transform3) + key6 = tm.add_transformation_edge(transform6) # cs3 -> cs5 + tm.add_transformation_edge(transform7) # cs3 -> cs5 + + expected_intermediate_edges = [key6] + + sequences = tm.get_all_shortest_transformation_sequences( + cs1, cs5, expected_intermediate_edge_keys=expected_intermediate_edges + ) + assert_sequences_equivalent(sequences, [[transform1, transform3, transform6], [transform4, transform2, transform6]]) + + +def test_get_all_shortest_transformation_sequences_multiple_paths_multiple_edges_failure(five_point_graph): + """Test getting all shortest transformation sequences, with multiple path and multiple edges between nodes.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, _transform5, transform6, transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + # tm.add_transformation(cs1, cs3, transform3) + tm.add_transformation_edge(transform6) + tm.add_transformation_edge(transform7) + + expected_msg = f"Transformation Path ambiguous from {cs1.name} to {cs5.name}." + with pytest.raises(TransformationPathAmbiguousError, match=expected_msg): + tm.get_all_shortest_transformation_sequences(cs1, cs5) + + +def test_get_all_transformation_sequences(four_point_graph): + """Test getting all transformation sequences.""" + tm = TransformationGraph() + [cs1, cs2, cs3, cs4], [transform1, transform2, transform3, transform4] = four_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + + sequences = tm.get_all_transformation_sequences(cs1, cs4) + assert_sequences_equivalent(sequences, [[transform1, transform2, transform3], [transform4, transform3]]) + + +def test_get_all_transformation_sequences_multiple_paths(five_point_graph): + """Test getting all transformation sequences, with multiple path and multiple edges between nodes.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, transform5, _transform6, _transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform5) + + sequences = tm.get_all_transformation_sequences(cs1, cs3) + assert_sequences_equivalent(sequences, [[transform1, transform3], [transform4, transform2], [transform5]]) + + +def test_get_all_transformation_sequences_multiple_paths_multiple_edges_success(five_point_graph): + """Test getting all transformation sequences, with multiple paths and multiple edges between nodes.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, transform5, transform6, transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform5) + key6 = tm.add_transformation_edge(transform6) + tm.add_transformation_edge(transform7) + + expected_intermediate_edges = [key6] + sequences = tm.get_all_transformation_sequences( + cs1, cs5, expected_intermediate_edge_keys=expected_intermediate_edges + ) + assert_sequences_equivalent( + sequences, + [[transform1, transform3, transform6], [transform4, transform2, transform6], [transform5, transform6]], + ) + + +def test_get_all_transformation_sequences_multiple_paths_multiple_edges_no_edge_expected(five_point_graph): + """Test getting all transformation sequences with multiple paths and multiple edges, no edge expected.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, transform5, transform6, transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform5) + tm.add_transformation_edge(transform6) + tm.add_transformation_edge(transform7) + + expected_msg = f"Transformation Path ambiguous from {cs3.name} to {cs5.name}." + with pytest.raises(TransformationPathAmbiguousNoEdgeExpectedError, match=expected_msg): + tm.get_all_transformation_sequences(cs1, cs5) + + +def test_get_all_transformation_sequences_multiple_paths_multiple_edges_multiple_expected(five_point_graph): + """Test getting all transformation sequences with multiple paths and multiple edges, multiple edges expected.""" + tm = TransformationGraph() + + coord_systems, transform_edges = five_point_graph + [cs1, cs2, cs3, cs4, cs5] = coord_systems + [transform1, transform2, transform3, transform4, transform5, transform6, transform7] = transform_edges + + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + tm.add_coordinate_system(cs5) + + tm.add_transformation_edge(transform1) + tm.add_transformation_edge(transform3) + tm.add_transformation_edge(transform4) + tm.add_transformation_edge(transform2) + tm.add_transformation_edge(transform5) + key6 = tm.add_transformation_edge(transform6) + key7 = tm.add_transformation_edge(transform7) + + expected_intermediate_edges = [key6, key7] + + expected_msg = f"Transformation Path ambiguous from {cs3.name} to {cs5.name}." + with pytest.raises(TransformationPathAmbiguousMultipleEdgeExpectedError, match=expected_msg): + tm.get_all_transformation_sequences(cs1, cs5, expected_intermediate_edge_keys=expected_intermediate_edges) + + +def test_get_transformation_sequences_from_simple_paths_after_disambiguation_with_non_simple_paths( + fully_connected_two_point_graph, +): + """Test _get_transformation_sequences_from_simple_paths_after_disambiguation with non-simple paths.""" + tm = TransformationGraph() + [cs1, cs2], _ = fully_connected_two_point_graph + + tm.add_coordinate_system(cs1) + + paths = [[cs1, cs2, cs1, cs2]] + + with pytest.raises(TransformationPathNotSimple): + tm._get_transformation_sequences_from_simple_paths_after_disambiguation(paths, ()) + + +def test_get_transformation_sequences_from_simple_paths_after_disambiguation_with_single_node_path( + fully_connected_two_point_graph, +): + """Test _get_transformation_sequences_from_simple_paths_after_disambiguation with a path with only one node.""" + tm = TransformationGraph() + [cs1, cs2], [transformation] = fully_connected_two_point_graph + + tm.add_coordinate_system(cs1) + + paths = [[cs1]] + with pytest.raises(InvalidPathError): + tm._get_transformation_sequences_from_simple_paths_after_disambiguation(paths, ()) + + +def test_get_all_transformation_sequences_no_path(four_point_graph): + """Test getting all transformation sequences.""" + tm = TransformationGraph() + [cs1, cs2, cs3, cs4], [transform1, transform2, transform3, transform4] = four_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + tm.add_coordinate_system(cs3) + tm.add_coordinate_system(cs4) + + expected_error_msg = f"No transformation path found from {cs1.name} to {cs4.name}" + with pytest.raises(TransformationPathNotFoundError, match=expected_error_msg): + tm.get_all_shortest_transformation_sequences(cs1, cs4) + + +def test_get_element_coordinate_system_success(one_point_graph): + """Test successfully getting an element's coordinate system.""" + tm = TransformationGraph() + [cs1] = one_point_graph + tm.add_coordinate_system(cs1) + tm.add_element("image1", cs1) + + # This should successfully return the coordinate system + result = tm.get_element_coordinate_system("image1") + assert result == cs1 + + +def test_repr_method_comprehensive(fully_connected_two_point_graph): + """Test TransformationManager.__repr__ method comprehensively.""" + + # Test empty TransformationManager + tm_empty = TransformationGraph() + repr_empty = repr(tm_empty) + assert "TransformationManager" in repr_empty + assert "coordinate_systems=[]" in repr_empty + assert "coordinate_transforms=[]" in repr_empty + assert "elements=[]" in repr_empty + + # Test TransformationManager with data + tm = TransformationGraph() + [cs1, cs2], [transform_edge] = fully_connected_two_point_graph + tm.add_coordinate_system(cs1) + tm.add_coordinate_system(cs2) + + tm.add_transformation_edge(transform_edge) + tm.add_element("image1", cs1) + tm.add_element("image2", cs2) + + repr_str = repr(tm) + assert "TransformationManager" in repr_str + assert "cs1" in repr_str + assert "cs2" in repr_str + assert repr_str.find(repr(transform_edge)) + assert "image1" in repr_str + assert "image2" in repr_str diff --git a/tests/transformations/graph/utils/__init__.py b/tests/transformations/graph/utils/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/transformations/graph/utils/test_get_union_of_axes.py b/tests/transformations/graph/utils/test_get_union_of_axes.py new file mode 100644 index 000000000..b23d48e53 --- /dev/null +++ b/tests/transformations/graph/utils/test_get_union_of_axes.py @@ -0,0 +1,170 @@ +""" +Test suite for ``get_union_of_axes`` that uses the reusable ``Axis`` +fixtures defined in the top‑level ``conftest.py``. +Each logical behavior is expressed as a separate test function; when the +behavior is identical we use ``@pytest.mark.parametrize`` to drive the +different input combinations. +""" + +from __future__ import annotations + +import pytest + +from spatialdata.literals import NgffAxisType +from spatialdata.transformations._graph import Axis, get_union_of_axes + + +# ---------------------------------------------------------------------- +# 1. Reverse‑alphabetical ordering (basic union & mixed channel/space) +# ---------------------------------------------------------------------- +@pytest.mark.parametrize( + "axes_fixture_names,expected_names", + [ + # three pure spatial axes + (("axis_x_space", "axis_y_space", "axis_z_space"), ("z", "y", "x")), + # channel + full 3‑D spatial axes + (("axis_c_channel", "axis_x_space", "axis_y_space", "axis_z_space"), ("c", "z", "y", "x")), + # mixed groups: xy + full 3‑D (duplicates should be removed) + (("axes_xy", "axes_zyx"), ("z", "y", "x")), + # mixed groups: cyx + full 3‑D (duplicates removed, c kept) + (("axes_cyx", "axes_zyx"), ("c", "z", "y", "x")), + ], + ids=[ + "basic_spatial", + "channel_plus_spatial", + "xy_plus_zyx", + "cyx_plus_zyx", + ], +) +def test_ordering(axes_fixture_names, expected_names, request): + """ + Verify ordering of the resulting axis names. + The fixture names are resolved to the actual ``Axis`` objects. + """ + # Resolve each name (whether a single Axis fixture or a composite tuple fixture) + resolved = [] + for name in axes_fixture_names: + obj = request.getfixturevalue(name) + # Composite fixtures return a tuple; flatten them + if isinstance(obj, tuple): + resolved.extend(obj) + else: + resolved.append(obj) + + result = get_union_of_axes(tuple(resolved)) + assert tuple(ax.name for ax in result) == expected_names + + +# ---------------------------------------------------------------------- +# 2. Duplicate removal when name and type are identical +# ---------------------------------------------------------------------- +def test_duplicate_removal(axis_x_space, axis_y_space): + """ + Two identical ``x`` axes should be deduplicated, leaving ``y`` untouched. + """ + result = get_union_of_axes((axis_x_space, axis_x_space, axis_y_space)) + assert len(result) == 2 + assert result[0].name == "y" + assert result[1].name == "x" + + +# ---------------------------------------------------------------------- +# 3. Same name, different type – both must be kept +# ---------------------------------------------------------------------- +def test_same_name_different_type(axis_x_space): + """ + ``x`` as a spatial axis and ``x`` as a channel axis are distinct. + """ + axis_x_channel = Axis(name="x", type=NgffAxisType.CHANNEL) # ad‑hoc, no fixture needed + result = get_union_of_axes((axis_x_space, axis_x_channel)) + assert len(result) == 2 + # Reverse alphabetical: spatial ``x`` comes before channel ``x`` + assert result == (axis_x_channel, axis_x_space) + + +# ---------------------------------------------------------------------- +# 4. Place axis with unit before an equivalent axis without unit +# ---------------------------------------------------------------------- +def test_prefer_non_none_unit(axis_x_space, axis_x_space_with_unit): + """ + Between an axis with ``unit=None`` and one with a concrete unit, + the latter must survive duplicate removal. + """ + result = get_union_of_axes((axis_x_space, axis_x_space_with_unit)) + assert len(result) == 2 + assert result[0].unit == "micrometer" + + +# ---------------------------------------------------------------------- +# 5. Place axis with long_name before an equivalent axis without long_name +# ---------------------------------------------------------------------- +def test_prefer_non_none_long_name( + axis_x_space, + axis_x_space_with_long_name, +): + """ + When only one duplicate defines ``long_name``, that value must be kept. + """ + result = get_union_of_axes((axis_x_space, axis_x_space_with_long_name)) + assert len(result) == 2 + assert result[0].long_name == "X axis" + + +# ---------------------------------------------------------------------- +# 6. ``unit`` outranks ``long_name`` when they are on different duplicates +# ---------------------------------------------------------------------- +def test_unit_beats_long_name( + axis_x_space_with_unit, + axis_x_space_with_long_name, +): + """ + The sorting key places ``unit_not_none`` before ``long_name_not_none``, + therefore the axis that defines a unit must be chosen even if another + duplicate defines a long name. + """ + result = get_union_of_axes((axis_x_space_with_long_name, axis_x_space_with_unit)) + assert len(result) == 2 + assert result[0].unit == "micrometer" + assert result[0].long_name is None + assert result[1].unit is None + assert result[1].long_name == "X axis" + + +# ---------------------------------------------------------------------- +# 7. Full hierarchy – axis with both optional fields is placed first +# ---------------------------------------------------------------------- +def test_full_hierarchy( + axis_x_space_full, + axis_x_space_with_unit, + axis_x_space, +): + """ + Three duplicate ``x`` axes: + * full (unit + long_name) + * unit only + * nothing + The most complete axis must be kept. + """ + result = get_union_of_axes((axis_x_space, axis_x_space_with_unit, axis_x_space_full)) + assert len(result) == 3 + assert result[0].unit == "micrometer" + assert result[0].long_name == "X axis" + assert result[1].unit == "micrometer" + assert result[1].long_name is None + assert result[2].unit is None + assert result[2].long_name is None + + +# ---------------------------------------------------------------------- +# 8. Empty input & single‑axis sanity checks +# ---------------------------------------------------------------------- +def test_empty_input(): + """Empty tuple must return an empty tuple.""" + assert get_union_of_axes(()) == () + + +def test_single_axis(axis_x_space): + """A single axis should be returned unchanged.""" + result = get_union_of_axes((axis_x_space,)) + assert len(result) == 1 + assert result[0] is axis_x_space diff --git a/tests/transformations/ngff/test_ngff_transformations.py b/tests/transformations/ngff/test_ngff_transformations.py index 7051f3646..ef3803e15 100644 --- a/tests/transformations/ngff/test_ngff_transformations.py +++ b/tests/transformations/ngff/test_ngff_transformations.py @@ -8,7 +8,7 @@ import pytest from spatialdata._types import ArrayLike -from spatialdata.models import C, X, Y, Z +from spatialdata.literals import C, X, Y, Z from spatialdata.transformations.ngff._utils import get_default_coordinate_system from spatialdata.transformations.ngff.ngff_coordinate_system import NgffCoordinateSystem from spatialdata.transformations.ngff.ngff_transformations import ( diff --git a/tests/transformations/test_transformations.py b/tests/transformations/test_transformations.py index dd5b52dba..63fcf7ecf 100644 --- a/tests/transformations/test_transformations.py +++ b/tests/transformations/test_transformations.py @@ -10,8 +10,9 @@ from spatialdata import transform from spatialdata.datasets import blobs -from spatialdata.models import Image2DModel, PointsModel -from spatialdata.models._utils import DEFAULT_COORDINATE_SYSTEM, ValidAxis_t, get_channel_names +from spatialdata.literals import DEFAULT_COORDINATE_SYSTEM +from spatialdata.models import Image2DModel, PointsModel, get_channel_names +from spatialdata.transformations._graph.types import ValidAxis_t from spatialdata.transformations.ngff._utils import get_default_coordinate_system from spatialdata.transformations.ngff.ngff_coordinate_system import NgffCoordinateSystem from spatialdata.transformations.ngff.ngff_transformations import (