From 469de0df49fc2dff682377c1fa3ac00ef8af94e0 Mon Sep 17 00:00:00 2001 From: Marco Tezzele Date: Mon, 14 Jan 2019 11:38:34 +0100 Subject: [PATCH 1/2] reflection method for FFD lattice --- pygem/params/ffdparams.py | 60 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/pygem/params/ffdparams.py b/pygem/params/ffdparams.py index 82a7886..b074d5c 100644 --- a/pygem/params/ffdparams.py +++ b/pygem/params/ffdparams.py @@ -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. """ @@ -127,6 +127,64 @@ def position_vertices(self): 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.sum(self.array_mu_x[-1, :, :]) != 0) or \ + (axis == 1 and np.sum(self.array_mu_y[:, -1, :]) != 0) or \ + (axis == 2 and np.sum(self.array_mu_z[:, :, -1]) != 0): + raise RuntimeError("If you want to reflect the FFD " + \ + "bounding box along axis {} ".format(axis) + \ + "you can not diplace the control points in the " + \ + "symmetry plane along that axis.") + + # double the control points in the given axis minus 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 need to 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. From 4df6b55d697ec875a55fb757bd34c619d82b236f Mon Sep 17 00:00:00 2001 From: Marco Tezzele Date: Mon, 14 Jan 2019 14:08:03 +0100 Subject: [PATCH 2/2] test for FFD lattice reflection --- pygem/params/ffdparams.py | 102 +++++++++++++++++++------------------- tests/test_ffdparams.py | 100 ++++++++++++++++++++++++++++++++++++- 2 files changed, 149 insertions(+), 53 deletions(-) diff --git a/pygem/params/ffdparams.py b/pygem/params/ffdparams.py index b074d5c..463d682 100644 --- a/pygem/params/ffdparams.py +++ b/pygem/params/ffdparams.py @@ -123,11 +123,10 @@ 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 @@ -146,44 +145,44 @@ def reflect(self, axis=0): """ # 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)) + 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.sum(self.array_mu_x[-1, :, :]) != 0) or \ - (axis == 1 and np.sum(self.array_mu_y[:, -1, :]) != 0) or \ - (axis == 2 and np.sum(self.array_mu_z[:, :, -1]) != 0): - raise RuntimeError("If you want to reflect the FFD " + \ - "bounding box along axis {} ".format(axis) + \ - "you can not diplace the control points in the " + \ - "symmetry plane along that axis.") - - # double the control points in the given axis minus 1 (the symmetry plane) + 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 need to 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:] + # 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) - + 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'): """ @@ -261,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 ' @@ -340,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, ' @@ -353,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, ' @@ -366,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: @@ -422,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() diff --git a/tests/test_ffdparams.py b/tests/test_ffdparams.py index 601b524..750c35a 100644 --- a/tests/test_ffdparams.py +++ b/tests/test_ffdparams.py @@ -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') @@ -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)