"""This module provides a complicated algorithm for making voxels out of regularly
gridded points. Considering that this algorithm is rather complex, we are keeping
it in its own module until we can simplify it, clean up the code, and make it
capable of handling non-uniformly gridded points
"""
__all__ = [
'VoxelizePoints',
]
__displayname__ = 'Voxelize'
import numpy as np
import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.util import numpy_support as nps
from .. import _helpers, interface
from ..base import FilterBase
from ..version import check_numpy
from .xyz import RotationTool
###############################################################################
[docs]
class VoxelizePoints(FilterBase):
"""This makes a ``vtkUnstructuredGrid`` of scattered points given voxel
sizes as input arrays. This assumes that the data is at least 2-Dimensional
on the XY Plane.
"""
__displayname__ = 'Voxelize Points'
__category__ = 'filter'
def __init__(self, **kwargs):
FilterBase.__init__(
self,
nInputPorts=1,
inputType='vtkPointSet',
nOutputPorts=1,
outputType='vtkUnstructuredGrid',
)
self.__dx = kwargs.get('dx', None)
self.__dy = kwargs.get('dy', None)
self.__dz = kwargs.get('dz', None)
self.__estimate_grid = kwargs.get('estimate', True)
self.__safe = kwargs.get('safe', 10.0)
self.__unique = kwargs.get('unique', True)
self.__tolerance = kwargs.get('tolerance', None)
self.__angle = kwargs.get('angle', 0.0)
[docs]
def add_field_data(self, grid):
"""An internal helper to add the recovered information as field data"""
# Add angle
a = vtk.vtkDoubleArray()
a.SetName('Recovered Angle (Deg.)')
a.SetNumberOfValues(1)
a.SetValue(0, np.rad2deg(self.__angle))
grid.GetFieldData().AddArray(a)
# Add cell sizes
s = vtk.vtkDoubleArray()
s.SetName('Recovered Cell Sizes')
s.SetNumberOfComponents(3)
s.InsertNextTuple3(self.__dx, self.__dy, self.__dz)
grid.GetFieldData().AddArray(s)
return grid
[docs]
@staticmethod
def add_cell_data(grid, arr, name):
"""Add a NumPy array as cell data to the given grid input"""
c = interface.convert_array(arr, name=name)
grid.GetCellData().AddArray(c)
return grid
[docs]
def points_to_grid(self, xo, yo, zo, dx, dy, dz, grid=None):
"""Convert XYZ points to a ``vtkUnstructuredGrid``."""
if not check_numpy(alert='warn'):
return grid
if grid is None:
grid = vtk.vtkUnstructuredGrid()
# TODO: Check dtypes on all arrays. Need to be floats
if self.__estimate_grid:
x, y, z = self.estimate_uniform_spacing(xo, yo, zo)
else:
x, y, z = xo, yo, zo
dx, dy, dz = self.__dx, self.__dy, self.__dz
if isinstance(dx, np.ndarray) and len(dx) != len(x):
raise _helpers.PVGeoError(
'X-Cell spacings are not properly defined for all points.'
)
if isinstance(dy, np.ndarray) and len(dy) != len(y):
raise _helpers.PVGeoError(
'Y-Cell spacings are not properly defined for all points.'
)
if isinstance(dz, np.ndarray) and len(dz) != len(z):
raise _helpers.PVGeoError(
'Z-Cell spacings are not properly defined for all points.'
)
n_cells = len(x)
# Generate cell nodes for all points in data set
# - Bottom
c_n1 = np.stack(((x - dx / 2), (y - dy / 2), (z - dz / 2)), axis=1)
c_n2 = np.stack(((x + dx / 2), (y - dy / 2), (z - dz / 2)), axis=1)
c_n3 = np.stack(((x - dx / 2), (y + dy / 2), (z - dz / 2)), axis=1)
c_n4 = np.stack(((x + dx / 2), (y + dy / 2), (z - dz / 2)), axis=1)
# - Top
c_n5 = np.stack(((x - dx / 2), (y - dy / 2), (z + dz / 2)), axis=1)
c_n6 = np.stack(((x + dx / 2), (y - dy / 2), (z + dz / 2)), axis=1)
c_n7 = np.stack(((x - dx / 2), (y + dy / 2), (z + dz / 2)), axis=1)
c_n8 = np.stack(((x + dx / 2), (y + dy / 2), (z + dz / 2)), axis=1)
# - Concatenate
all_nodes = np.concatenate(
(c_n1, c_n2, c_n3, c_n4, c_n5, c_n6, c_n7, c_n8), axis=0
)
pts = vtk.vtkPoints()
cells = vtk.vtkCellArray()
if self.__unique:
# Search for unique nodes and use the min cell size as the tolerance
if self.__tolerance is None:
TOLERANCE = np.min([dx, dy]) / 2.0
else:
TOLERANCE = self.__tolerance
# Round XY plane by the tolerance
txy = np.around(all_nodes[:, 0:2] / TOLERANCE)
all_nodes[:, 0:2] = txy
unique_nodes, ind_nodes = np.unique(all_nodes, return_inverse=True, axis=0)
unique_nodes[:, 0:2] *= TOLERANCE
all_nodes = unique_nodes
else:
ind_nodes = np.arange(0, len(all_nodes), dtype=int)
all_nodes[:, 0:2] = RotationTool.rotate(all_nodes[:, 0:2], -self.__angle)
if self.__estimate_grid:
self.add_field_data(grid)
# Add unique nodes as points in output
pts.SetData(interface.convert_array(all_nodes))
# Add cell vertices
j = np.multiply(np.tile(np.arange(0, 8, 1), n_cells), n_cells)
arridx = np.add(j, np.repeat(np.arange(0, n_cells, 1, dtype=int), 8))
ids = ind_nodes[arridx].reshape((n_cells, 8))
cells_mat = np.concatenate(
(np.ones((ids.shape[0], 1), dtype=np.int_) * ids.shape[1], ids), axis=1
).ravel()
cells = vtk.vtkCellArray()
cells.SetNumberOfCells(n_cells)
cells.SetCells(
n_cells, nps.numpy_to_vtk(cells_mat, deep=True, array_type=vtk.VTK_ID_TYPE)
)
# Set the output
grid.SetPoints(pts)
grid.SetCells(vtk.VTK_VOXEL, cells)
return grid
[docs]
@staticmethod
def _copy_arrays(pdi, pdo):
"""internal helper to copy arrays from point data to cell data in the voxels."""
for i in range(pdi.GetPointData().GetNumberOfArrays()):
arr = pdi.GetPointData().GetArray(i)
_helpers.add_array(pdo, 1, arr) # adds to CELL data
return pdo
[docs]
def RequestData(self, request, inInfo, outInfo):
"""Used by pipeline to generate output"""
# Get input/output of Proxy
pdi = self.GetInputData(inInfo, 0, 0)
pdo = self.GetOutputData(outInfo, 0)
# Perform task
wpdi = dsa.WrapDataObject(pdi)
pts = wpdi.Points
x, y, z = pts[:, 0], pts[:, 1], pts[:, 2]
self.points_to_grid(x, y, z, self.__dx, self.__dy, self.__dz, grid=pdo)
# Now append data to grid
self._copy_arrays(pdi, pdo)
return 1
#### Setters and Getters ####
[docs]
def set_safe_size(self, safe):
"""A voxel size to use if a spacing cannot be determined for an axis"""
if self.__safe != safe:
self.__safe = safe
self.Modified()
[docs]
def set_delta_x(self, dx):
"""Set the X cells spacing
Args:
dx (float or np.array(floats)): the spacing(s) for the cells in
the X-direction
"""
self.__dx = dx
self.Modified()
[docs]
def set_delta_y(self, dy):
"""Set the Y cells spacing
Args:
dy (float or np.array(floats)): the spacing(s) for the cells in
the Y-direction
"""
self.__dy = dy
self.Modified()
[docs]
def set_delta_z(self, dz):
"""Set the Z cells spacing
Args:
dz (float or np.array(floats)): the spacing(s) for the cells in
the Z-direction
"""
self.__dz = dz
self.set_safe_size(np.min(dz))
self.Modified()
[docs]
def set_deltas(self, dx, dy, dz):
"""Set the cell spacings for each axial direction
Args:
dx (float or np.array(floats)): the spacing(s) for the cells in
the X-direction
dy (float or np.array(floats)): the spacing(s) for the cells in
the Y-direction
dz (float or np.array(floats)): the spacing(s) for the cells in
the Z-direction
"""
self.set_delta_x(dx)
self.set_delta_y(dy)
self.set_delta_z(dz)
[docs]
def set_estimate_grid(self, flag):
"""Set a flag on whether or not to estimate the grid spacing/rotation"""
if self.__estimate_grid != flag:
self.__estimate_grid = flag
self.Modified()
[docs]
def set_unique(self, flag):
"""Set a flag on whether or not to try to elimate non unique elements"""
if self.__unique != flag:
self.__unique = flag
self.Modified()
[docs]
def get_angle(self, degrees=True):
"""Returns the recovered angle if set to recover the input grid. If the
input points are rotated, then this angle will reflect a close
approximation of that rotation.
Args:
degrees (bool): A flag on to return decimal degrees or radians.
"""
if degrees:
return np.rad2deg(self.__angle)
return self.__angle
[docs]
def get_recovered_angle(self, degrees=True):
"""DEPRECATED: use `get_angle`"""
return self.get_angle(degrees=degrees)
[docs]
def set_angle(self, angle):
"""Set the rotation angle manually"""
if self.__angle != angle:
self.__angle = angle
self.Modified()
[docs]
def get_spacing(self):
"""Get the cell spacings"""
return (self.__dx, self.__dy, self.__dz)
###############################################################################