completed transition from indexes to maskes

Author: guglielmopadula

pygem/cffd.py

@@ -60,20 +60,37 @@ class CFFD(FFD):
         >>> 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):
+    def __init__(self, n_control_points=None, fun=None, fixval=None, M=None, mask=None ):
         super().__init__(n_control_points)
-        self.fun = None
-        self.fixval = None
-        self.indices = None
-        self.M = None
+
+        if mask==None:
+            self.mask=np.full((*self.n_control_points,3), True, dtype=bool)
+        else:
+            self.mask=mask
+
+        if fixval==None:
+            self.fixval=np.array([1.])    
+        else:
+            self.fixval=fixval
+        
+        if fun==None:
+            self.fun=lambda x: self.fixval
+        
+        else:
+            self.fun=fun
+
+        if M==None:
+            self.M=np.eye(np.sum(self.mask.astype(int)))
 
     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.fixval - d)
-        saved_parameters[self.indices] = saved_parameters[self.indices] + deltax
+        indices=np.arange(np.prod(self.n_control_points)*3)[self.mask.reshape(-1)]
+        A, b = self._compute_linear_map(src_pts, saved_parameters.copy(),indices)
+        d = A @ saved_parameters[indices] + b
+        invM=np.linalg.inv(self.M)
+        deltax = invM @ A.T @ np.linalg.inv(
+            (A @ invM @ A.T)) @ (self.fixval - d)
+        saved_parameters[indices] = saved_parameters[indices] + deltax
         self._load_parameters(saved_parameters)
         return self.ffd(src_pts)
 
@@ -112,10 +129,14 @@ def _load_parameters(self, tmp):
         self.array_mu_y = tmp[:, :, :, 1]
         self.array_mu_z = tmp[:, :, :, 2]
 
+    def read_parameters(self,filename='parameters.prm'):
+        super().read_parameters(filename)
+        self.mask=np.full((*self.n_control_points,3), True, dtype=bool)
+        self.M=np.eye(np.sum(self.mask.astype(int)))
 
 # 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):
+    def _compute_linear_map(self, src_pts, saved_parameters,indices):
         '''
         Computes the coefficient and the intercept of the linear map from the control points to the output.
         
@@ -124,7 +145,7 @@ def _compute_linear_map(self, src_pts, saved_parameters):
         :return: a tuple containing the coefficient and the intercept.
         :rtype: tuple(np.ndarray,np.ndarray)
         '''
-        n_indices = len(self.indices)
+        n_indices = len(indices)
         inputs = np.zeros([n_indices + 1, n_indices + 1])
         outputs = np.zeros([n_indices + 1, self.fixval.shape[0]])
         np.random.seed(0)
@@ -134,7 +155,7 @@ def _compute_linear_map(self, src_pts, saved_parameters):
             tmp = tmp.reshape(1, -1)
             inputs[i] = np.hstack([tmp, np.ones(
                 (tmp.shape[0], 1))])  #dependent variable
-            saved_parameters[self.indices] = tmp
+            saved_parameters[indices] = tmp
             self._load_parameters(
                 saved_parameters
             )  #loading the depent variable as a control point
@@ -146,3 +167,22 @@ def _compute_linear_map(self, src_pts, saved_parameters):
         A = sol[0].T[:, :-1]  #coefficient
         b = sol[0].T[:, -1]  #intercept
         return A, b
+
+np.random.seed(0)
+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.fun = fun
+cffd.fixval = b
+new_mesh_points = cffd(original_mesh_points)
+assert np.isclose(np.linalg.norm(fun(new_mesh_points) - b),
+                    np.array([0.0]))

pygem/vffd.py

@@ -73,19 +73,15 @@ def volume(x):
 
     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]
+        mask_bak=self.mask.copy()        
         diffvolume = self.fixval - self.fun(self.ffd(src_pts))
         for i in range(3):
-            self.indices = indexes[i]
-            self.M = np.eye(len(self.indices))
+            self.mask=np.full((*self.n_control_points,3), False, dtype=bool)
+            self.mask[:,:,:,i]=mask_bak[:,:,:,i].copy()
+            self.M = np.eye(np.sum(self.mask.astype(int)))
             self.fixval = self.fun(
                 self.ffd(src_pts)) + self.vweight[i] * (diffvolume)
             _ = super().__call__(src_pts)
         tmp = super().__call__(src_pts)
-        self.indices = indices_bak
+        self.mask = mask_bak
         return tmp

tests/test_cffd.py

@@ -20,8 +20,6 @@ def fun(x):
         b = fun(original_mesh_points)
         cffd.fun = fun
         cffd.fixval = 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.array_equal(original_mesh_points, new_mesh_points)
 
@@ -41,8 +39,6 @@ def fun(x):
         b = fun(original_mesh_points)
         cffd.fun = fun
         cffd.fixval = 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.0]))
@@ -60,13 +56,12 @@ def fun(x):
         b = fun(original_mesh_points)
         cffd.fixval = b
         cffd.fun = fun
-        cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
-        cffd.M = np.eye(len(cffd.indices))
         save_par = cffd._save_parameters()
-        C, d = cffd._compute_linear_map(original_mesh_points, save_par.copy())
-        for i in range(2 * len(cffd.indices)):
-            tmp = np.random.rand(len(cffd.indices))
-            save_par[cffd.indices] = tmp
+        indices=np.arange(np.prod(cffd.n_control_points)*3)[cffd.mask.reshape(-1)]
+        C, d = cffd._compute_linear_map(original_mesh_points, save_par.copy(),indices)
+        for i in range(2 * len(indices)):
+            tmp = np.random.rand(len(indices))
+            save_par[indices] = tmp
             cffd._load_parameters(tmp)
             assert np.allclose(
                 np.linalg.norm(C @ tmp + d -

tests/test_vffd.py

@@ -38,8 +38,6 @@ def test_nothing_happens(self):
         cffd.vweight = np.array([1 / 3, 1 / 3, 1 / 3])
         b = cffd.fun(points)
         cffd.fixval = np.array([b])
-        cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
-        cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(points)
         assert np.allclose(np.linalg.norm(points - new_mesh_points), 0.0)
 
@@ -76,8 +74,6 @@ def test_constraint(self):
         cffd.read_parameters(
             "tests/test_datasets/parameters_test_ffd_sphere.prm")
         cffd.fixval = np.array([b])
-        cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
-        cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(points)
         assert np.isclose(np.linalg.norm(cffd.fun(new_mesh_points) - b),
                           np.array([0.0]))