# -*- coding: utf-8 -*-
"""
MeasureRadius
=============
Measures intensity decays at specified points in the data.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>, Charly Rousseau <charly.rousseau@icm-institute.org>'
__license__ = 'GPLv3 - GNU General Public License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'https://idisco.info'
__download__ = 'https://www.github.com/ChristophKirst/ClearMap2'
import warnings
import numpy as np
import ClearMap.IO.IO as io
import ClearMap.ParallelProcessing.DataProcessing.MeasurePointList as mpl
import ClearMap.Utils.Timer as tmr
###############################################################################
### Measure Radius
###############################################################################
[docs]
def measure_radius(source, points, fraction=None, value=None,
max_radius=100, method='sphere', default=np.inf, scale=None,
return_radii=True, return_radii_as_scalar=True,
return_indices=False, processes=None,
verbose=False) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
"""
Measures a radius via decay of intensity values for a list of points.
Arguments
---------
source : array
Source for measurement.
points : array
List of indices to measure radis for.
fraction : float or None
Fraction of center intensity that needs to be reached to determine the
radius.
If None, value needs to be given.
value : array or float or None:
The value below which the intensity has to fall to measure the radius
from the center pixel. For array, it has to be the same size as the
points. If None, fraction has to be given.
max_radius : int or tuple of ints
The maximal pixel radius to consider in each dimension. The larger, the
slower the measurement.
default : number or None
Default value to use if no radius was detected.
scale: tuple or float
An optional scale in each direction to determine the distance.
return_radii : bool
If True, return the radii measured.
return_radii_as_scalar : bool
If True, return the radii as single floats, otherwise a radius for
each dimension.
return_indices : bool
If True, return the indices of the search which allows identifying the
pixel at which the radius condition was met.
processes : int or None
Number of processes to use.
verbose : bool
If True, print progress info.
Returns
-------
radii : array
Array of measured radii if return_radii is True.
indices : array
Array of measured indices at which the radius detection condition is met.
"""
if scale is None:
warnings.warn(f'No scale given, using 1 for all dimensions. '
f'Please note this behavior is deprecated. Future versions will require a scale.',
DeprecationWarning, stacklevel=2)
source = io.read(source)
if verbose:
timer = tmr.Timer()
print(f'Measuring radii of {points.shape[0]} points in array of shape {source.shape}.')
ndim = source.ndim
if not hasattr(max_radius, '__len__'):
max_radius = [max_radius] * ndim
if len(max_radius) != ndim:
raise ValueError(f'The maximal search radius {max_radius} has wrong dimension!')
if method == 'sphere':
search = search_indices_sphere(max_radius)
elif method == 'rectangle':
search = search_indices_rectangle(max_radius)
else:
raise ValueError(f'The method is not "sphere" or "rectangle" but {method}!')
worker_args = [source, points, search] + [value if value is not None else fraction]
worker_kwargs = {'sink': None, 'processes': processes, 'verbose': verbose}
if value is not None:
if hasattr(value, '__len__'):
measured = mpl.find_smaller_than_values(*worker_args, **worker_kwargs)
else:
measured = mpl.find_smaller_than_value(*worker_args, **worker_kwargs)
elif fraction is not None:
measured = mpl.find_smaller_than_fraction(*worker_args, **worker_kwargs)
else:
raise ValueError('fraction or value cannot both be None!')
if verbose:
timer.print_elapsed_time('Measuring radii done')
result = ()
if return_radii:
if scale is None:
scale = 1
if not hasattr(scale, '__len__'): # scalar provided -> broadcast
scale = np.full(ndim, scale, dtype=float)
distance_table = np.abs(search) * scale
if return_radii_as_scalar:
distance_table = np.sqrt(np.sum(distance_table * distance_table, axis=1))
distance_table = np.hstack([distance_table, default])
else:
distance_table = np.vstack([distance_table, [default] * ndim])
radii = distance_table[measured]
result += (radii,)
if return_indices:
search = np.vstack([search, [np.max(search, axis=0) + 1]])
indices = search[measured]
result += (indices,)
if len(result) == 1:
result = result[0]
return result
# ##############################################################################
# ## Search indices
# ##############################################################################
[docs]
def search_indices_sphere(radius):
"""
Creates all relative indices within a sphere of specified radius in an array with specified strides.
Arguments
---------
radius : tuple of int
Radius of the sphere of the search index list.
Returns
-------
indices : np.ndarray
Array of ints of relative indices for the search area voxels.
"""
grid = [np.arange(-r,r+1, dtype=float)/np.maximum(1,r) for r in radius]
grid = np.array(np.meshgrid(*grid, indexing = 'ij'))
#sort indices by radius
dist = np.sum(grid*grid, axis = 0)
dist_shape = dist.shape
dist = dist.reshape(-1)
dist_index = np.argsort(dist)
dist_sorted = dist[dist_index]
keep = dist_sorted <= 1
dist_index = dist_index[keep]
# convert to relative coordinates
indices = np.array(np.unravel_index(dist_index, dist_shape)).T
indices -= radius
return indices
[docs]
def search_indices_rectangle(radius):
"""
Creates all relative indices within a rectangle.
Arguments
---------
radius : tuple or float
Radius of the sphere of the search index list.
Returns
-------
indices : np.ndarray
Array of ints of relative indices for the search area voxels.
"""
# create coordinate grid
grid = [np.arange(-r,r+1, dtype=int) for r in radius]
grid = np.array(np.meshgrid(*grid, indexing = 'ij'))
# sort indices by radius
dist = np.sum(grid*grid, axis = 0)
dist_shape = dist.shape
dist = dist.reshape(-1)
dist_index = np.argsort(dist)
# convert to relative coordinates
indices = np.array(np.unravel_index(dist_index, dist_shape)).T
indices -= radius
return indices
###############################################################################
### Tests
###############################################################################
[docs]
def test():
import numpy as np
import ClearMap.IO.IO as io
# noinspection PyUnresolvedReferences
import ClearMap.Analysis.Measurements.radius_measurements as mr
data = 10-np.abs(10-np.arange(0,21))
search = mr.search_indices_sphere(radius=[10,10,10])
print(search)
points = np.array([10])
d,i = mr.measure_radius(data, points, fraction=0.75, max_radius=10, scale=2, verbose=True, processes=4, return_indices=True)
data = np.random.rand(*(30,40,50))
io.write('data.npy', data)
points = np.array([np.random.randint(0,s, size=10) for s in data.shape]).T
d,i = mr.measure_radius(data, points, value=0.5, max_radius=10, scale=2, verbose=True, processes=4, return_indices=True)
data = np.zeros((30,40,50), dtype=int)
data[10:20, 15:25,10:20] = 1
data[15,20,15] = 2
import ClearMap.Visualization.Plot3d as p3d
p3d.plot(data)
points = np.array([[15, 20, 15],[4,4,4]])
d,i = mr.measure_radius(data, points, value=0.0, max_radius=10, scale=None, verbose=True, processes=None, return_indices=True)
