#!/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 timagetk.components.spatial_image import SpatialImage
from timagetk.synthetic_data.util import point_position_optimization
[docs]
def nuclei_image_from_point_positions(points, extent=50., voxelsize=(0.5, 0.25, 0.25), nuclei_radius=1.5,
nuclei_intensity=10000., point_signals=None, **kwargs):
"""Create a 3D nuclei image from a dictionary of 3D points.
Parameters
----------
points : dict
Dictionary with X,Y,Z 3D positions as values
extent : float, optional
The extent of the output image (in µm)
voxelsize : list(float), optional
Length-3 list of image Z,Y,X voxel-size (in µm)
nuclei_radius : float, optional
Radius of the spherical blobs representing nuclei (in µm)
nuclei_intensity : float, optional
The intensity value of the spherical blobs
point_signals : dict, optional
Dictionary with [0,1] intensity factors for each point
Other Parameters
----------------
dtype : {"uint8", "uint16"}
The bit-depth of the intensity image.
Returns
-------
timagetk.SpatialImage
The generated nuclei image
"""
image_size = [int(extent / v) for v in voxelsize]
x, y, z = np.ogrid[0:image_size[2] * voxelsize[2]:voxelsize[2], 0:image_size[1] * voxelsize[1]:voxelsize[1],
0:image_size[0] * voxelsize[0]:voxelsize[0]]
points_array = np.array(list(points.values()))
cell_distances = np.power(x[np.newaxis] - points_array[:, 0, np.newaxis, np.newaxis, np.newaxis], 2) + \
np.power(y[np.newaxis] - points_array[:, 1, np.newaxis, np.newaxis, np.newaxis], 2) + \
np.power(z[np.newaxis] - points_array[:, 2, np.newaxis, np.newaxis, np.newaxis], 2)
density_k = 2. / nuclei_radius
density_potential = (1. - np.tanh(density_k * (cell_distances - nuclei_radius))) / 2.
density = np.zeros_like(density_potential[0])
for i, p in enumerate(points.keys()):
if point_signals is None:
density += density_potential[i]
else:
if p in point_signals:
density += point_signals[p] * density_potential[i]
img = (nuclei_intensity * density).astype(kwargs.get('dtype', "uint16"))
img = np.transpose(img, (2, 1, 0))
from timagetk.io.image import default_image_attributes
from timagetk.util import now
im_kwargs = default_image_attributes(img, voxelsize=voxelsize)
im_kwargs['metadata'] = {'acquisition_date': now()}
return SpatialImage(img, **im_kwargs)
[docs]
def example_nuclei_image(n_points=100, extent=50., voxelsize=(0.5, 0.25, 0.25), nuclei_radius=1.5,
nuclei_intensity=10000., **kwargs):
"""Generate a synthetic 3D nuclei image.
The function creates a 3D image over a cubic volume of a given extent and populates it with randomly scattered nuclei.
Nuclei are represented by a spherical blobs of intensity of a same radius.
Parameters
----------
n_points : int, optional
Number of nuclei in the image.
Defaults to ``100``.
extent : float, optional
The extent of the output image (in µm).
Defaults to ``50.``.
voxelsize : list(float), optional
Length-3 list of image Z,Y,X voxelsize (in µm).
Defaults to ``(0.5, 0.25, 0.25)``.
nuclei_radius : float, optional
Radius of the spherical blobs representing nuclei (in µm).
Defaults to ``1.5``.
nuclei_intensity : float, optional
The intensity value of the spherical blobs.
Defaults to ``10000``.
Other Parameters
----------------
dtype : {"uint8", "uint16"}
The bit-depth of the intensity image.
Defaults to ``'uint16'``.
Returns
-------
timagetk.SpatialImage
The generated nuclei image.
dict
The XYZ-ordered 3D point positions used to generate the nuclei image.
Examples
--------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from timagetk.synthetic_data.nuclei_image import example_nuclei_image
>>> from timagetk.visu.mplt import image_plot
>>> nuclei_img, points = example_nuclei_image(n_points=10, extent=15., nuclei_radius=1.5, nuclei_intensity=10000)
>>> points = np.array(list(points.values()))
>>> plt.figure(figsize=(5, 5))
>>> plt.scatter(points[:,0], points[:,1], color='r') # Show X & Y coordinates
>>> image_plot(nuclei_img.max(axis=0),plt.gca(),val_range='auto')
>>> plt.show()
"""
image_size = [int(extent / v) for v in voxelsize]
center = (np.array(image_size) * np.array(voxelsize) / 2.)
points = {}
rng = np.random.default_rng()
for i in range(n_points):
points[i] = center + rng.random(3, dtype=float)
point_target_distance = np.power(n_points, 1 / 3.) * nuclei_radius * 1.5
sigma_deformation = nuclei_radius / 5.
points = point_position_optimization(points, target_distance=point_target_distance,
sigma_deformation=sigma_deformation, force_centering=True, center=center)
img = nuclei_image_from_point_positions(points, extent, voxelsize, nuclei_radius, nuclei_intensity,
dtype=kwargs.get('dtype', "uint16"))
return img, points
[docs]
def example_nuclei_signal_images(n_points=100, extent=50., voxelsize=(0.5, 0.25, 0.25), nuclei_radius=1.5,
signal_type='random', signal_intensity=10000., **kwargs):
"""Generate a synthetic 3D nuclei image with a second signal channel.
The function creates a 3D image over a cubic volume of a given extent and populates it with randomly scattered nuclei.
Nuclei are represented by a spherical blobs of intensity of a same radius.
In the second channel, a signal value is associated with each nucleus and multiplied with the intensity of the corresponding spherical blob.
The signal distribution can either be random or radially decreasing.
Parameters
----------
n_points : int, optional
Number of nuclei in the image.
Defaults to ``100``.
extent : float, optional
The extent of the output image (in µm).
Defaults to ``50.``.
voxelsize : list(float), optional
Len-3 list of image Z,Y,X voxelsize (in µm).
Defaults to ``(0.5, 0.25, 0.25)``.
nuclei_radius : float, optional
Radius of the spherical blobs representing nuclei (in µm).
Defaults to ``1.5``.
signal_type: {'random', 'center'}
Whether to use a radially decreasing ('center') or a random signal ('random').
Defaults to ``'random'``.
signal_intensity : float, optional
The intensity value of the spherical blobs.
Defaults to ``10000``.
Other Parameters
----------------
dtype : {"uint8", "uint16"}
The bit-depth of the intensity image.
Defaults to ``'uint16'``.
Returns
-------
timagetk.SpatialImage
The generated nuclei reference image
timagetk.SpatialImage
The generated signal intensity image
dict
The random X,Y,Z 3D point positions used to generate the nuclei
dict
The signal values associated with each nucleus
Examples
--------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from timagetk.synthetic_data.nuclei_image import example_nuclei_signal_images
>>> from timagetk.visu.mplt import image_plot
>>> nuclei_img, signal_img, points, point_signals = example_nuclei_signal_images(n_points=10, extent=15., nuclei_radius=1.5, signal_intensity=10000)
>>> points = np.array(list(points.values()))
>>> point_signals = np.array(list(point_signals.values()))
>>> fig, axes = plt.subplots(2, 2)
>>> fig.set_size_inches(w=4*5., h=5.)
>>> axes[0, 0].scatter(points[:,0], points[:,1], color='r') # Show X & Y coordinates
>>> axes[0, 1].scatter(points[:,0], points[:,1], c=point_signals, cmap='viridis') # Show X & Y coordinates
>>> image_plot(nuclei_img.max(axis=0),axes[1, 0],'gray',val_range='auto')
>>> image_plot(signal_img.max(axis=0),axes[1, 1],'gray',val_range='auto')
>>> plt.show()
"""
image_size = [int(extent / v) for v in voxelsize]
center = (np.array(image_size) * np.array(voxelsize) / 2.)
img, points = example_nuclei_image(n_points, extent, voxelsize, nuclei_radius, signal_intensity, **kwargs)
rng = np.random.default_rng()
if signal_type == 'random':
point_signals = {p: rng.random() for p in points.keys()}
elif signal_type == 'center':
point_signals = {p: 1. - np.linalg.norm(points[p] - center) / float(extent / 2) for p in points.keys()}
else:
raise ValueError(
f"Unknown `signal_type` value ('{signal_type}'), valid options are: {', '.join(['random', 'center'])}")
signal_img = nuclei_image_from_point_positions(points, extent, voxelsize, nuclei_radius, signal_intensity,
point_signals, **kwargs)
from timagetk.io.image import default_image_attributes
from timagetk.util import now
im_kwargs = default_image_attributes(img, voxelsize=voxelsize)
im_kwargs['metadata'] = {'acquisition_date': now()}
sim_kwargs = default_image_attributes(signal_img, voxelsize=voxelsize)
sim_kwargs['metadata'] = {'acquisition_date': now()}
return SpatialImage(img, **im_kwargs), SpatialImage(signal_img, **sim_kwargs), points, point_signals
