+"""
+Utilities for performing Free Form Deformation (FFD).
+
+:Theoretical Insight:
+
+ Free Form Deformation is a technique for the efficient, smooth and accurate
+ geometrical parametrization. It has been proposed the first time in
+ *Sederberg, Thomas W., and Scott R. Parry. "Free-form deformation of solid
+ geometric models." ACM SIGGRAPH computer graphics 20.4 (1986): 151-160*. It
+ consists in three different step:
+
+ - Mapping the physical domain to the reference one with map
+ :math:`\\boldsymbol{\\psi}`. In the code it is named *transformation*.
+
+ - Moving some control points to deform the lattice with :math:`\\hat{T}`.
+ The movement of the control points is basically the weight (or
+ displacement) :math:`\\boldsymbol{\\mu}` we set in the *parameters file*.
+
+ - Mapping back to the physical domain with map
+ :math:`\\boldsymbol{\\psi}^{-1}`. In the code it is named
+ *inverse_transformation*.
+
+ FFD map (:math:`T`) is the composition of the three maps, that is
+
+ .. math:: T(\\cdot, \\boldsymbol{\\mu}) = (\\Psi^{-1} \\circ \\hat{T} \\circ
+ \\Psi) (\\cdot, \\boldsymbol{\\mu})
+
+ In this way, every point inside the FFD box changes it position according to
+
+ .. math:: \\boldsymbol{P} = \\boldsymbol{\\psi}^{-1} \\left( \\sum_{l=0}^L
+ \\sum_{m=0}^M \\sum_{n=0}^N
+ \\mathsf{b}_{lmn}(\\boldsymbol{\\psi}(\\boldsymbol{P}_0))
+ \\boldsymbol{\\mu}_{lmn} \\right)
+
+ where :math:`\\mathsf{b}_{lmn}` are Bernstein polynomials. We improve the
+ traditional version by allowing a rotation of the FFD lattice in order to
+ give more flexibility to the tool.
+
+ You can try to add more shapes to the lattice to allow more and more
+ involved transformations.
+
+"""
+try:
+ import configparser as configparser
+except ImportError:
+ import ConfigParser as configparser
+import os
+import copy
+import numpy as np
+from scipy import special
+
+from pygem import Deformation
+from pygem.utils import fit_affine_transformation, angles2matrix
+
+
+
+
[docs]
+
class FFD(Deformation):
+
"""
+
Class that handles the Free Form Deformation on the mesh points.
+
+
:param list n_control_points: number of control points in the x, y, and z
+
direction. Default is [2, 2, 2].
+
+
:cvar numpy.ndarray box_length: dimension of the FFD bounding box, in the
+
x, y and z direction (local coordinate system).
+
:cvar numpy.ndarray box_origin: the x, y and z coordinates of the origin of
+
the FFD bounding box.
+
:cvar numpy.ndarray rot_angle: rotation angle around x, y and z axis of the
+
FFD bounding box.
+
:cvar numpy.ndarray n_control_points: the number of control points in the
+
x, y, and z direction.
+
:cvar numpy.ndarray array_mu_x: collects the displacements (weights) along
+
x, normalized with the box length x.
+
:cvar numpy.ndarray array_mu_y: collects the displacements (weights) along
+
y, normalized with the box length y.
+
:cvar numpy.ndarray array_mu_z: collects the displacements (weights) along
+
z, normalized with the box length z.
+
+
:Example:
+
+
>>> from pygem import FFD
+
>>> import numpy as np
+
>>> ffd = FFD()
+
>>> ffd.read_parameters(
+
>>> 'tests/test_datasets/parameters_test_ffd_sphere.prm')
+
>>> original_mesh_points = np.load(
+
>>> 'tests/test_datasets/meshpoints_sphere_orig.npy')
+
>>> new_mesh_points = ffd(original_mesh_points)
+
"""
+
reference_frame = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]])
+
+
def __init__(self, n_control_points=None):
+
self.conversion_unit = 1. # TODO: unused at the moment
+
+
self.box_length = np.array([1., 1., 1.])
+
self.box_origin = np.array([0., 0., 0.])
+
self.rot_angle = np.array([0., 0., 0.])
+
+
self.array_mu_x = None
+
self.array_mu_y = None
+
self.array_mu_z = None
+
+
if n_control_points is None:
+
n_control_points = [2, 2, 2]
+
self.n_control_points = n_control_points
+
+
@property
+
def n_control_points(self):
+
"""
+
The number of control points in X, Y and Z directions
+
+
:rtype: numpy.ndarray
+
"""
+
return self._n_control_points
+
+
@n_control_points.setter
+
def n_control_points(self, npts):
+
self._n_control_points = np.array(npts)
+
self.array_mu_x = np.zeros(self.n_control_points)
+
self.array_mu_y = np.zeros(self.n_control_points)
+
self.array_mu_z = np.zeros(self.n_control_points)
+
+
@property
+
def psi(self):
+
"""
+
Return the function that map the physical domain to the reference
+
domain.
+
+
:rtype: callable
+
"""
+
physical_frame = self.position_vertices - self.box_origin
+
return fit_affine_transformation(physical_frame, self.reference_frame)
+
+
@property
+
def inverse_psi(self):
+
"""
+
Return the function that map the reference domain to the physical
+
domain.
+
+
:rtype: callable
+
"""
+
physical_frame = self.position_vertices - self.box_origin
+
return fit_affine_transformation(self.reference_frame, physical_frame)
+
+
@property
+
def T(self):
+
"""
+
Return the function that deforms the points within the unit cube.
+
+
:rtype: callable
+
"""
+
def T_mapping(points):
+
(n_rows, n_cols) = points.shape
+
(dim_n_mu, dim_m_mu, dim_t_mu) = self.array_mu_x.shape
+
+
# Initialization. In order to exploit the contiguity in memory the
+
# following are transposed
+
bernstein_x = np.zeros((dim_n_mu, n_rows))
+
bernstein_y = np.zeros((dim_m_mu, n_rows))
+
bernstein_z = np.zeros((dim_t_mu, n_rows))
+
shift_points = np.zeros((n_cols, n_rows))
+
+
# TODO check no-loop implementation
+
#bernstein_x = (
+
# np.power(mesh_points[:, 0][:, None], range(dim_n_mu)) *
+
# np.power(1 - mesh_points[:, 0][:, None], range(dim_n_mu-1, -1, -1)) *
+
# special.binom(np.array([dim_n_mu-1]*dim_n_mu), np.arange(dim_n_mu))
+
#)
+
for i in range(0, dim_n_mu):
+
aux1 = np.power((1 - points[:, 0]), dim_n_mu - 1 - i)
+
aux2 = np.power(points[:, 0], i)
+
bernstein_x[i, :] = (special.binom(dim_n_mu - 1, i) *
+
np.multiply(aux1, aux2))
+
+
for i in range(0, dim_m_mu):
+
aux1 = np.power((1 - points[:, 1]), dim_m_mu - 1 - i)
+
aux2 = np.power(points[:, 1], i)
+
bernstein_y[i, :] = special.binom(dim_m_mu - 1,
+
i) * np.multiply(aux1, aux2)
+
+
for i in range(0, dim_t_mu):
+
aux1 = np.power((1 - points[:, 2]), dim_t_mu - 1 - i)
+
aux2 = np.power(points[:, 2], i)
+
bernstein_z[i, :] = special.binom(dim_t_mu - 1,
+
i) * np.multiply(aux1, aux2)
+
+
aux_x = 0.
+
aux_y = 0.
+
aux_z = 0.
+
+
for j in range(0, dim_m_mu):
+
for k in range(0, dim_t_mu):
+
bernstein_yz = np.multiply(bernstein_y[j, :],
+
bernstein_z[k, :])
+
for i in range(0, dim_n_mu):
+
aux = np.multiply(bernstein_x[i, :], bernstein_yz)
+
aux_x += aux * self.array_mu_x[i, j, k]
+
aux_y += aux * self.array_mu_y[i, j, k]
+
aux_z += aux * self.array_mu_z[i, j, k]
+
+
shift_points[0, :] += aux_x
+
shift_points[1, :] += aux_y
+
shift_points[2, :] += aux_z
+
return shift_points.T + points
+
+
return T_mapping
+
+
@property
+
def rotation_matrix(self):
+
"""
+
The rotation matrix (according to rot_angle_x, rot_angle_y,
+
rot_angle_z).
+
+
:rtype: numpy.ndarray
+
"""
+
return angles2matrix(np.radians(self.rot_angle[2]),
+
np.radians(self.rot_angle[1]),
+
np.radians(self.rot_angle[0]))
+
+
@property
+
def position_vertices(self):
+
"""
+
The position of the vertices of the FFD bounding box.
+
+
:rtype: numpy.ndarray
+
"""
+
return self.box_origin + np.vstack([
+
np.zeros((1, 3)),
+
self.rotation_matrix.dot(np.diag(self.box_length)).T
+
])
+
+
+
[docs]
+
def reset_weights(self):
+
"""
+
Set transformation parameters to arrays of zeros.
+
"""
+
self.array_mu_x.fill(0.0)
+
self.array_mu_y.fill(0.0)
+
self.array_mu_z.fill(0.0)
+
+
+
+
[docs]
+
def read_parameters(self, filename='parameters.prm'):
+
"""
+
Reads in the parameters file and fill the self structure.
+
+
:param string filename: parameters file to be read in.
+
"""
+
if not isinstance(filename, str):
+
raise TypeError("filename must be a string")
+
+
# Checks if the parameters file exists. If not it writes the default
+
# class into filename.
+
if not os.path.isfile(filename):
+
self.write_parameters(filename)
+
return
+
+
config = configparser.RawConfigParser()
+
config.read(filename)
+
+
self.n_control_points[0] = config.getint('Box info',
+
'n control points x')
+
self.n_control_points[1] = config.getint('Box info',
+
'n control points y')
+
self.n_control_points[2] = config.getint('Box info',
+
'n control points z')
+
+
self.box_length[0] = config.getfloat('Box info', 'box length x')
+
self.box_length[1] = config.getfloat('Box info', 'box length y')
+
self.box_length[2] = config.getfloat('Box info', 'box length z')
+
+
self.box_origin[0] = config.getfloat('Box info', 'box origin x')
+
self.box_origin[1] = config.getfloat('Box info', 'box origin y')
+
self.box_origin[2] = config.getfloat('Box info', 'box origin z')
+
+
self.rot_angle[0] = config.getfloat('Box info', 'rotation angle x')
+
self.rot_angle[1] = config.getfloat('Box info', 'rotation angle y')
+
self.rot_angle[2] = config.getfloat('Box info', 'rotation angle z')
+
+
self.array_mu_x = np.zeros(self.n_control_points)
+
self.array_mu_y = np.zeros(self.n_control_points)
+
self.array_mu_z = np.zeros(self.n_control_points)
+
+
mux = config.get('Parameters weights', 'parameter x')
+
muy = config.get('Parameters weights', 'parameter y')
+
muz = config.get('Parameters weights', 'parameter z')
+
+
for line in mux.split('\n'):
+
values = np.array(line.split())
+
self.array_mu_x[tuple(map(int, values[0:3]))] = float(values[3])
+
+
for line in muy.split('\n'):
+
values = line.split()
+
self.array_mu_y[tuple(map(int, values[0:3]))] = float(values[3])
+
+
for line in muz.split('\n'):
+
values = line.split()
+
self.array_mu_z[tuple(map(int, values[0:3]))] = float(values[3])
+
+
+
+
[docs]
+
def write_parameters(self, filename='parameters.prm'):
+
"""
+
This method writes a parameters file (.prm) called `filename` and fills
+
it with all the parameters class members.
+
+
:param string filename: parameters file to be written out.
+
"""
+
if not isinstance(filename, str):
+
raise TypeError("filename must be a string")
+
+
output_string = ""
+
output_string += '\n[Box info]\n'
+
output_string += '# This section collects all the properties of the'
+
output_string += ' FFD bounding box.\n'
+
+
output_string += '\n# n control points indicates the number of control'
+
output_string += ' points in each direction (x, y, z).\n'
+
output_string += '# For example, to create a 2 x 3 x 2 grid, use the'
+
output_string += ' following: n control points: 2, 3, 2\n'
+
output_string += 'n control points x: ' + str(
+
self.n_control_points[0]) + '\n'
+
output_string += 'n control points y: ' + str(
+
self.n_control_points[1]) + '\n'
+
output_string += 'n control points z: ' + str(
+
self.n_control_points[2]) + '\n'
+
+
output_string += '\n# box length indicates the length of the FFD '
+
output_string += 'bounding box along the three canonical directions '
+
output_string += '(x, y, z).\n'
+
+
output_string += '# It uses the local coordinate system.\n'
+
output_string += '# For example to create a 2 x 1.5 x 3 meters box '
+
output_string += 'use the following: box length: 2.0, 1.5, 3.0\n'
+
+
output_string += 'box length x: ' + str(self.box_length[0]) + '\n'
+
output_string += 'box length y: ' + str(self.box_length[1]) + '\n'
+
output_string += 'box length z: ' + str(self.box_length[2]) + '\n'
+
+
output_string += '\n# box origin indicates the x, y, and z coordinates '
+
output_string += 'of the origin of the FFD bounding box. That is '
+
output_string += 'center of\n'
+
+
output_string += '# rotation of the bounding box. It corresponds to '
+
output_string += 'the point coordinates with position [0][0][0].\n'
+
+
output_string += '# See section "Parameters weights" for more '
+
output_string += 'details.\n'
+
output_string += '# For example, if the origin is equal to 0., 0., 0., '
+
output_string += 'use the following: box origin: 0., 0., 0.\n'
+
+
output_string += 'box origin x: ' + str(self.box_origin[0]) + '\n'
+
output_string += 'box origin y: ' + str(self.box_origin[1]) + '\n'
+
output_string += 'box origin z: ' + str(self.box_origin[2]) + '\n'
+
+
output_string += '\n# rotation angle indicates the rotation angle '
+
output_string += 'around the x, y, and z axis of the FFD bounding box '
+
output_string += 'in degrees.\n'
+
+
output_string += '# The rotation is done with respect to the box '
+
output_string += 'origin.\n'
+
output_string += '# For example, to rotate the box by 2 deg along '
+
output_string += 'the z '
+
output_string += 'direction, use the following: rotation angle: '
+
output_string += '0., 0., 2.\n'
+
+
output_string += 'rotation angle x: ' + str(self.rot_angle[0]) + '\n'
+
output_string += 'rotation angle y: ' + str(self.rot_angle[1]) + '\n'
+
output_string += 'rotation angle z: ' + str(self.rot_angle[2]) + '\n'
+
+
output_string += '\n\n[Parameters weights]\n'
+
output_string += '# This section describes the weights of the FFD '
+
output_string += 'control points.\n'
+
+
output_string += '# We adopt the following convention:\n'
+
output_string += '# For example with a 2x2x2 grid of control points we '
+
output_string += 'have to fill a 2x2x2 matrix of weights.\n'
+
+
output_string += '# If a weight is equal to zero you can discard the '
+
output_string += 'line since the default is zero.\n'
+
+
output_string += '#\n'
+
output_string += '# | x index | y index | z index | weight |\n'
+
output_string += '# --------------------------------------\n'
+
output_string += '# | 0 | 0 | 0 | 1.0 |\n'
+
output_string += '# | 0 | 1 | 1 | 0.0 | --> you '
+
output_string += 'can erase this line without effects\n'
+
output_string += '# | 0 | 1 | 0 | -2.1 |\n'
+
output_string += '# | 0 | 0 | 1 | 3.4 |\n'
+
+
output_string += '\n# parameter x collects the displacements along x, '
+
output_string += 'normalized with the box length x.'
+
+
output_string += '\nparameter x:'
+
offset = 1
+
for i in range(0, self.n_control_points[0]):
+
for j in range(0, self.n_control_points[1]):
+
for k in range(0, self.n_control_points[2]):
+
output_string += offset * ' ' + str(i) + ' ' + str(
+
j) + ' ' + str(k) + ' ' + str(
+
self.array_mu_x[i][j][k]) + '\n'
+
offset = 13
+
+
output_string += '\n# parameter y collects the displacements along y, '
+
output_string += 'normalized with the box length y.'
+
+
output_string += '\nparameter y:'
+
offset = 1
+
for i in range(0, self.n_control_points[0]):
+
for j in range(0, self.n_control_points[1]):
+
for k in range(0, self.n_control_points[2]):
+
output_string += offset * ' ' + str(i) + ' ' + str(
+
j) + ' ' + str(k) + ' ' + str(
+
self.array_mu_y[i][j][k]) + '\n'
+
offset = 13
+
+
output_string += '\n# parameter z collects the displacements along z, '
+
output_string += 'normalized with the box length z.'
+
+
output_string += '\nparameter z:'
+
offset = 1
+
for i in range(0, self.n_control_points[0]):
+
for j in range(0, self.n_control_points[1]):
+
for k in range(0, self.n_control_points[2]):
+
output_string += offset * ' ' + str(i) + ' ' + str(
+
j) + ' ' + str(k) + ' ' + str(
+
self.array_mu_z[i][j][k]) + '\n'
+
offset = 13
+
+
with open(filename, 'w') as f:
+
f.write(output_string)
+
+
+
def __str__(self):
+
"""
+
This method prints all the FFD parameters on the screen. Its purpose is
+
for debugging.
+
"""
+
string = ""
+
string += 'conversion_unit = {}\n'.format(self.conversion_unit)
+
string += 'n_control_points = {}\n\n'.format(self.n_control_points)
+
string += 'box_length = {}\n'.format(self.box_length)
+
string += 'box_origin = {}\n'.format(self.box_origin)
+
string += 'rot_angle = {}\n'.format(self.rot_angle)
+
string += '\narray_mu_x =\n{}\n'.format(self.array_mu_x)
+
string += '\narray_mu_y =\n{}\n'.format(self.array_mu_y)
+
string += '\narray_mu_z =\n{}\n'.format(self.array_mu_z)
+
string += '\nrotation_matrix = \n{}\n'.format(self.rotation_matrix)
+
string += '\nposition_vertices = {}\n'.format(self.position_vertices)
+
return string
+
+
+
[docs]
+
def control_points(self, deformed=True):
+
"""
+
Method that returns the FFD control points. If the `deformed` flag is
+
set to True the method returns the deformed lattice, otherwise it
+
returns the original undeformed lattice.
+
+
:param bool deformed: flag to select the original or modified FFD
+
control lattice. The default is True.
+
:return: the FFD control points (by row).
+
:rtype: numpy.ndarray
+
"""
+
x = np.linspace(0, self.box_length[0], self.n_control_points[0])
+
y = np.linspace(0, self.box_length[1], self.n_control_points[1])
+
z = np.linspace(0, self.box_length[2], self.n_control_points[2])
+
+
y_coords, x_coords, z_coords = np.meshgrid(y, x, z)
+
+
box_points = np.array(
+
[x_coords.ravel(),
+
y_coords.ravel(),
+
z_coords.ravel()])
+
+
if deformed:
+
box_points += np.array([
+
self.array_mu_x.ravel() * self.box_length[0],
+
self.array_mu_y.ravel() * self.box_length[1],
+
self.array_mu_z.ravel() * self.box_length[2]
+
])
+
+
n_rows = box_points.shape[1]
+
+
box_points = np.dot(self.rotation_matrix, box_points) + np.transpose(
+
np.tile(self.box_origin, (n_rows, 1)))
+
+
return box_points.T
+
+
+
+
[docs]
+
def reflect(self, axis=0, in_place=True):
+
"""
+
Reflect the lattice of control points along the direction defined
+
by `axis`. In particular the origin point of the lattice is preserved.
+
So, for instance, the reflection along x, is made with respect to the
+
face of the lattice in the yz plane that is opposite to the origin.
+
Same for the other directions. Only the weights (mu) along the chosen
+
axis are reflected, while the others are preserved. The symmetry plane
+
can not present deformations along the chosen axis.
+
After the refletcion there will be 2n-1 control points along `axis`,
+
witha doubled box length.
+
+
:param int axis: axis along which the reflection is performed.
+
Default is 0. Possible values are 0, 1, or 2, corresponding
+
to x, y, and z respectively.
+
:param bool in_place: if True, the object attributes are modified in
+
place; if False, a new object is return with the reflected lattice.
+
Default is True.
+
:return: a new object with the same parameters and the reflected
+
lattice if `in_place` is False, otherwise NoneType.
+
+
:Example:
+
+
>>> ffd.reflect(axis=0, in_place=True) # irreversible
+
>>> # or ...
+
>>> refle_ffd = ffd.reflect(axis=0, in_place=False)
+
"""
+
# check axis value
+
if axis not in (0, 1, 2):
+
raise ValueError(
+
"The axis has to be 0, 1, or 2. Current value {}.".format(axis))
+
+
# check that the plane of symmetry is undeformed
+
if (axis == 0 and np.count_nonzero(self.array_mu_x[-1, :, :]) != 0) or (
+
axis == 1 and np.count_nonzero(self.array_mu_y[:, -1, :]) != 0
+
) or (axis == 2 and np.count_nonzero(self.array_mu_z[:, :, -1]) != 0):
+
raise RuntimeError(
+
"If you want to reflect the FFD bounding box along axis " + \
+
"{} you can not diplace the control ".format(axis) + \
+
"points in the symmetry plane along that axis."
+
)
+
+
if in_place is False:
+
self = copy.deepcopy(self)
+
+
# double the control points in the given axis -1 (the symmetry plane)
+
self.n_control_points[axis] = 2 * self.n_control_points[axis] - 1
+
# double the box length
+
self.box_length[axis] *= 2
+
+
# we have to reflect the dispacements only along the correct axis
+
reflection = np.ones(3)
+
reflection[axis] = -1
+
+
# we select all the indeces but the ones in the plane of symmetry
+
indeces = [slice(None), slice(None), slice(None)] # = [:, :, :]
+
indeces[axis] = slice(1, None) # = [1:]
+
indeces = tuple(indeces)
+
+
# we append along the given axis all the displacements reflected
+
# and in the reverse order
+
self.array_mu_x = np.append(self.array_mu_x,
+
reflection[0] *
+
np.flip(self.array_mu_x, axis)[indeces],
+
axis=axis)
+
self.array_mu_y = np.append(self.array_mu_y,
+
reflection[1] *
+
np.flip(self.array_mu_y, axis)[indeces],
+
axis=axis)
+
self.array_mu_z = np.append(self.array_mu_z,
+
reflection[2] *
+
np.flip(self.array_mu_z, axis)[indeces],
+
axis=axis)
+
if in_place is False:
+
return self
+
+
+
+
+
+
[docs]
+
def __call__(self, src_pts):
+
"""
+
This method performs the FFD to `src_pts` and return the deformed
+
points.
+
+
:param numpy.ndarray src_pts: the array of dimensions (*n_points*, *3*)
+
containing the points to deform. The points have to be arranged by
+
row.
+
:return: the deformed points
+
:rtype: numpy.ndarray (with shape = (*n_points*, *3*))
+
"""
+
def is_inside(pts, boundaries):
+
"""
+
Check is `pts` is inside the ranges provided by `boundaries`.
+
"""
+
return np.all(np.logical_and(pts >= boundaries[0],
+
pts <= boundaries[1]),
+
axis=1)
+
+
# map to the reference domain
+
src_reference_frame_pts = self.psi(src_pts - self.box_origin)
+
+
# apply deformation for all the pts in the unit cube
+
index_pts_inside = is_inside(src_reference_frame_pts,
+
np.array([[0., 0., 0.], [1., 1., 1.]]))
+
shifted_reference_frame_pts = self.T(
+
src_reference_frame_pts[index_pts_inside])
+
+
# map to the physical domain
+
shifted_pts = self.inverse_psi(
+
shifted_reference_frame_pts) + self.box_origin
+
+
dst_pts = src_pts.copy()
+
dst_pts[index_pts_inside] = shifted_pts
+
return dst_pts
+
- 
- 
- 
- 
- 
-
-
-
+