"""
Vector
======
Module to compute filters on vector/direction fields
Used for filtering orientation and fiber data.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Public License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2023 by Christoph Kirst'
import numpy as np
from ..StructureElement import sphere
from .VectorStructureElement import vectors_to_neighbours
import pyximport
pyximport.install(setup_args={"include_dirs": np.get_include()}, reload_support=True)
from . import VectorCode as code
########################################################################################################################
# Vector filters
########################################################################################################################
[docs]
def similarity(vectors, structural_element, structural_element_offset=None, weights=None, sink=None):
"""Compute sum of scalar products of the center vector with neighbouring vectors defined by the structural element.
Arguments
---------
vectors : array
Input vector field array.
structural_element: array or str
Structural element specification.
structural_element_offset : tuple or None
Structural element offset
weights : None or array
Weights for each vector in the vector field.
sink : array
Output, if None, a new array is allocated.
Returns
-------
sink : array
Filtered output.
Algorithm
---------
For each center voxel compute the scalar products of the center vector with the neighbouring vectors defined by the
structural element and sum.
"""
if weights is None:
return _apply_code(code.similarity, vectors, structural_element, structural_element_offset, weights, sink)
else:
return _apply_code(code.similarity_weighted, vectors, structural_element, structural_element_offset, weights,
sink)
[docs]
def alignment(vectors, structural_element, structural_element_offset=None, weights=None, sink=None):
"""Compute alignment of central vector with neighboring vectors.
Arguments
---------
vectors : array
Input vector field array.
structural_element: array or str
Structural element specification.
structural_element_offset : tuple or None
Structural element offset
weights : None or array
Weights for each vector in the vector field.
sink : array
Output, if None, a new array is allocated.
Returns
-------
sink : array
Filtered output.
Algorithm
---------
For each central voxel
1. compute a set of weights to neighboring vectors via the scalar product of the
central vector with the structural element vectors
2. compute the scalar product of the central vector with the neighbouring vectors
3. from the weighed sum
"""
if weights is None:
return _apply_code(code.alignment, vectors,
structural_element, structural_element_offset,
weights, sink)
else:
return _apply_code(code.alignment_weighted, vectors,
structural_element, structural_element_offset,
weights, sink)
[docs]
def orientation(vectors, structural_element, structural_element_offset=None, weights=None, sink=None):
"""Compute orientation alignment of central vector with neighboring vectors.
Arguments
---------
vectors : array
Input vector field array.
structural_element: array or str
Structural element specification.
structural_element_offset : tuple or None
Structural element offset
weights : None or array
Weights for each vector in the vector field.
sink : array
Output, if None, a new array is allocated.
Returns
-------
sink : array
Filtered output.
Algorithm
---------
For each central voxel
1. compute a set of matched orientation weights to neighboring vectors via the absolute value scalar product
of the central vector with the structural element vectors (e.g. direction vectors from center), and further
multiply by the input weights for the vectors.
2. compute the absolute value of the scalar product of the central vector with the neighbouring vectors
3. from the weighed sum
"""
if weights is None:
return _apply_code(code.orientation, vectors,
structural_element, structural_element_offset,
weights, sink)
else:
return _apply_code(code.orientation_weighted, vectors,
structural_element, structural_element_offset,
weights,
sink)
[docs]
def orientation_ellipsoid(vectors, weights=None, max_radius=5, scales=(2,1), sink=None):
"""Compute orientation alignment of central vector with neighboring vectors.
Arguments
---------
vectors : array
Input vector field array.
weights : None or array
Weights for each vector in the vector field.
max_radius : float
The maximal radius of the ellipse
scales : tuple or list
A tuple (s1,s2) with s1 for the scale of the mayor axis and s2 the scale for the 2 minor axes of the ellipsoid.
sink : array
Output, if None, a new array is allocated.
Returns
-------
sink : array
Filtered output.
Algorithm
---------
For each central voxel
1. compute a set of matched orientation weights to neighboring vectors via the absolute value scalar product
of the central vector with the structural element vectors (e.g. direction vectors from center), and further
multiply by the input weights for the vectors.
2. compute the absolute value of the scalar product of the central vector with the neighbouring vectors
3. from the weighed sum
"""
if not hasattr(scales, '__len__') or len(scales) != 2:
raise ValueError('scales must be a tuple or list of two elements, found %r' % (scales,))
radius = max_radius / np.max(scales)
parameter = tuple(scales) + (radius,)
# compute structural element that fits ellipsoid
shape = (int(radius * np.max(scales))-1,) * 3
structural_element = sphere(shape)
structural_element_offset = None
if weights is None:
return _apply_code(code.orientation_ellipsoid, vectors,
structural_element, structural_element_offset,
weights, sink, parameter=parameter)
else:
return _apply_code(code.orientation_ellipsoid_weighted, vectors, structural_element, structural_element_offset,
weights, sink, parameter=parameter)
########################################################################################################################
# Helpers
########################################################################################################################
def _apply_code(function, vectors, structural_element, structural_element_offset, weights,
sink=None, sink_dtype=None, sink_shape_per_pixel=None, parameter=None):
"""Helper to apply the core functions"""
if vectors.ndim != 4:
raise ValueError('The vector field is assumed to be 4d array, found %dd!' % (vectors.ndim,))
shape = vectors.shape[:-1]
if isinstance(structural_element, tuple) and len(structural_element) == 2 and \
structural_element[0] == 'relative_distance':
structural_element_shape = (structural_element[1],) * 3
structural_element, _, _ = vectors_to_neighbours(structural_element_shape, reshape=False,
rescale=True, max_radius='sphere')
if structural_element.ndim == 3:
structural_element = structural_element[..., np.newaxis]
structural_element_index = np.array(np.array(np.where(np.sum(structural_element != 0, axis=-1))).T)
structural_element_vector = structural_element[structural_element_index[:, 0],
structural_element_index[:, 1],
structural_element_index[:, 2]]
if structural_element_offset is None:
structural_element_offset = np.array(structural_element.shape[:3]) // 2
structural_element_offset = np.array(structural_element_offset, dtype=int)
if weights is not None:
if weights.shape != shape:
raise ValueError(
'The weights do not match the shape of the vector field %r != %r!' % (weights.shape, shape))
if sink_shape_per_pixel is None:
shape_per_pixel = (1,)
else:
shape_per_pixel = sink_shape_per_pixel
if sink is None:
if sink_dtype is None:
sink_dtype = float
sink = np.zeros(shape + shape_per_pixel, dtype=sink_dtype, order='F')
else:
if shape_per_pixel != (1,):
if sink.shape != shape + shape_per_pixel:
raise ValueError(
'The sink of shape %r does not have expected shape %r!' % (sink.shape, shape + shape_per_pixel))
if sink.shape != shape + shape_per_pixel:
sink = sink.reshape(shape + shape_per_pixel)
if structural_element_vector.dtype != sink_dtype:
structural_element_vector = np.asarray(structural_element_vector, dtype=sink_dtype)
if sink.dtype == bool:
sink_view = sink.view('uint8')
else:
sink_view = sink
if parameter is None:
parameter = np.zeros(0)
parameter = np.asarray([parameter], dtype=float).flatten()
if weights is None:
function(vectors=vectors,
selem_index=structural_element_index,
selem_vector=structural_element_vector,
selem_offset=structural_element_offset,
sink=sink_view,
parameter=parameter)
else:
function(vectors=vectors,
weights=weights,
selem_index=structural_element_index,
selem_vector=structural_element_vector,
selem_offset=structural_element_offset,
sink=sink_view,
parameter=parameter)
if sink_shape_per_pixel is None:
sink = sink.reshape(sink.shape[:-1])
return sink
########################################################################################################################
# Tests
########################################################################################################################
def _test():
# %gui qt
import numpy as np
import ClearMap.ImageProcessing.Filter.Vector.Vector as vec
import ClearMap.ImageProcessing.Filter.Vector.VectorStructureElement as vse
import ClearMap.ImageProcessing.Filter.StructureElement as se
import ClearMap.Visualization.PyVista.plot as pvp
import ClearMap.Visualization.Plot3d as p3d
print('done')
from importlib import reload
reload(vec.code)
reload(vec)
# axon along x axis
vectors = np.zeros((10, 10, 10, 3))
vectors[:, 4:7, 4:7, 0] = 1
vectors[:5] *= -1
# similarity
selem = se.sphere((3,3,3)) > 0
filtered = vec.similarity(vectors, selem)
#pvp.plot_3d_volume(filtered)
p3d.plot({'source' : filtered, 'vectors': vectors})
# orientation ellipsoid filter
filtered = vec.orientation_ellipsoid(vectors, scales=(5,0.5), max_radius=10)
p3d.plot({'source': filtered, 'vectors': vectors})
p3d.close(all)
# real data
filtered = vec.orientation_ellipsoid(o, scales=(4, 1), max_radius=6, weights=a3)
p3d.plot([a3, filtered, {'source': filtered, 'orientations': o}])
# structural_element
def rescale(x):
return 1.0/x**2
structure, lengths, valid = vse.vectors_to_neighbours((7, 7, 7), reshape=False, rescale=rescale, max_radius='sphere')
p = pvp.plot_volumetric_vector_field(structure)
p.show_grid()
# orientation weights
central_vector = np.array([1,1,0])
weights = np.abs(structure @ central_vector) / np.linalg.norm(central_vector)
nrm = np.linalg.norm(structure, axis=-1)
valid = nrm > 0
weights[valid] = weights[valid] / nrm[valid]
weights = weights**5
weights = weights * nrm
pvp.plot_3d_volume(weights)
p3d.plot(weights)
p3d.plot(nrm)
filtered = vec.alignment(vectors=vectors, selem=selem)
print('done')
# %gui qt
p = p3d.plot_volumetric_vector_field(selem)
p.show_grid()
p = p3d.plot_volumetric_vector_field(vectors)
p.show_grid()
# %gui qt
p3d.plot_3d_volume(filtered)
weights = np.random.rand(*vectors.shape[:-1])
filtered = vec.alignment(vectors=vectors, selem=selem, weights=weights)
print('done')
