#!/usr/bin/env python
# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Copyright (c) 2022 Univ. Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAe, Inria
# All rights reserved.
# This file is part of the TimageTK library, and is released under the "GPLv3"
# license. Please see the LICENSE.md file that should have been included as
# part of this package.
# ------------------------------------------------------------------------------
import numpy as np
from scipy.cluster.vq import vq
from timagetk.components.labelled_image import LabelledImage
from timagetk.synthetic_data.util import point_position_optimization
[docs]
def example_layered_sphere_labelled_image(extent=60., voxelsize=(0.6, 0.6, 0.6), n_points=12, n_layers=1,
return_points=False, **kwargs):
"""Generate a synthetic 3D labelled image representing a spherical tissue.
The function creates a 3D labelled image over a cubic volume of a given extent.
It consists of successive spherical cell layers of equal height around a central spherical cell.
The number of cells in each layer is set so that cells in every layer have similar volumes.
Parameters
----------
extent : float, optional
The extent of the output image (in µm).
Defaults to ``60.``.
voxelsize : list(float), optional
Length-3 list of image Z,Y,X voxelsize (in µm).
Defaults to ``(0.6, 0.6, 0.6)``.
n_points : int, optional
The number of cells on innermost spherical layer.
Defaults to ``12``.
n_layers : int, optional
The number of layer around the central cell.
Defaults to ``1``.
return_points : bool, optional
Whether to return the points representing the cell centers.
Defaults to ``False``.
Returns
-------
timagetk.LabelledImage
The generated labelled cell image
dict, optional
The 3D X,Y,Z point positions used to generate the cells
Examples
--------
>>> from timagetk.synthetic_data.labelled_image import example_layered_sphere_labelled_image
>>> seg_img = example_layered_sphere_labelled_image()
>>> print(seg_img)
LabelledImage object with following metadata:
- shape: (100, 100, 100)
- ndim: 3
- dtype: uint16
- axes_order: ZYX
- voxelsize: [0.6, 0.6, 0.6]
- unit: 1e-06
- origin: [0, 0, 0]
- extent: [59.4, 59.4, 59.4]
- acquisition_date: None
- not_a_label: 0
>>> from timagetk.visu.stack import orthogonal_view
>>> v = orthogonal_view(seg_img, cmap='glasbey', val_range='auto')
"""
rng = np.random.default_rng(2021)
size = np.round(extent / np.array(voxelsize)).astype(int)
center = np.array([extent, extent, extent], float) / 2.
coords = np.transpose(np.mgrid[0:size[0], 0:size[1], 0:size[2]],
(1, 2, 3, 0)).reshape((np.prod(size), 3)).astype(int)
coords_points = np.array(coords, float) * np.array(voxelsize)
coords_distances = np.linalg.norm(coords_points - center, axis=1)
points = {2: center}
layer_img = {}
for layer in range(n_layers):
layer_n_points = n_points * np.power(layer + 1, 2)
radius = (layer + 1) * extent / float(2 * n_layers + 1)
layer_points = {}
max_point_id = np.max(list(points.keys()))
for p in range(layer_n_points):
theta = rng.random() * 2. * np.pi
phi = rng.random() * np.pi - np.pi / 2.
sphere_point = np.array([np.cos(theta) * np.cos(phi), np.sin(theta) * np.cos(phi), np.sin(phi)])
layer_points[max_point_id + 1 + p] = center + radius * sphere_point
point_target_area = 4. * np.pi * np.power(radius, 2.) / float(layer_n_points)
point_target_distance = np.power(point_target_area / np.pi, 0.5)
sigma_deformation = (extent / 100.) * (20. / layer_n_points)
omega_attraction = 0.1 * extent / 100.
omega_repulsion = 100.0 * np.power(extent / 100., 2)
layer_points = point_position_optimization(layer_points,
omega_attraction=omega_attraction,
omega_repulsion=omega_repulsion,
target_distance=point_target_distance,
sigma_deformation=sigma_deformation,
constraint_radius=True,
center=center,
radius=radius)
for p in layer_points.keys():
points[p] = layer_points[p]
layer_points_labels = np.array(list(layer_points.keys()))
layer_points_array = np.array(list(layer_points.values()))
labels = layer_points_labels[vq(coords_points, layer_points_array)[0]]
layer_img[layer + 1] = np.ones(size, np.uint16)
layer_img[layer + 1][tuple(np.transpose(coords))] = labels
img = np.ones(size, np.uint16)
for layer in range(n_layers):
layer_coords = coords[(coords_distances > (2 * layer + 1) * extent / float(4 * (n_layers + 1))) & (
coords_distances <= (2 * layer + 3) * extent / float(4 * (n_layers + 1)))]
img[tuple(np.transpose(layer_coords))] = layer_img[layer + 1][tuple(np.transpose(layer_coords))]
center_coords = coords[coords_distances <= extent / float(4 * (n_layers + 1))]
img[tuple(np.transpose(center_coords))] = 2
ext_coords = coords[coords_distances > (n_layers + 1) * extent / float(2 * (n_layers + 2))]
img[tuple(np.transpose(ext_coords))] = 1
from timagetk.io.image import default_image_attributes
from timagetk.util import now
img_kwargs = default_image_attributes(img, voxelsize=voxelsize)
img_kwargs["acquisition_date"] = now()
img = LabelledImage(img, not_a_label=0, **img_kwargs)
for p in points.keys():
# go from ZYX image convention to XYZ coordinates
points[p] = points[p][[2, 1, 0]]
_return = (img,)
if return_points:
_return += (points,)
if len(_return) == 1:
return _return[0]
return _return
