Skip to content

Commit f7405da

Browse files
mtezzelemtezzele
authored
Merge pull request #161 from mtezzele/symmetry
Symmetry
2 parents 3676a3a + 4df6b55 commit f7405da

File tree

2 files changed

+181
-27
lines changed

2 files changed

+181
-27
lines changed

‎pygem/params/ffdparams.py‎

Lines changed: 82 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ class FFDParameters(object):
6666
Four vertex (non coplanar) are sufficient to uniquely identify a
6767
parallelepiped.
6868
If the four vertex are coplanar, an assert is thrown when
69-
`affine_points_fit is used.
69+
`affine_points_fit` is used.
7070
7171
"""
7272

@@ -123,10 +123,67 @@ def position_vertices(self):
123123
:rtype: numpy.ndarray
124124
"""
125125
return self.origin_box + np.vstack([
126-
np.zeros((1, 3)),
127-
self.rotation_matrix.dot(np.diag(self.lenght_box)).T
126+
np.zeros(
127+
(1, 3)), self.rotation_matrix.dot(np.diag(self.lenght_box)).T
128128
])
129129

130+
def reflect(self, axis=0):
131+
"""
132+
Reflect the lattice of control points along the direction defined
133+
by `axis`. In particular the origin point of the lattice is preserved.
134+
So, for instance, the reflection along x, is made with respect to the
135+
face of the lattice in the yz plane that is opposite to the origin.
136+
Same for the other directions. Only the weights (mu) along the chosen
137+
axis are reflected, while the others are preserved. The symmetry plane
138+
can not present deformations along the chosen axis.
139+
After the refletcion there will be 2n-1 control points along `axis`,
140+
witha doubled box length.
141+
142+
:param int axis: axis along which the reflection is performed.
143+
Default is 0. Possible values are 0, 1, or 2, corresponding
144+
to x, y, and z respectively.
145+
"""
146+
# check axis value
147+
if axis not in (0, 1, 2):
148+
raise ValueError(
149+
"The axis has to be 0, 1, or 2. Current value {}.".format(axis))
150+
151+
# check that the plane of symmetry is undeformed
152+
if (axis == 0 and np.count_nonzero(self.array_mu_x[-1, :, :]) != 0) or (
153+
axis == 1 and np.count_nonzero(self.array_mu_y[:, -1, :]) != 0
154+
) or (axis == 2 and np.count_nonzero(self.array_mu_z[:, :, -1]) != 0):
155+
raise RuntimeError(
156+
"If you want to reflect the FFD bounding box along axis " + \
157+
"{} you can not diplace the control ".format(axis) + \
158+
"points in the symmetry plane along that axis."
159+
)
160+
161+
# double the control points in the given axis -1 (the symmetry plane)
162+
self.n_control_points[axis] = 2 * self.n_control_points[axis] - 1
163+
# double the box length
164+
self.lenght_box[axis] *= 2
165+
166+
# we have to reflect the dispacements only along the correct axis
167+
reflection = np.ones(3)
168+
reflection[axis] = -1
169+
170+
# we select all the indeces but the ones in the plane of symmetry
171+
indeces = [slice(None), slice(None), slice(None)] # = [:, :, :]
172+
indeces[axis] = slice(1, None) # = [1:]
173+
indeces = tuple(indeces)
174+
175+
# we append along the given axis all the displacements reflected
176+
# and in the reverse order
177+
self.array_mu_x = np.append(
178+
self.array_mu_x,
179+
reflection[0] * np.flip(self.array_mu_x, axis)[indeces], axis=axis)
180+
self.array_mu_y = np.append(
181+
self.array_mu_y,
182+
reflection[1] * np.flip(self.array_mu_y, axis)[indeces], axis=axis)
183+
self.array_mu_z = np.append(
184+
self.array_mu_z,
185+
reflection[2] * np.flip(self.array_mu_z, axis)[indeces], axis=axis)
186+
130187
def read_parameters(self, filename='parameters.prm'):
131188
"""
132189
Reads in the parameters file and fill the self structure.
@@ -203,12 +260,12 @@ def write_parameters(self, filename='parameters.prm'):
203260
output_string += ' points in each direction (x, y, z).\n'
204261
output_string += '# For example, to create a 2 x 3 x 2 grid, use the'
205262
output_string += ' following: n control points: 2, 3, 2\n'
206-
output_string += 'n control points x: ' + str(
207-
self.n_control_points[0]) + '\n'
208-
output_string += 'n control points y: ' + str(
209-
self.n_control_points[1]) + '\n'
210-
output_string += 'n control points z: ' + str(
211-
self.n_control_points[2]) + '\n'
263+
output_string += 'n control points x: ' + str(self.n_control_points[
264+
0]) + '\n'
265+
output_string += 'n control points y: ' + str(self.n_control_points[
266+
1]) + '\n'
267+
output_string += 'n control points z: ' + str(self.n_control_points[
268+
2]) + '\n'
212269

213270
output_string += '\n# box lenght indicates the length of the FFD '
214271
output_string += 'bounding box along the three canonical directions '
@@ -282,8 +339,8 @@ def write_parameters(self, filename='parameters.prm'):
282339
for j in range(0, self.n_control_points[1]):
283340
for k in range(0, self.n_control_points[2]):
284341
output_string += offset * ' ' + str(i) + ' ' + str(
285-
j) + ' ' + str(k) + ' ' + str(
286-
self.array_mu_x[i][j][k]) + '\n'
342+
j) + ' ' + str(k) + ' ' + str(self.array_mu_x[i][j][
343+
k]) + '\n'
287344
offset = 13
288345

289346
output_string += '\n# parameter y collects the displacements along y, '
@@ -295,8 +352,8 @@ def write_parameters(self, filename='parameters.prm'):
295352
for j in range(0, self.n_control_points[1]):
296353
for k in range(0, self.n_control_points[2]):
297354
output_string += offset * ' ' + str(i) + ' ' + str(
298-
j) + ' ' + str(k) + ' ' + str(
299-
self.array_mu_y[i][j][k]) + '\n'
355+
j) + ' ' + str(k) + ' ' + str(self.array_mu_y[i][j][
356+
k]) + '\n'
300357
offset = 13
301358

302359
output_string += '\n# parameter z collects the displacements along z, '
@@ -308,8 +365,8 @@ def write_parameters(self, filename='parameters.prm'):
308365
for j in range(0, self.n_control_points[1]):
309366
for k in range(0, self.n_control_points[2]):
310367
output_string += offset * ' ' + str(i) + ' ' + str(
311-
j) + ' ' + str(k) + ' ' + str(
312-
self.array_mu_z[i][j][k]) + '\n'
368+
j) + ' ' + str(k) + ' ' + str(self.array_mu_z[i][j][
369+
k]) + '\n'
313370
offset = 13
314371

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

367-
lattice_y_coords, lattice_x_coords, lattice_z_coords = np.meshgrid(
368-
y, x, z)
424+
lattice_y_coords, lattice_x_coords, lattice_z_coords = np.meshgrid(y, x,
425+
z)
369426

370427
if write_deformed:
371428
box_points = np.array([
372-
lattice_x_coords.ravel() +
373-
self.array_mu_x.ravel() * self.lenght_box[0],
374-
lattice_y_coords.ravel() +
429+
lattice_x_coords.ravel() + self.array_mu_x.ravel() *
430+
self.lenght_box[0], lattice_y_coords.ravel() +
375431
self.array_mu_y.ravel() * self.lenght_box[1],
376-
lattice_z_coords.ravel() +
377-
self.array_mu_z.ravel() * self.lenght_box[2]
432+
lattice_z_coords.ravel() + self.array_mu_z.ravel() *
433+
self.lenght_box[2]
378434
])
379435
else:
380436
box_points = np.array([
381-
lattice_x_coords.ravel(),
382-
lattice_y_coords.ravel(),
437+
lattice_x_coords.ravel(), lattice_y_coords.ravel(),
383438
lattice_z_coords.ravel()
384439
])
385440

386441
n_rows = box_points.shape[1]
387442

388-
box_points = np.dot(self.rotation_matrix, box_points) + np.transpose(
389-
np.tile(self.origin_box, (n_rows, 1)))
443+
box_points = np.dot(
444+
self.rotation_matrix,
445+
box_points) + np.transpose(np.tile(self.origin_box, (n_rows, 1)))
390446

391447
points = vtk.vtkPoints()
392448

‎tests/test_ffdparams.py‎

Lines changed: 99 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,105 @@ def test_class_members_generic_array_mu_z(self):
8686
np.testing.assert_array_almost_equal(params.array_mu_z,
8787
np.zeros((2, 3, 5)))
8888

89+
def test_reflect_n_control_points_1(self):
90+
params = FFDParameters([2, 3, 5])
91+
params.reflect(axis=0)
92+
assert np.array_equal(params.n_control_points, [3, 3, 5])
93+
94+
def test_reflect_n_control_points_2(self):
95+
params = FFDParameters([2, 3, 5])
96+
params.reflect(axis=1)
97+
assert np.array_equal(params.n_control_points, [2, 5, 5])
98+
99+
def test_reflect_n_control_points_3(self):
100+
params = FFDParameters([2, 3, 5])
101+
params.reflect(axis=2)
102+
assert np.array_equal(params.n_control_points, [2, 3, 9])
103+
104+
def test_reflect_box_length_1(self):
105+
params = FFDParameters([2, 3, 5])
106+
params.reflect(axis=0)
107+
assert params.lenght_box[0] == 2
108+
109+
def test_reflect_box_length_2(self):
110+
params = FFDParameters([2, 3, 5])
111+
params.reflect(axis=1)
112+
assert params.lenght_box[1] == 2
113+
114+
def test_reflect_box_length_3(self):
115+
params = FFDParameters([2, 3, 5])
116+
params.reflect(axis=2)
117+
assert params.lenght_box[2] == 2
118+
119+
def test_reflect_wrong_axis(self):
120+
params = FFDParameters([2, 3, 5])
121+
with self.assertRaises(ValueError):
122+
params.reflect(axis=4)
123+
124+
def test_reflect_wrong_symmetry_plane_1(self):
125+
params = FFDParameters([3, 2, 2])
126+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
127+
params.array_mu_x = np.array(
128+
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.]).reshape((3, 2,
129+
2))
130+
with self.assertRaises(RuntimeError):
131+
params.reflect(axis=0)
132+
133+
def test_reflect_wrong_symmetry_plane_2(self):
134+
params = FFDParameters([3, 2, 2])
135+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
136+
params.array_mu_y = np.array(
137+
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.]).reshape((3, 2,
138+
2))
139+
with self.assertRaises(RuntimeError):
140+
params.reflect(axis=1)
141+
142+
def test_reflect_wrong_symmetry_plane_3(self):
143+
params = FFDParameters([3, 2, 2])
144+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
145+
params.array_mu_z = np.array(
146+
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0.3, 0.1]).reshape((3, 2,
147+
2))
148+
with self.assertRaises(RuntimeError):
149+
params.reflect(axis=2)
150+
151+
def test_reflect_axis_0(self):
152+
params = FFDParameters([3, 2, 2])
153+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
154+
params.array_mu_x = np.array(
155+
[0.2, 0., 0., 0., 0.5, 0., 0., .2, 0., 0., 0., 0.]).reshape((3, 2,
156+
2))
157+
params.reflect(axis=0)
158+
array_mu_x_exact = np.array([0.2, 0., 0., 0., 0.5, 0., 0., 0.2, 0.,
159+
0., 0., 0., -0.5, -0., -0., -0.2, -0.2, -0., -0., -0.]).reshape((5, 2,
160+
2))
161+
np.testing.assert_array_almost_equal(params.array_mu_x,
162+
array_mu_x_exact)
163+
164+
def test_reflect_axis_1(self):
165+
params = FFDParameters([3, 2, 2])
166+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
167+
params.array_mu_y = np.array(
168+
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 0., 0., 0., 0.]).reshape((3, 2,
169+
2))
170+
params.reflect(axis=1)
171+
array_mu_y_exact = np.array([0.2, 0., 0., 0., -0.2, -0., 0.5, 0., 0., 0.,
172+
-0.5, -0., 0., 0., 0., 0., 0., 0.]).reshape((3, 3, 2))
173+
np.testing.assert_array_almost_equal(params.array_mu_y,
174+
array_mu_y_exact)
175+
176+
def test_reflect_axis_2(self):
177+
params = FFDParameters([3, 2, 2])
178+
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
179+
params.array_mu_z = np.array(
180+
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 0., 0., 0., 0.]).reshape((3, 2,
181+
2))
182+
params.reflect(axis=2)
183+
array_mu_z_exact = np.array([0.2, 0., -0.2, 0., 0., 0., 0.5, 0., -0.5,
184+
0., 0., -0., 0., 0., -0., 0., 0., -0.]).reshape((3, 2, 3))
185+
np.testing.assert_array_almost_equal(params.array_mu_z,
186+
array_mu_z_exact)
187+
89188
def test_read_parameters_conversion_unit(self):
90189
params = FFDParameters(n_control_points=[3, 2, 2])
91190
params.read_parameters('tests/test_datasets/parameters_sphere.prm')
@@ -119,7 +218,6 @@ def test_read_parameters_array_mu_x(self):
119218
array_mu_x_exact = np.array(
120219
[0.2, 0., 0., 0., 0.5, 0., 0., 0., 1., 0., 0., 0.]).reshape((3, 2,
121220
2))
122-
print(params.array_mu_x)
123221
np.testing.assert_array_almost_equal(params.array_mu_x,
124222
array_mu_x_exact)
125223

0 commit comments

Comments
 (0)