1313 :math:`\\ boldsymbol{\\ psi}`. In the code it is named *transformation*.
1414
1515 - Moving some control points to deform the lattice with :math:`\\ hat{T}`.
16- The movement of the control points is basically the weight (or displacement)
17- :math:`\\ boldsymbol{\\ mu}` we set in the *parameters file*.
16+ The movement of the control points is basically the weight (or
17+ displacement) :math:`\\ boldsymbol{\\ mu}` we set in the *parameters file*.
1818
1919 - Mapping back to the physical domain with map
2020 :math:`\\ boldsymbol{\\ psi}^{-1}`. In the code it is named
3535 where :math:`\\ mathsf{b}_{lmn}` are Bernstein polynomials. We improve the
3636 traditional version by allowing a rotation of the FFD lattice in order to
3737 give more flexibility to the tool.
38-
38+
3939 You can try to add more shapes to the lattice to allow more and more
4040 involved transformations.
4141
@@ -81,15 +81,14 @@ class FFD(Deformation):
8181
8282 :Example:
8383
84- >>> import pygem.freeform as ffd
85- >>> import pygem.params as ffdp
84+ >>> from pygem import FFD
8685 >>> import numpy as np
87- >>> ffd_params = ffdp.FFDParameters ()
88- >>> ffd_params .read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
89- >>> original_mesh_points = np.load( 'tests/test_datasets/meshpoints_sphere_orig.npy ')
90- >>> free_form = ffd.FFD(ffd_params, original_mesh_points)
91- >>> free_form.perform( )
92- >>> new_mesh_points = free_form.modified_mesh_points
86+ >>> ffd = FFD ()
87+ >>> ffd .read_parameters(
88+ >>> 'tests/test_datasets/parameters_test_ffd_sphere.prm ')
89+ >>> original_mesh_points = np.load(
90+ >>> 'tests/test_datasets/meshpoints_sphere_orig.npy' )
91+ >>> new_mesh_points = ffd(original_mesh_points)
9392 """
9493 reference_frame = np .array ([[0 , 0 , 0 ], [1 , 0 , 0 ], [0 , 1 , 0 ], [0 , 0 , 1 ]])
9594
@@ -100,6 +99,10 @@ def __init__(self, n_control_points=None):
10099 self .box_origin = np .array ([0. , 0. , 0. ])
101100 self .rot_angle = np .array ([0. , 0. , 0. ])
102101
102+ self .array_mu_x = None
103+ self .array_mu_y = None
104+ self .array_mu_z = None
105+
103106 if n_control_points is None :
104107 n_control_points = [2 , 2 , 2 ]
105108 self .n_control_points = n_control_points
@@ -120,11 +123,11 @@ def n_control_points(self, npts):
120123 self .array_mu_y = np .zeros (self .n_control_points )
121124 self .array_mu_z = np .zeros (self .n_control_points )
122125
123-
124126 @property
125127 def psi (self ):
126128 """
127- Return the function that map the physical domain to the reference domain.
129+ Return the function that map the physical domain to the reference
130+ domain.
128131
129132 :rtype: callable
130133 """
@@ -134,7 +137,8 @@ def psi(self):
134137 @property
135138 def inverse_psi (self ):
136139 """
137- Return the function that map the reference domain to the physical domain.
140+ Return the function that map the reference domain to the physical
141+ domain.
138142
139143 :rtype: callable
140144 """
@@ -144,11 +148,10 @@ def inverse_psi(self):
144148 @property
145149 def T (self ):
146150 """
147- Return the function that deforms the points within the unit cube, using the
151+ Return the function that deforms the points within the unit cube.
148152
149153 :rtype: callable
150154 """
151-
152155 def T_mapping (points ):
153156 (n_rows , n_cols ) = points .shape
154157 (dim_n_mu , dim_m_mu , dim_t_mu ) = self .array_mu_x .shape
@@ -162,35 +165,36 @@ def T_mapping(points):
162165
163166 # TODO check no-loop implementation
164167 #bernstein_x = (
165- # np.power(mesh_points[:, 0][:, None], range(dim_n_mu)) *
168+ # np.power(mesh_points[:, 0][:, None], range(dim_n_mu)) *
166169 # np.power(1 - mesh_points[:, 0][:, None], range(dim_n_mu-1, -1, -1)) *
167170 # special.binom(np.array([dim_n_mu-1]*dim_n_mu), np.arange(dim_n_mu))
168171 #)
169172 for i in range (0 , dim_n_mu ):
170173 aux1 = np .power ((1 - points [:, 0 ]), dim_n_mu - 1 - i )
171174 aux2 = np .power (points [:, 0 ], i )
172- bernstein_x [i , :] = (
173- special . binom ( dim_n_mu - 1 , i ) * np .multiply (aux1 , aux2 ))
175+ bernstein_x [i , :] = (special . binom ( dim_n_mu - 1 , i ) *
176+ np .multiply (aux1 , aux2 ))
174177
175178 for i in range (0 , dim_m_mu ):
176179 aux1 = np .power ((1 - points [:, 1 ]), dim_m_mu - 1 - i )
177180 aux2 = np .power (points [:, 1 ], i )
178- bernstein_y [i , :] = special .binom (dim_m_mu - 1 , i ) * np . multiply (
179- aux1 , aux2 )
181+ bernstein_y [i , :] = special .binom (dim_m_mu - 1 ,
182+ i ) * np . multiply ( aux1 , aux2 )
180183
181184 for i in range (0 , dim_t_mu ):
182185 aux1 = np .power ((1 - points [:, 2 ]), dim_t_mu - 1 - i )
183186 aux2 = np .power (points [:, 2 ], i )
184- bernstein_z [i , :] = special .binom (dim_t_mu - 1 , i ) * np . multiply (
185- aux1 , aux2 )
187+ bernstein_z [i , :] = special .binom (dim_t_mu - 1 ,
188+ i ) * np . multiply ( aux1 , aux2 )
186189
187190 aux_x = 0.
188191 aux_y = 0.
189192 aux_z = 0.
190193
191194 for j in range (0 , dim_m_mu ):
192195 for k in range (0 , dim_t_mu ):
193- bernstein_yz = np .multiply (bernstein_y [j , :], bernstein_z [k , :])
196+ bernstein_yz = np .multiply (bernstein_y [j , :],
197+ bernstein_z [k , :])
194198 for i in range (0 , dim_n_mu ):
195199 aux = np .multiply (bernstein_x [i , :], bernstein_yz )
196200 aux_x += aux * self .array_mu_x [i , j , k ]
@@ -201,8 +205,8 @@ def T_mapping(points):
201205 shift_points [1 , :] += aux_y
202206 shift_points [2 , :] += aux_z
203207 return shift_points .T + points
204- return T_mapping
205208
209+ return T_mapping
206210
207211 @property
208212 def rotation_matrix (self ):
@@ -212,9 +216,9 @@ def rotation_matrix(self):
212216
213217 :rtype: numpy.ndarray
214218 """
215- return angles2matrix (
216- np . radians ( self . rot_angle [ 2 ]), np .radians (self .rot_angle [1 ]),
217- np .radians (self .rot_angle [0 ]))
219+ return angles2matrix (np . radians ( self . rot_angle [ 2 ]),
220+ np .radians (self .rot_angle [1 ]),
221+ np .radians (self .rot_angle [0 ]))
218222
219223 @property
220224 def position_vertices (self ):
@@ -224,8 +228,8 @@ def position_vertices(self):
224228 :rtype: numpy.ndarray
225229 """
226230 return self .box_origin + np .vstack ([
227- np .zeros (
228- ( 1 , 3 )), self .rotation_matrix .dot (np .diag (self .box_length )).T
231+ np .zeros (( 1 , 3 )),
232+ self .rotation_matrix .dot (np .diag (self .box_length )).T
229233 ])
230234
231235 def reset_weights (self ):
@@ -312,12 +316,12 @@ def write_parameters(self, filename='parameters.prm'):
312316 output_string += ' points in each direction (x, y, z).\n '
313317 output_string += '# For example, to create a 2 x 3 x 2 grid, use the'
314318 output_string += ' following: n control points: 2, 3, 2\n '
315- output_string += 'n control points x: ' + str (self . n_control_points [
316- 0 ]) + '\n '
317- output_string += 'n control points y: ' + str (self . n_control_points [
318- 1 ]) + '\n '
319- output_string += 'n control points z: ' + str (self . n_control_points [
320- 2 ]) + '\n '
319+ output_string += 'n control points x: ' + str (
320+ self . n_control_points [ 0 ]) + '\n '
321+ output_string += 'n control points y: ' + str (
322+ self . n_control_points [ 1 ]) + '\n '
323+ output_string += 'n control points z: ' + str (
324+ self . n_control_points [ 2 ]) + '\n '
321325
322326 output_string += '\n # box length indicates the length of the FFD '
323327 output_string += 'bounding box along the three canonical directions '
@@ -391,8 +395,8 @@ def write_parameters(self, filename='parameters.prm'):
391395 for j in range (0 , self .n_control_points [1 ]):
392396 for k in range (0 , self .n_control_points [2 ]):
393397 output_string += offset * ' ' + str (i ) + ' ' + str (
394- j ) + ' ' + str (k ) + ' ' + str (self . array_mu_x [ i ][ j ][
395- k ]) + '\n '
398+ j ) + ' ' + str (k ) + ' ' + str (
399+ self . array_mu_x [ i ][ j ][ k ]) + '\n '
396400 offset = 13
397401
398402 output_string += '\n # parameter y collects the displacements along y, '
@@ -404,8 +408,8 @@ def write_parameters(self, filename='parameters.prm'):
404408 for j in range (0 , self .n_control_points [1 ]):
405409 for k in range (0 , self .n_control_points [2 ]):
406410 output_string += offset * ' ' + str (i ) + ' ' + str (
407- j ) + ' ' + str (k ) + ' ' + str (self . array_mu_y [ i ][ j ][
408- k ]) + '\n '
411+ j ) + ' ' + str (k ) + ' ' + str (
412+ self . array_mu_y [ i ][ j ][ k ]) + '\n '
409413 offset = 13
410414
411415 output_string += '\n # parameter z collects the displacements along z, '
@@ -417,8 +421,8 @@ def write_parameters(self, filename='parameters.prm'):
417421 for j in range (0 , self .n_control_points [1 ]):
418422 for k in range (0 , self .n_control_points [2 ]):
419423 output_string += offset * ' ' + str (i ) + ' ' + str (
420- j ) + ' ' + str (k ) + ' ' + str (self . array_mu_z [ i ][ j ][
421- k ]) + '\n '
424+ j ) + ' ' + str (k ) + ' ' + str (
425+ self . array_mu_z [ i ][ j ][ k ]) + '\n '
422426 offset = 13
423427
424428 with open (filename , 'w' ) as f :
@@ -442,7 +446,6 @@ def __str__(self):
442446 string += '\n position_vertices = {}\n ' .format (self .position_vertices )
443447 return string
444448
445-
446449 def control_points (self , deformed = True ):
447450 """
448451 Method that returns the FFD control points. If the `deformed` flag is
@@ -460,8 +463,10 @@ def control_points(self, deformed=True):
460463
461464 y_coords , x_coords , z_coords = np .meshgrid (y , x , z )
462465
463- box_points = np .array ([
464- x_coords .ravel (), y_coords .ravel (), z_coords .ravel ()])
466+ box_points = np .array (
467+ [x_coords .ravel (),
468+ y_coords .ravel (),
469+ z_coords .ravel ()])
465470
466471 if deformed :
467472 box_points += np .array ([
@@ -472,13 +477,11 @@ def control_points(self, deformed=True):
472477
473478 n_rows = box_points .shape [1 ]
474479
475- box_points = np .dot (
476- self .rotation_matrix ,
477- box_points ) + np .transpose (np .tile (self .box_origin , (n_rows , 1 )))
480+ box_points = np .dot (self .rotation_matrix , box_points ) + np .transpose (
481+ np .tile (self .box_origin , (n_rows , 1 )))
478482
479483 return box_points .T
480484
481-
482485 def reflect (self , axis = 0 ):
483486 """
484487 Reflect the lattice of control points along the direction defined
@@ -526,37 +529,44 @@ def reflect(self, axis=0):
526529
527530 # we append along the given axis all the displacements reflected
528531 # and in the reverse order
529- self .array_mu_x = np .append (
530- self .array_mu_x ,
531- reflection [0 ] * np .flip (self .array_mu_x , axis )[indeces ], axis = axis )
532- self .array_mu_y = np .append (
533- self .array_mu_y ,
534- reflection [1 ] * np .flip (self .array_mu_y , axis )[indeces ], axis = axis )
535- self .array_mu_z = np .append (
536- self .array_mu_z ,
537- reflection [2 ] * np .flip (self .array_mu_z , axis )[indeces ], axis = axis )
538-
539-
532+ self .array_mu_x = np .append (self .array_mu_x ,
533+ reflection [0 ] *
534+ np .flip (self .array_mu_x , axis )[indeces ],
535+ axis = axis )
536+ self .array_mu_y = np .append (self .array_mu_y ,
537+ reflection [1 ] *
538+ np .flip (self .array_mu_y , axis )[indeces ],
539+ axis = axis )
540+ self .array_mu_z = np .append (self .array_mu_z ,
541+ reflection [2 ] *
542+ np .flip (self .array_mu_z , axis )[indeces ],
543+ axis = axis )
540544
541545 def __call__ (self , src_pts ):
542546 """
543547 This method performs the deformation on the mesh pts. After the
544548 execution it sets `self.modified_mesh_pts`.
545549 """
546550 def is_inside (pts , boundaries ):
547- return np .all (np .logical_and (
548- pts >= boundaries [0 ], pts <= boundaries [1 ]), axis = 1 )
551+ """
552+ Check is `pts` is inside the ranges provided by `boundaries`.
553+ """
554+ return np .all (np .logical_and (pts >= boundaries [0 ],
555+ pts <= boundaries [1 ]),
556+ axis = 1 )
549557
550558 # map to the reference domain
551559 src_reference_frame_pts = self .psi (src_pts - self .box_origin )
552560
553561 # apply deformation for all the pts in the unit cube
554- index_pts_inside = is_inside (
555- src_reference_frame_pts , np .array ([[0. , 0. , 0. ], [1. , 1. , 1. ]]))
556- shifted_reference_frame_pts = self .T (src_reference_frame_pts [index_pts_inside ])
562+ index_pts_inside = is_inside (src_reference_frame_pts ,
563+ np .array ([[0. , 0. , 0. ], [1. , 1. , 1. ]]))
564+ shifted_reference_frame_pts = self .T (
565+ src_reference_frame_pts [index_pts_inside ])
557566
558567 # map to the physical domain
559- shifted_pts = self .inverse_psi (shifted_reference_frame_pts ) + self .box_origin
568+ shifted_pts = self .inverse_psi (
569+ shifted_reference_frame_pts ) + self .box_origin
560570
561571 dst_pts = src_pts .copy ()
562572 dst_pts [index_pts_inside ] = shifted_pts
0 commit comments