reformat code, fixed sphinx, add test, simplied interpolation function with comments

Author: guglielmopadula

docs/source/_tutorials/tutorial-8-cffd.html

@@ -0,0 +1,18 @@
+Barycenter Free-form Deformation
+=====================
+
+.. currentmodule:: pygem.bffd
+
+.. automodule:: pygem.bffd
+
+.. autosummary::
+	:toctree: _summaries
+	:nosignatures:
+
+.. autoclass:: pygem.bffd.BFFD
+	:members:
+	:special-members: __call__
+	:private-members:
+	:exclude-members: _abd_impl
+	:show-inheritance:
+	:noindex:

docs/source/bffd.rst

@@ -0,0 +1,18 @@
+Constrained Free-form Deformation
+=====================
+
+.. currentmodule:: pygem.cffd
+
+.. automodule:: pygem.cffd
+
+.. autosummary::
+	:toctree: _summaries
+	:nosignatures:
+
+.. autoclass:: pygem.cffd.CFFD
+	:members:
+	:special-members: __call__
+	:private-members:
+	:exclude-members: _abd_impl
+	:show-inheritance:
+	:noindex:

docs/source/cffd.rst

@@ -5,6 +5,9 @@ Code Documentation
    :maxdepth: 2
 
    ffd
+   cffd
+   bffd
+   vffd   
    rbf
    rbf_factory
    idw

docs/source/code.rst

@@ -4,9 +4,9 @@ Tutorials
 We made some tutorial examples:
 
 - `Tutorial 1 <tutorial-1-ffd.html>`_ shows how to apply the free-form deformation to mesh nodes.
-- `Turorial 2 <tutorial-2-iges.html>`_ shows how to deal with iges (CAD) files and FFD.
-- `Turorial 3 <tutorial-3-rbf.html>`_ shows how to apply the radial basis function to mesh nodes.
-- `Turorial 4 <tutorial-4-idw.html>`_ shows how to perform a deformation using the inverse distance weighting method to mesh nodes.
-- `Turorial 5 <tutorial-5-file.html>`_ shows how to perform a deformation to an object stored in a file.
-- `Turorial 6 <tutorial-6-ffd-rbf.html>`_ shows deforming a computational mesh.
-
+- `Tutorial 2 <tutorial-2-iges.html>`_ shows how to deal with iges (CAD) files and FFD.
+- `Tutorial 3 <tutorial-3-rbf.html>`_ shows how to apply the radial basis function to mesh nodes.
+- `Tutorial 4 <tutorial-4-idw.html>`_ shows how to perform a deformation using the inverse distance weighting method to mesh nodes.
+- `Tutorial 5 <tutorial-5-file.html>`_ shows how to perform a deformation to an object stored in a file.
+- `Tutorial 6 <tutorial-6-ffd-rbf.html>`_ shows deforming a computational mesh.
+- `Tutorial 8 <tutorial-8-cffd.html>`_ shows deforming a computational mesh.

docs/source/tutorials.rst

@@ -0,0 +1,18 @@
+Volume Free-form Deformation
+=====================
+
+.. currentmodule:: pygem.vffd
+
+.. automodule:: pygem.vffd
+
+.. autosummary::
+	:toctree: _summaries
+	:nosignatures:
+
+.. autoclass:: pygem.vffd.VFFD
+	:members:
+	:special-members: __call__
+	:private-members:
+	:exclude-members: _abd_impl
+	:show-inheritance:
+	:noindex:

docs/source/vffd.rst

@@ -1,23 +1,45 @@
-from pygem.cffd import CFFD
-import numpy as np
-
 """
 Utilities for performing Barycenter Free Form Deformation (BFFD).
 
 :Theoretical Insight:
     It performs Free Form Deformation while trying to enforce the barycenter to be a certain value specified by the user.
     The constraint is enforced exactly (up to numerical errors).
-    For details see the mother and the grandmother clases.
+    For details see the mother and the grandmother classes.
 
- 
 """
 
+from pygem.cffd import CFFD
+import numpy as np
 
 
 class BFFD(CFFD):
     '''
     Class that handles the Barycenter 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.
+    :cvar callable linconstraint: it defines the F of the constraint F(x)=c.
+    :cvar numpy.ndarray valconstraint: it defines the c of the constraint F(x)=c.
+    :cvar list indices: it defines the indices of the control points 
+        that are moved to enforce the constraint. The control index is obtained by doing:
+        all_indices=np.arange(n_x*n_y*n_z*3).reshape(n_x,n_y,n_z,3).tolist().
+    :cvar numpy.ndarray M: a SDP weigth matrix. It must be of size len(indices) x len(indices).
+
     :Example:
 
         >>> from pygem import BFFD
@@ -31,8 +53,6 @@ class BFFD(CFFD):
         >>> bffd.M=np.eye(len(bffd.indices))
         >>> new_mesh_points = bffd(original_mesh_points)
         >>> assert np.isclose(np.linalg.norm(bffd.linconstraint(new_mesh_points)-b),np.array([0.]))
-
-    
     '''
     def __init__(self, n_control_points=None):
         super().__init__(n_control_points)

pygem/bffd.py

@@ -1,31 +1,46 @@
-from pygem.ffd import FFD
-import numpy as np
-
 """
 Utilities for performing Constrained Free Form Deformation (CFFD).
 
 :Theoretical Insight:
+    
     It performs Free Form Deformation while trying to enforce a costraint of the form F(x)=c. 
     The constraint is enforced exactly (up to numerical errors) if and only if the function is linear.
     For details on Free Form Deformation see the mother class.
 
- 
 """
 
+from pygem.ffd import FFD
+import numpy as np
+
 
 class CFFD(FFD):
     """
     Class that handles the Constrained 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.
     :cvar callable linconstraint: it defines the F of the constraint F(x)=c.
-    
     :cvar numpy.ndarray valconstraint: it defines the c of the constraint F(x)=c.
-    :param list indices: it defines the indices of the control points 
+    :cvar list indices: it defines the indices of the control points 
         that are moved to enforce the constraint. The control index is obtained by doing:
         all_indices=np.arange(n_x*n_y*n_z*3).reshape(n_x,n_y,n_z,3).tolist().
     :cvar numpy.ndarray M: a SDP weigth matrix. It must be of size len(indices) x len(indices).
 
-    
     :Example:
 
         >>> from pygem import CFFD
@@ -44,93 +59,86 @@ class CFFD(FFD):
         >>> cffd.M=np.eye(len(cffd.indices))
         >>> new_mesh_points = cffd(original_mesh_points)
         >>> assert np.isclose(np.linalg.norm(fun(new_mesh_points)-b),np.array([0.]))
-
     """
+
     def __init__(self, n_control_points=None):
         super().__init__(n_control_points)
-        self.linconstraint=None 
-        self.valconstraint=None
-        self.indices=None
-        self.M=None
+        self.linconstraint = None 
+        self.valconstraint = None
+        self.indices = None
+        self.M = None
 
     def __call__(self, src_pts):
         saved_parameters=self._save_parameters()
-        A,b=self._compute_linear_map(src_pts,saved_parameters.copy())
-        d=A@saved_parameters[self.indices]+b
-        deltax=np.linalg.inv(self.M)@A.T@np.linalg.inv((A@np.linalg.inv(self.M)@A.T))@(self.valconstraint-d)
-        saved_parameters[self.indices]=saved_parameters[self.indices]+deltax
+        A, b = self._compute_linear_map(src_pts, saved_parameters.copy())
+        d = A @ saved_parameters[self.indices] + b
+        deltax = np.linalg.inv(self.M) @ A.T @ np.linalg.inv((A @ np.linalg.inv(self.M)@ A.T)) @ (self.valconstraint - d)
+        saved_parameters[self.indices] = saved_parameters[self.indices] + deltax
         self._load_parameters(saved_parameters)
         return self.ffd(src_pts)
 
     def ffd(self,src_pts):
         '''
-            Performs Classic Free Form Deformation.
+        Performs Classic Free Form Deformation.
     
-            :param np.ndarray src_pts
-            :return the deformed points
-            :rtype: numpy.ndarray
+        :param np.ndarray src_pts: the points to deform.
+        :return: the deformed points.
+        :rtype: numpy.ndarray
         '''
         return super().__call__(src_pts)
 
 
     def _save_parameters(self):
         '''
-            Saves the FFD control points in an array of shape [n_x,ny,nz,3]
-            :return the FFD control points in an array of shape [n_x,ny,nz,3]
-            :rtype: numpy.ndarray
+        Saves the FFD control points in an array of shape [n_x,ny,nz,3].
+
+        :return: the FFD control points in an array of shape [n_x,ny,nz,3].
+        :rtype: numpy.ndarray
         '''
-        tmp=np.zeros([*self.n_control_points,3])
-        tmp[:,:,:,0]=self.array_mu_x
-        tmp[:,:,:,1]=self.array_mu_y
-        tmp[:,:,:,2]=self.array_mu_z
+        tmp = np.zeros([*self.n_control_points,3])
+        tmp[:, :, :, 0] = self.array_mu_x
+        tmp[:, :, :, 1] = self.array_mu_y
+        tmp[:, :, :, 2] = self.array_mu_z
         return tmp.reshape(-1)
 
     def _load_parameters(self,tmp):
         '''
-            Loads the FFD control points from an array of shape [n_x,ny,nz,3]
-            :param np.ndarray tmp
-            :rtype: None
+        Loads the FFD control points from an array of shape [n_x,ny,nz,3].
+
+        :param np.ndarray tmp: the array of FFD control points.
+        :rtype: None
         '''
-        tmp=tmp.reshape(*self.n_control_points,3)
-        self.array_mu_x=tmp[:,:,:,0]
-        self.array_mu_y=tmp[:,:,:,1]
-        self.array_mu_z=tmp[:,:,:,2]
+        tmp = tmp.reshape(*self.n_control_points,3)
+        self.array_mu_x = tmp[:, :, :, 0]
+        self.array_mu_y = tmp[:, :, :, 1]
+        self.array_mu_z = tmp[:, :, :, 2]
 
 
 # I see that a similar function already exists in pygem.utils, but it does not work for inputs and outputs of different dimensions
     def _compute_linear_map(self,src_pts,saved_parameters):
         '''
-        Computes the coefficient and the intercept of the linear map from the control points to the output
+        Computes the coefficient and the intercept of the linear map from the control points to the output.
         
-        :param np.ndarray src_pts
-        :param np.ndarray saved_parameters
-        :return a tuple containing the coefficient and the intercept
-        :rtype tuple(np.ndarray,np.ndarray)
-
+        :param np.ndarray src_pts: the points to deform.
+        :param np.ndarray saved_parameters: the array of FFD control points.
+        :return: a tuple containing the coefficient and the intercept.
+        :rtype: tuple(np.ndarray,np.ndarray)
         '''
-        
-        
-        saved_parameters_bak=saved_parameters.copy()
+        saved_parameters_bak=saved_parameters.copy() #saving ffd parameters
         n_indices=len(self.indices)
         inputs=np.zeros([n_indices+1,n_indices+1])
         outputs=np.zeros([n_indices+1,self.valconstraint.shape[0]])
-        tmp=saved_parameters_bak[self.indices].reshape(1,-1)
-        inputs[0]=np.hstack([tmp, np.ones((tmp.shape[0], 1))])
-        saved_parameters[self.indices]=tmp
-        self._load_parameters(saved_parameters)
-        def_pts=super().__call__(src_pts)
-        outputs[0]=self.linconstraint(def_pts)
-        for i in range(1,n_indices+1):
-            tmp=np.eye(n_indices)[i%n_indices]
+        np.random.seed(0)
+        for i in range(n_indices+1): ##now we generate the interpolation points
+            tmp=np.random.rand(1,n_indices)
             tmp=tmp.reshape(1,-1)
-            inputs[i]=np.hstack([tmp, np.ones((tmp.shape[0], 1))])
-            saved_parameters[self.indices]=tmp
-            self._load_parameters(saved_parameters)
-            def_pts=super().__call__(src_pts)
-            outputs[i]=self.linconstraint(def_pts)
-        self._load_parameters(saved_parameters_bak)
-        sol=np.linalg.lstsq(inputs,outputs,rcond=None)
-        A=sol[0].T[:,:-1]
-        b=sol[0].T[:,-1]
-        self._load_parameters(saved_parameters_bak)
+            inputs[i]=np.hstack([tmp, np.ones((tmp.shape[0], 1))]) #dependent variable
+            saved_parameters[self.indices]=tmp 
+            self._load_parameters(saved_parameters) #loading the depent variable as a control point
+            def_pts=super().__call__(src_pts) #computing the deformation with the dependent variable
+            outputs[i]=self.linconstraint(def_pts) #computing the independent variable
+        sol=np.linalg.lstsq(inputs,outputs,rcond=None) #computation of the linear map
+        A=sol[0].T[:,:-1] #coefficient
+        b=sol[0].T[:,-1] #intercept
+        self._load_parameters(saved_parameters_bak) #restoring the original FFD parameters
         return A,b