added cffd and variations

Author: guglielmopadula

pygem/__init__.py

@@ -8,10 +8,16 @@
     "idw",
     "rbf_factory",
     "custom_deformation",
+    "cffd"
+    "bffd"
+    "vffd"
 ]
 
 from .deformation import Deformation
 from .ffd import FFD
+from .cffd import CFFD
+from .bffd import BFFD
+from .vffd import VFFD
 from .rbf import RBF
 from .idw import IDW
 from .rbf_factory import RBFFactory

pygem/bffd.py

@@ -0,0 +1,58 @@
+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.
+
+ 
+"""
+
+
+
+class BFFD(CFFD):
+    '''
+    Class that handles the Barycenter Free Form Deformation on the mesh points.
+ 
+    :Example:
+
+        >>> from pygem import BFFD
+        >>> import numpy as np
+        >>> bffd = BFFD()
+        >>> bffd.read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
+        >>> original_mesh_points = np.load('tests/test_datasets/meshpoints_sphere_orig.npy')
+        >>> b=bffd.linconstraint(original_mesh_points)
+        >>> bffd.valconstraint=b
+        >>> bffd.indices=np.arange(np.prod(bffd.n_control_points)*3).tolist()
+        >>> 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)
+        def linfun(x):
+            return np.mean(x.reshape(-1,3),axis=0)
+
+        self.linconstraint=linfun
+
+    def __call__(self, src_pts):
+        return super().__call__(src_pts)
+
+if __name__ == "__main__":
+        from pygem import BFFD
+        import numpy as np
+        bffd = BFFD()
+        bffd.read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
+        original_mesh_points = np.load('tests/test_datasets/meshpoints_sphere_orig.npy')
+        b=bffd.linconstraint(original_mesh_points)
+        bffd.valconstraint=b
+        bffd.indices=np.arange(np.prod(bffd.n_control_points)*3).tolist()
+        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.]))

pygem/cffd.py

@@ -0,0 +1,136 @@
+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.
+
+ 
+"""
+
+
+class CFFD(FFD):
+    """
+    Class that handles the Constrained Free Form Deformation on the mesh points.
+
+    :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 
+        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
+        >>> import numpy as np
+        >>> cffd = CFFD()
+        >>> cffd.read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
+        >>> original_mesh_points = np.load('tests/test_datasets/meshpoints_sphere_orig.npy')
+        >>> A=np.random.rand(3,original_mesh_points.reshape(-1).shape[0])
+        >>> def fun(x):
+        >>>     x=x.reshape(-1)
+        >>>     return A@x
+        >>> b=fun(original_mesh_points)
+        >>> cffd.linconstraint=fun
+        >>> cffd.valconstraint=b
+        >>> cffd.indices=np.arange(np.prod(cffd.n_control_points)*3).tolist()
+        >>> 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
+
+    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
+        self._load_parameters(saved_parameters)
+        return self.ffd(src_pts)
+
+    def ffd(self,src_pts):
+        '''
+            Performs Classic Free Form Deformation.
+    
+            :param np.ndarray src_pts
+            :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
+        '''
+        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
+        '''
+        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
+        
+        :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)
+
+        '''
+        
+        
+        saved_parameters_bak=saved_parameters.copy()
+        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]
+            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)
+        return A,b

pygem/vffd.py

@@ -0,0 +1,67 @@
+from pygem.cffd import CFFD
+import numpy as np
+from copy import deepcopy
+
+"""
+Utilities for performing Volume Free Form Deformation (BFFD).
+
+:Theoretical Insight:
+    It performs Free Form Deformation while trying to enforce the volume 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.
+
+"""
+
+
+class VFFD(CFFD):
+    '''
+    Class that handles the Volumetric Free Form Deformation on the mesh points.
+ 
+    cvar np.ndarray vweight: specifies the weight of every step of the volume enforcement.
+
+    :Example:
+
+        >>> from pygem import VFFD
+        >>> import numpy as np
+        >>> import meshio
+        >>> mesh = meshio.read('tests/test_datasets/test_sphere.stl')  
+        >>> original_mesh_points=mesh.points
+        >>> triangles = mesh.cells_dict["triangle"]
+        >>> vffd = VFFD(triangles,[2,2,2])
+        >>> vffd.read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
+        >>> b=vffd.linconstraint(original_mesh_points)
+        >>> vffd.valconstraint=np.array([b])
+        >>> vffd.indices=np.arange(np.prod(vffd.n_control_points)*3).tolist()
+        >>> vffd.M=np.eye(len(vffd.indices))
+        >>> new_mesh_points = vffd(original_mesh_points)
+        >>> assert np.isclose(np.linalg.norm(vffd.linconstraint(new_mesh_points)-b),np.array([0.]))
+
+    '''
+    def __init__(self,triangles ,n_control_points=None):
+        super().__init__(n_control_points)
+        self.triangles=triangles
+        self.vweight=[1/3,1/3,1/3]
+        def volume(x):
+            x=x.reshape(-1,3)
+            mesh=x[self.triangles]
+            return np.sum(np.linalg.det(mesh)) 
+        self.linconstraint=volume
+
+
+    def __call__(self,src_pts):
+        self.vweight=np.abs(self.vweight)/np.sum(np.abs(self.vweight))
+        indices_bak=deepcopy(self.indices)
+        self.indices=np.array(self.indices)
+        indices_x=self.indices[self.indices%3==0].tolist()
+        indices_y=self.indices[self.indices%3==1].tolist()
+        indices_z=self.indices[self.indices%3==2].tolist()
+        indexes=[indices_x,indices_y,indices_z]
+        diffvolume=self.valconstraint-self.linconstraint(self.ffd(src_pts))
+        for i in range(3):
+            self.indices=indexes[i]
+            self.M=np.eye(len(self.indices))
+            self.valconstraint=self.linconstraint(self.ffd(src_pts))+self.vweight[i]*(diffvolume)
+            _=super().__call__(src_pts)
+        tmp=super().__call__(src_pts)
+        self.indices=indices_bak
+        return tmp