Source code for ClearMap.Analysis.Measurements.Voxelization

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

Converts point data into a voxelized image for visulalization and analysis.
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'

import numpy as np

import ClearMap.IO.IO as io

import ClearMap.ParallelProcessing.DataProcessing.DevolvePointList as dpl


###############################################################################
### Voxelization
###############################################################################

[docs] def voxelize(source, sink = None, shape = None, dtype = None, weights = None, method = 'sphere', radius = (1,1,1), kernel = None, processes = None, verbose = False): """Converts a list of points into an volumetric image array Arguments --------- source : str, array or Source Source of point of nxd coordinates. sink : str, array or None The sink for the voxelized image, if None return array. shape : tuple or None Shape of the final voxelized data. If None, deterimine from points. dtype : dtype or None Optional data type of the sink. weights : array or None Weight array of length n for each point. If None, use uniform weights. method : str Method for voxelization: 'sphere', 'rectangle' or 'pixel'. radius : tuple Radius of the voxel region to integrate over. kernel : function Optional function of distance to set weights in the voxelization. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : str, array Volumetric data of voxelied point data. """ points = io.read(source); points_shape = points.shape; if len(points_shape) > 1: ndim = points_shape[1]; else: ndim = 1; if not hasattr(radius, '__len__'): radius = [radius] * ndim; if len(radius) != ndim: raise ValueError('Radius %r and points with shape %r do not match in dimension!' % (radius, points_shape)); if method == 'sphere': indices, kernel = search_indices_sphere(radius, kernel) elif method == 'rectangle': indices, kernel = search_indices_rectangle(radius, kernel) elif method == 'pixel': indices = np.array(0, dtype=int); if kernel is not None: kernel = np.array([kernel(0)]) else: raise ValueError("method not 'sphere', 'rectangle', or 'pixel', but %r!" % method) return dpl.devolve(points, sink=sink, shape=shape, dtype=dtype, weights=weights, indices=indices, kernel=kernel, processes=processes, verbose=verbose);
############################################################################### ### Search indices ###############################################################################
[docs] def search_indices_sphere(radius, kernel = None): """Creates all relative indices within a sphere of specified radius. Arguments --------- radius : tuple or int Radius of the sphere of the search index list. Returns ------- indices : array Array of ints of relative indices for the search area voxels. """ #create coordiante grid 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]; if kernel is not None: dist_sorted = np.sqrt(dist_sorted[keep]); kernel = np.array([kernel(d) for d in dist_sorted], dtype=float); else: kernel = None; # convert to relative coordinates indices = np.array(np.unravel_index(dist_index, dist_shape)).T; indices -= radius; return indices, kernel;
[docs] def search_indices_rectangle(radius, kernel = None): """Creates all relative indices within a rectangle. Arguments --------- radius : tuple or float Radius of the sphere of the search index list. Returns ------- indices : array Array of ints of relative indices for the search area voxels. """ #create coordiante grid grid = [np.arange(-r,r+1, dtype=int) for r in radius]; grid = np.array(np.meshgrid(*grid, indexing = 'ij')); if kernel is not None: dist = np.sqrt(np.sum(grid*grid, axis =0)); kernel = np.array([kernel(d) for d in dist], dtype=float); else: kernel = None; indices = grid.reshape((len(radius),-1)).T; return indices, kernel;
############################################################################### ### Tests ############################################################################### def _test(): """Tests""" import numpy as np import ClearMap.Visualization.Plot3d as p3d import ClearMap.Analysis.Measurements.Voxelization as vox #points = np.random.rand(20,3) * 15; #points = np.asarray(points, dtype=int); points = np.array([[10,1,1]]) v = vox.voxelize(points, shape = (20,20,20), radius = (0,0,0)); w = vox.voxelize(points, shape = (20,20,20), radius = (2,2,2)); p3d.plot([[v,w]]) indices, kernel = vox.search_indices_sphere((2,2,2), kernel=None) print(indices) points = np.array([[10,10,18]]) v = vox.voxelize(points, shape = (20,20,20), weights=None, radius=(5,7,10)); p3d.plot(v) points = np.array([[10,10,10]]); weights = np.random.rand(len(points)); v = vox.voxelize(points, shape = (20,20,20), weights=weights, radius=(2,2,2), method = 'rectangle'); p3d.plot(v) def kernel(d): return np.exp(-d); v = vox.voxelize(points, shape = (20,20,20), weights=None, radius=(5,5,8), kernel=kernel, method = 'sphere'); p3d.plot(v)