Merge pull request #34 from lucas-diedrich/dvpio-writer-1

[Feature] DVPIO Writer (pyLMD shapes)

.github/workflows/test.yaml

@@ -30,9 +30,9 @@ jobs:
           - os: ubuntu-latest
             python: "3.10"
           - os: ubuntu-latest
-            python: "3.12"
+            python: "3.11"
           - os: ubuntu-latest
-            python: "3.12"
+            python: "3.11"
             pip-flags: "--pre"
             name: PRE-RELEASE DEPENDENCIES
 

pyproject.toml

@@ -27,7 +27,7 @@ dependencies = [
   "napari-spatialdata",
   "openslide-bin",
   "openslide-python",
-  "py-lmd",
+  "py-lmd @ git+https://github.com/MannLabs/py-lmd.git",
   "pylibczirw",
   "spatialdata",
   "spatialdata-plot",
@@ -77,6 +77,9 @@ scripts.clean = "git clean -fdX -- {args:docs}"
 [tool.hatch.envs.hatch-test]
 features = [ "test" ]
 
+[tool.hatch.metadata]
+allow-direct-references = true
+
 [tool.ruff]
 line-length = 120
 src = [ "src" ]

src/dvpio/read/shapes/__init__.py

@@ -1,4 +1,3 @@
-from .geometry import transform_shapes
-from .lmd_reader import read_lmd
+from .lmd_reader import read_lmd, transform_shapes
 
 __all__ = ["read_lmd", "transform_shapes"]

src/dvpio/read/shapes/geometry.py

@@ -1,20 +1,13 @@
-import geopandas as gpd
 import numpy as np
-import shapely
 from numpy.typing import NDArray
-from spatialdata.models import ShapesModel
-
-
-def _polygon_to_array(polygon):
-    return np.array(polygon.exterior.coords)
 
 
 def compute_affine_transformation(
     query_points: NDArray[np.float64], reference_points: NDArray[np.float64], precision: int | None = None
 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Computes the affine transformation mapping query_points to reference_points.
 
-    .. math:
+    .. math::
         Aq = r
 
     Parameters
@@ -40,23 +33,22 @@ def compute_affine_transformation(
         raise ValueError("At least three points are required to compute the transformation.")
 
     query_points = np.concatenate([query_points, np.ones(shape=(query_points.shape[0], 1))], axis=1)
+    reference_points = np.concatenate([reference_points, np.ones(shape=(reference_points.shape[0], 1))], axis=1)
     affine_matrix, _, _, _ = np.linalg.lstsq(query_points, reference_points, rcond=None)
 
     if precision is not None:
         affine_matrix = np.around(affine_matrix, precision)
 
-    rotation, translation = affine_matrix[:2, :], affine_matrix[2, :].reshape(1, -1)
-    return rotation, translation
+    return affine_matrix
 
 
 def apply_affine_transformation(
     shape: NDArray[np.float64],
-    rotation: NDArray[np.float64],
-    translation: NDArray[np.float64] | None = None,
+    affine_transformation: NDArray[np.float64],
 ) -> NDArray[np.float64]:
     """Transform shapes between coordinate systems
 
-    Applies rotation, translation, and switch of coordinate systems to a shape,
+    Applies affine transformation to a shape,
     in this order.
 
     Parameters
@@ -71,64 +63,11 @@ def apply_affine_transformation(
     Returns
     -------
     NDArray[np.float64]
-        Shape after affine transformation.
+        Shape (N, 2) after affine transformation.
     """
-    if translation is None:
-        # Identity translation
-        translation = np.zeros(shape=(1, 2))
-
-    return shape @ rotation + translation
-
-
-def transform_shapes(
-    shapes: gpd.GeoDataFrame | ShapesModel,
-    calibration_points_target: gpd.GeoDataFrame | ShapesModel,
-    calibration_points_source: gpd.GeoDataFrame | ShapesModel,
-    precision: int = 3,
-) -> ShapesModel:
-    """Apply coordinate transformation to shapes based on calibration points from a target and a source
-
-    Computes transformation between source and target coordinates.
-
-    Parameters
-    ----------
-    shapes
-        Shapes in source coordinate system (usually LMD coordinates)
-    calibration_points_target
-        3 Calibration points in target coordinate system (usually image/pixel coordinates)
-        Expects :class:`geopandas.GeoDataFrame` with calibration points in `geometry` column
-    calibration_points_source
-        3 Calibration points, matched to `calibration_points_target` in source coordinate system (usually LMD coordinates)
-        Expects :class:`geopandas.GeoDataFrame` with calibration points in `geometry` column
-    precision
-        Precision of affine transformation
-
-    Returns
-    -------
-    ShapesModel
-        Transformed shapes in target coordinate system
-    """
-    # Convert to numpy arrays
-    calibration_points_source = np.array(
-        calibration_points_source["geometry"].apply(lambda point: [point.x, point.y]).tolist()
-    )
-    calibration_points_target = np.array(
-        calibration_points_target["geometry"].apply(lambda point: [point.x, point.y]).tolist()
-    )
-
-    # Compute rotation (2x2) and translation (2x1) matrices
-    rotation, translation = compute_affine_transformation(
-        calibration_points_source, calibration_points_target, precision=precision
-    )
-    # Transform shapes
-    # Iterate through shapes and apply affine transformation
-    transformed_shapes = shapes["geometry"].apply(lambda shape: _polygon_to_array(shape))
-    transformed_shapes = transformed_shapes.apply(
-        lambda shape: apply_affine_transformation(shape, rotation=rotation, translation=translation)
-    )
-    transformed_shapes = transformed_shapes.apply(lambda shape: shapely.Polygon(shape))
-
-    # Reassign as DataFrame and parse with spatialdata
-    transformed_shapes = shapes.assign(geometry=transformed_shapes)
-
-    return ShapesModel.parse(transformed_shapes)
+    # Extend shape with ones
+    shape_mod = np.hstack([shape, np.ones(shape=(shape.shape[0], 1))])
+    # Apply affine transformation
+    shape_transformed = shape_mod @ affine_transformation
+    # Reuturn shape without padded ones
+    return shape_transformed[:, :-1]

src/dvpio/read/shapes/lmd_reader.py

@@ -1,10 +1,88 @@
-import geopandas as gpd
 import lmd.lib as pylmd
 import numpy as np
 import shapely
 from spatialdata.models import PointsModel, ShapesModel
+from spatialdata.transformations import Affine, set_transformation
 
-from .geometry import transform_shapes
+from .geometry import apply_affine_transformation, compute_affine_transformation
+
+
+def transform_shapes(
+    shapes: ShapesModel,
+    calibration_points_target: PointsModel,
+    calibration_points_source: PointsModel,
+    precision: int = 3,
+) -> ShapesModel:
+    """Apply coordinate transformation to shapes based on calibration points from a target and a source
+
+    Computes transformation between source and target coordinates.
+
+    Parameters
+    ----------
+    shapes
+        Shapes in source coordinate system (usually LMD coordinates)
+    calibration_points_target
+        3 Calibration points in target coordinate system (usually image/pixel coordinates)
+        Expects :class:`spatialdata.models.PointsModel` with calibration points in `x`/`y` column
+    calibration_points_source
+        3 Calibration points, matched to `calibration_points_target` in source coordinate system (usually LMD coordinates)
+        Expects :class:`spatialdata.models.PointsModel` with calibration points in `x`/`y` column
+    precision
+        Precision of affine transformation
+
+    Returns
+    -------
+    `class`::`spatialdata.models.ShapesModel`
+        Transformed shapes in target coordinate system
+
+    Object has special attributes
+
+    - `ShapesModel.attrs.transformation`
+        - global: (image coordinates)
+        - leica_micro_dissection: Leica coordinate system transformation
+
+    Raises
+    ------
+    AttributeError
+        Checks validity of shapes and calibration points data formats
+    """
+    ShapesModel.validate(shapes)
+    PointsModel.validate(calibration_points_source)
+    PointsModel.validate(calibration_points_target)
+
+    # Convert to numpy arrays
+    calibration_points_source = calibration_points_source[["x", "y"]].to_dask_array().compute()
+    calibration_points_target = calibration_points_target[["x", "y"]].to_dask_array().compute()
+
+    # Compute rotation (2x2) and translation (2x1) matrices
+    affine_transformation = compute_affine_transformation(
+        calibration_points_source, calibration_points_target, precision=precision
+    )
+
+    affine_transformation_inverse = np.around(np.linalg.inv(affine_transformation), precision)
+
+    # Transform shapes
+    # Iterate through shapes and apply affine transformation
+    transformed_shapes = shapes["geometry"].apply(
+        lambda shape: shapely.transform(
+            shape, transformation=lambda geom: apply_affine_transformation(geom, affine_transformation)
+        )
+    )
+
+    # Reassign as DataFrame and parse with spatialdata
+    transformed_shapes = ShapesModel.parse(shapes.assign(geometry=transformed_shapes))
+
+    # Set inverse transformation as transformation to leica coordinate system
+    set_transformation(
+        transformed_shapes,
+        transformation=Affine(affine_transformation_inverse.T, input_axes=("x", "y"), output_axes=("x", "y")),
+        to_coordinate_system="to_lmd",
+    )
+
+    # Store original calibration points
+    transformed_shapes.attrs["lmd_calibration_points"] = calibration_points_source
+
+    return transformed_shapes
 
 
 def read_lmd(path: str, calibration_points_image: PointsModel, switch_orientation: bool = False) -> ShapesModel:
@@ -26,31 +104,31 @@ def read_lmd(path: str, calibration_points_image: PointsModel, switch_orientatio
 
     Returns
     -------
-    ShapesModel
-        Transformed shapes in image coordinats
+    `class`::`spatialdata.models.ShapesModel`
+        Transformed shapes in image coordinates.
+
+        Object has special attributes
+
+            - `ShapesModel.attrs.transformation`
+                - global: (image coordinates)
+                - to_lmd: Transformation back to leica coordinate system
     """
+    PointsModel.validate(calibration_points_image)
+
     # Load LMD shapes with pyLMD
     lmd_shapes = pylmd.Collection()
     lmd_shapes.load(path)
+    shapes = lmd_shapes.to_geopandas("name", "well")
 
-    # Transform to geopandas
-    shapes = gpd.GeoDataFrame(geometry=[shapely.Polygon(shape.points) for shape in lmd_shapes.shapes])
-
-    calibration_points_lmd = gpd.GeoDataFrame(
-        data={"radius": np.ones(shape=len(lmd_shapes.calibration_points))},
-        geometry=[shapely.Point(point) for point in lmd_shapes.calibration_points],
-    )
-
-    calibration_points_image = gpd.GeoDataFrame(
-        data={"radius": np.ones(shape=len(calibration_points_image))},
-        geometry=[shapely.Point(row["x"], row["y"]) for _, row in calibration_points_image.iterrows()],
-    )
+    # Transform to spatialdata models
+    shapes = ShapesModel.parse(shapes)
+    calibration_points_lmd = PointsModel.parse(lmd_shapes.calibration_points)
 
     if len(calibration_points_lmd) < 3:
         raise ValueError(f"Require at least 3 calibration points, but only received {len(calibration_points_lmd)}")
     if len(calibration_points_lmd) != len(calibration_points_image):
         raise ValueError(
-            f"Number of calibration points in image ({len(calibration_points_image)})must be equal to number of calibration points in LMD file ({len(calibration_points_lmd)})"
+            f"Number of calibration points in image ({len(calibration_points_image)}) must be equal to number of calibration points in LMD file ({len(calibration_points_lmd)})"
         )
 
     transformed_shapes = transform_shapes(

src/dvpio/write/__init__.py

@@ -0,0 +1,3 @@
+from .lmd_writer import write_lmd
+
+__all__ = ["write_lmd"]