sopa.segmentation.shapes

sopa.segmentation.shapes.solve_conflicts(cells, threshold=0.5, patch_indices=None, return_indices=False)

Resolve segmentation conflicts (i.e. overlap) after running segmentation on patches

Parameters:

Name	Type	Description	Default
cells	list[Polygon]	List of cell polygons	required
threshold	float	When two cells are overlapping, we look at the area of intersection over the area of the smallest cell. If this value is higher than the threshold, the cells are merged	0.5
patch_indices	ndarray | None	Patch from which each cell belongs.	None
return_indices	bool	If True, returns also the cells indices. Merged cells have an index of -1.	False

Returns:

Type	Description
ndarray[Polygon] | tuple[ndarray[Polygon], ndarray]	Array of resolved cells polygons. If return_indices, it also returns an array of cell indices.

Source code in sopa/segmentation/shapes.py

def solve_conflicts(
    cells: list[Polygon],
    threshold: float = 0.5,
    patch_indices: np.ndarray | None = None,
    return_indices: bool = False,
) -> np.ndarray[Polygon] | tuple[np.ndarray[Polygon], np.ndarray]:
    """Resolve segmentation conflicts (i.e. overlap) after running segmentation on patches

    Args:
        cells: List of cell polygons
        threshold: When two cells are overlapping, we look at the area of intersection over the area of the smallest cell. If this value is higher than the `threshold`, the cells are merged
        patch_indices: Patch from which each cell belongs.
        return_indices: If `True`, returns also the cells indices. Merged cells have an index of -1.

    Returns:
        Array of resolved cells polygons. If `return_indices`, it also returns an array of cell indices.
    """
    cells = list(cells)
    n_cells = len(cells)
    resolved_indices = np.arange(n_cells)

    assert n_cells > 0, "No cells was segmented, cannot continue"

    tree = shapely.STRtree(cells)
    conflicts = tree.query(cells, predicate="intersects")

    if patch_indices is not None:
        conflicts = conflicts[:, patch_indices[conflicts[0]] != patch_indices[conflicts[1]]].T
    else:
        conflicts = conflicts[:, conflicts[0] != conflicts[1]].T

    for i1, i2 in tqdm(conflicts, desc="Resolving conflicts"):
        resolved_i1: int = resolved_indices[i1]
        resolved_i2: int = resolved_indices[i2]
        cell1, cell2 = cells[resolved_i1], cells[resolved_i2]

        intersection = cell1.intersection(cell2).area
        if intersection >= threshold * min(cell1.area, cell2.area):
            cell = _ensure_polygon(cell1.union(cell2))

            resolved_indices[np.isin(resolved_indices, [resolved_i1, resolved_i2])] = len(cells)
            cells.append(cell)

    unique_indices = np.unique(resolved_indices)
    unique_cells = np.array(cells)[unique_indices]

    if return_indices:
        return unique_cells, np.where(unique_indices < n_cells, unique_indices, -1)

    return unique_cells

sopa.segmentation.shapes.geometrize(mask, tolerance=None, smooth_radius_ratio=0.1)

Convert a cells mask to multiple shapely geometries. Inspired from https://github.com/Vizgen/vizgen-postprocessing

Parameters:

Name	Type	Description	Default
mask	ndarray	A cell mask. Non-null values correspond to cell ids	required
tolerance	float | None	Tolerance parameter used by shapely during simplification. By default, define the tolerance automatically.	None

Returns:

Type	Description
list[Polygon]	List of shapely polygons representing each cell ID of the mask

Source code in sopa/segmentation/shapes.py

def geometrize(mask: np.ndarray, tolerance: float | None = None, smooth_radius_ratio: float = 0.1) -> list[Polygon]:
    """Convert a cells mask to multiple `shapely` geometries. Inspired from https://github.com/Vizgen/vizgen-postprocessing

    Args:
        mask: A cell mask. Non-null values correspond to cell ids
        tolerance: Tolerance parameter used by `shapely` during simplification. By default, define the tolerance automatically.

    Returns:
        List of `shapely` polygons representing each cell ID of the mask
    """
    max_cells = mask.max()

    if max_cells == 0:
        log.warn("No cell was returned by the segmentation")
        return []

    cells = [_contours((mask == cell_id).astype("uint8")) for cell_id in range(1, max_cells + 1)]

    mean_radius = np.sqrt(np.array([cell.area for cell in cells]) / np.pi).mean()
    smooth_radius = mean_radius * smooth_radius_ratio

    if tolerance is None:
        tolerance = _default_tolerance(mean_radius)

    cells = [_smoothen_cell(cell, smooth_radius, tolerance) for cell in cells]
    cells = [cell for cell in cells if cell is not None]

    log.info(
        f"Percentage of non-geometrized cells: {(max_cells - len(cells)) / max_cells:.2%} (usually due to segmentation artefacts)"
    )

    return cells

sopa.segmentation.shapes.rasterize(cell, shape, xy_min=[0, 0])

Transform a cell polygon into a numpy array with value 1 where the polygon touches a pixel, else 0.

Parameters:

Name	Type	Description	Default
cell	Polygon | MultiPolygon	Cell polygon to rasterize.	required
shape	tuple[int, int]	Image shape as a tuple (y, x).	required
xy_min	tuple[int, int]	Tuple containing the origin of the image [x0, y0].	[0, 0]

Returns:

Type	Description
ndarray	The mask array.

Source code in sopa/segmentation/shapes.py

def rasterize(cell: Polygon | MultiPolygon, shape: tuple[int, int], xy_min: tuple[int, int] = [0, 0]) -> np.ndarray:
    """Transform a cell polygon into a numpy array with value 1 where the polygon touches a pixel, else 0.

    Args:
        cell: Cell polygon to rasterize.
        shape: Image shape as a tuple (y, x).
        xy_min: Tuple containing the origin of the image [x0, y0].

    Returns:
        The mask array.
    """
    import cv2

    xmin, ymin, xmax, ymax = [xy_min[0], xy_min[1], xy_min[0] + shape[1], xy_min[1] + shape[0]]

    cell_translated = shapely.affinity.translate(cell, -xmin, -ymin)
    geoms = cell_translated.geoms if isinstance(cell_translated, MultiPolygon) else [cell_translated]

    rasterized_image = np.zeros((ymax - ymin, xmax - xmin), dtype=np.int8)

    for geom in geoms:
        coords = np.array(geom.exterior.coords)[None, :].astype(np.int32)
        cv2.fillPoly(rasterized_image, coords, color=1)

    return rasterized_image