Skip to content

Symmetry #161

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Jan 15, 2019
Merged

Symmetry #161

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
469de0d
reflection method for FFD lattice
mtezzele Jan 14, 2019
4df6b55
test for FFD lattice reflection
mtezzele Jan 14, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 82 additions & 26 deletions pygem/params/ffdparams.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ class FFDParameters(object):
Four vertex (non coplanar) are sufficient to uniquely identify a
parallelepiped.
If the four vertex are coplanar, an assert is thrown when
`affine_points_fit is used.
`affine_points_fit` is used.

"""

Expand Down Expand Up @@ -123,10 +123,67 @@ def position_vertices(self):
:rtype: numpy.ndarray
"""
return self.origin_box + np.vstack([
np.zeros((1, 3)),
self.rotation_matrix.dot(np.diag(self.lenght_box)).T
np.zeros(
(1, 3)), self.rotation_matrix.dot(np.diag(self.lenght_box)).T
])

def reflect(self, axis=0):
"""
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.
"""
# 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."
)

# 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.lenght_box[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)

def read_parameters(self, filename='parameters.prm'):
"""
Reads in the parameters file and fill the self structure.
Expand Down Expand Up @@ -203,12 +260,12 @@ def write_parameters(self, filename='parameters.prm'):
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 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 lenght indicates the length of the FFD '
output_string += 'bounding box along the three canonical directions '
Expand Down Expand Up @@ -282,8 +339,8 @@ def write_parameters(self, filename='parameters.prm'):
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'
j) + ' ' + str(k) + ' ' + str(self.array_mu_x[i][j][
k]) + '\n'
offset = 13

output_string += '\n# parameter y collects the displacements along y, '
Expand All @@ -295,8 +352,8 @@ def write_parameters(self, filename='parameters.prm'):
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'
j) + ' ' + str(k) + ' ' + str(self.array_mu_y[i][j][
k]) + '\n'
offset = 13

output_string += '\n# parameter z collects the displacements along z, '
Expand All @@ -308,8 +365,8 @@ def write_parameters(self, filename='parameters.prm'):
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'
j) + ' ' + str(k) + ' ' + str(self.array_mu_z[i][j][
k]) + '\n'
offset = 13

with open(filename, 'w') as f:
Expand Down Expand Up @@ -364,29 +421,28 @@ def save_points(self, filename, write_deformed=True):
y = np.linspace(0, self.lenght_box[1], self.n_control_points[1])
z = np.linspace(0, self.lenght_box[2], self.n_control_points[2])

lattice_y_coords, lattice_x_coords, lattice_z_coords = np.meshgrid(
y, x, z)
lattice_y_coords, lattice_x_coords, lattice_z_coords = np.meshgrid(y, x,
z)

if write_deformed:
box_points = np.array([
lattice_x_coords.ravel() +
self.array_mu_x.ravel() * self.lenght_box[0],
lattice_y_coords.ravel() +
lattice_x_coords.ravel() + self.array_mu_x.ravel() *
self.lenght_box[0], lattice_y_coords.ravel() +
self.array_mu_y.ravel() * self.lenght_box[1],
lattice_z_coords.ravel() +
self.array_mu_z.ravel() * self.lenght_box[2]
lattice_z_coords.ravel() + self.array_mu_z.ravel() *
self.lenght_box[2]
])
else:
box_points = np.array([
lattice_x_coords.ravel(),
lattice_y_coords.ravel(),
lattice_x_coords.ravel(), lattice_y_coords.ravel(),
lattice_z_coords.ravel()
])

n_rows = box_points.shape[1]

box_points = np.dot(self.rotation_matrix, box_points) + np.transpose(
np.tile(self.origin_box, (n_rows, 1)))
box_points = np.dot(
self.rotation_matrix,
box_points) + np.transpose(np.tile(self.origin_box, (n_rows, 1)))

points = vtk.vtkPoints()

Expand Down
100 changes: 99 additions & 1 deletion tests/test_ffdparams.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,105 @@ def test_class_members_generic_array_mu_z(self):
np.testing.assert_array_almost_equal(params.array_mu_z,
np.zeros((2, 3, 5)))

def test_reflect_n_control_points_1(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=0)
assert np.array_equal(params.n_control_points, [3, 3, 5])

def test_reflect_n_control_points_2(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=1)
assert np.array_equal(params.n_control_points, [2, 5, 5])

def test_reflect_n_control_points_3(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=2)
assert np.array_equal(params.n_control_points, [2, 3, 9])

def test_reflect_box_length_1(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=0)
assert params.lenght_box[0] == 2

def test_reflect_box_length_2(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=1)
assert params.lenght_box[1] == 2

def test_reflect_box_length_3(self):
params = FFDParameters([2, 3, 5])
params.reflect(axis=2)
assert params.lenght_box[2] == 2

def test_reflect_wrong_axis(self):
params = FFDParameters([2, 3, 5])
with self.assertRaises(ValueError):
params.reflect(axis=4)

def test_reflect_wrong_symmetry_plane_1(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_x = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.]).reshape((3, 2,
2))
with self.assertRaises(RuntimeError):
params.reflect(axis=0)

def test_reflect_wrong_symmetry_plane_2(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_y = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.]).reshape((3, 2,
2))
with self.assertRaises(RuntimeError):
params.reflect(axis=1)

def test_reflect_wrong_symmetry_plane_3(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_z = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.1]).reshape((3, 2,
2))
with self.assertRaises(RuntimeError):
params.reflect(axis=2)

def test_reflect_axis_0(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_x = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., .2, 0., 0., 0., 0.]).reshape((3, 2,
2))
params.reflect(axis=0)
array_mu_x_exact = np.array([0.2, 0., 0., 0., 0.5, 0., 0., 0.2, 0.,
0., 0., 0., -0.5, -0., -0., -0.2, -0.2, -0., -0., -0.]).reshape((5, 2,
2))
np.testing.assert_array_almost_equal(params.array_mu_x,
array_mu_x_exact)

def test_reflect_axis_1(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_y = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 0., 0., 0., 0.]).reshape((3, 2,
2))
params.reflect(axis=1)
array_mu_y_exact = np.array([0.2, 0., 0., 0., -0.2, -0., 0.5, 0., 0., 0.,
-0.5, -0., 0., 0., 0., 0., 0., 0.]).reshape((3, 3, 2))
np.testing.assert_array_almost_equal(params.array_mu_y,
array_mu_y_exact)

def test_reflect_axis_2(self):
params = FFDParameters([3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
params.array_mu_z = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 0., 0., 0., 0.]).reshape((3, 2,
2))
params.reflect(axis=2)
array_mu_z_exact = np.array([0.2, 0., -0.2, 0., 0., 0., 0.5, 0., -0.5,
0., 0., -0., 0., 0., -0., 0., 0., -0.]).reshape((3, 2, 3))
np.testing.assert_array_almost_equal(params.array_mu_z,
array_mu_z_exact)

def test_read_parameters_conversion_unit(self):
params = FFDParameters(n_control_points=[3, 2, 2])
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
Expand Down Expand Up @@ -119,7 +218,6 @@ def test_read_parameters_array_mu_x(self):
array_mu_x_exact = np.array(
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0., 0.]).reshape((3, 2,
2))
print(params.array_mu_x)
np.testing.assert_array_almost_equal(params.array_mu_x,
array_mu_x_exact)

Expand Down