Source code for ClearMap.ImageProcessing.Filter.Vector.VectorStructureElement

# -*- coding: utf-8 -*-
"""
VectorStructureElement
======================

Routines to generate vector valued structure elements for vector filters.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Public License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'

import numpy as np




def vectors_to_neighbours(shape, max_radius=None, reshape=True, orientation=False, rescale=True):
    """Generate vectors to voxels in a neighborhood."""

    # center
    center = tuple(s // 2 - 1 if s % 2 == 0 else s // 2 for s in shape)
    ranges = (np.arange(-c, c + 1) for c in center)

    if max_radius == 'sphere':
        max_radius = np.min(center)

    # all vectors
    vectors = np.array(np.meshgrid(*ranges, indexing='ij'), dtype=float)

    if reshape:
        vectors = vectors.reshape(len(shape), -1).T
    else:
        vectors = vectors.transpose(tuple(d for d in range(1, len(shape) + 1)) + (0,))

    # length weights # is this a good measure?
    lengths = np.linalg.norm(vectors, axis=-1)

    # valid
    valid = lengths > 0
    if max_radius is not None:
        valid = np.logical_and(valid, lengths <= max_radius)
        vectors[np.logical_not(valid)] = 0

    # orientation
    if orientation:
        # invert z < 0
        vectors[vectors[..., 2] < 0, :] *= -1
        # invert z = 0, y < 0
        ids0 = vectors[..., 2] == 0
        vectors[np.logical_and(ids0, vectors[..., 1] < 0), :] *= -1
        # invert z = 0, y = 0, x < 0
        ids0 = np.logical_and(ids0, vectors[..., 1] == 0)
        vectors[np.logical_and(ids0, vectors[..., 0] < 0), :] *= -1

    # rescale vectors
    if rescale:
        if rescale is True:
            def inverse(x):
                return 1.0/x
            rescale = inverse
        vectors[valid] = (vectors[valid] * rescale(lengths[valid][:, np.newaxis]))

    return vectors, lengths, valid