[REBASED ] Symmetric Diffeomorphic registration workflow

bin/dipy_align_syn

@@ -0,0 +1,7 @@
+#!python
+
+from dipy.workflows.flow_runner import run_flow
+from dipy.workflows.align import SynRegistrationFlow
+
+if __name__ == "__main__":
+    run_flow(SynRegistrationFlow())

bin/dipy_syn_register

@@ -0,0 +1,6 @@
+#!python
+from dipy.workflows.flow_runner import run_flow
+from dipy.workflows.align import SynRegistrationFlow
+
+if __name__ == "__main__":
+    run_flow(SynRegistrationFlow())

dipy/io/image.py

@@ -40,7 +40,6 @@ def save_qa_metric(fname, xopt, fopt):
         image registration.
     fopt: int
         The distance between the registered images.
-
     """
     np.savetxt(fname, xopt, header="Optimal Parameter metric")
     with open(fname, 'a') as f:

dipy/workflows/align.py

@@ -6,6 +6,8 @@
 
 from dipy.align.imaffine import AffineMap, transform_centers_of_mass, \
     MutualInformationMetric, AffineRegistration
+from dipy.align.imwarp import SymmetricDiffeomorphicRegistration
+from dipy.align.metrics import CCMetric, SSDMetric, EMMetric
 from dipy.align.reslice import reslice
 from dipy.align.transforms import (TranslationTransform3D, RigidTransform3D,
                                    AffineTransform3D)
@@ -676,7 +678,6 @@ def run(self, static_image_file, moving_image_files, affine_matrix_file,
 
             # Loading the image data from the input files into object.
             static_image, static_grid2world = load_nifti(static_image_file)
-
             moving_image, moving_grid2world = load_nifti(moving_image_file)
 
             # Doing a sanity check for validating the dimensions of the input
@@ -698,3 +699,188 @@ def run(self, static_image_file, moving_image_files, affine_matrix_file,
             transformed = img_transformation.transform(moving_image)
 
             save_nifti(out_file, transformed, affine=static_grid2world)
+
+
+class SynRegistrationFlow(Workflow):
+
+    def run(self, static_image_file, moving_image_file, affine_matrix_file,
+            inv_static=False, level_iters=[10, 10, 5], metric="cc",
+            mopt_sigma_diff=2.0, mopt_radius=4, mopt_smooth=0.0,
+            mopt_inner_iter=0.0, mopt_q_levels=256, mopt_double_gradient=True,
+            mopt_step_type='', step_length=0.25,
+            ss_sigma_factor=0.2, opt_tol=1e-5, inv_iter=20,
+            inv_tol=1e-3, out_dir='', out_warped='warped_moved.nii.gz',
+            out_inv_static='inc_static.nii.gz',
+            out_field='displacefield.txt'):
+
+        """
+        Parameters
+        ----------
+        static_image_file : string
+            Path of the static image file.
+
+        moving_image_file : string
+            Path to the moving image file.
+
+        affine_matrix_file : string
+            The text file containing pre alignment information or the
+             affine matrix.
+
+        inv_static : boolean, optional
+            Apply the inverse mapping to the static image (default 'False').
+
+        level_iters : variable int, optional
+            The number of iterations at each level of the gaussian pyramid.
+             By default, a 3-level scale space with iterations
+             sequence equal to [10, 10, 5] will be used. The 0-th
+             level corresponds to the finest resolution.
+
+        metric : string, optional
+            The metric to be used (Default cc, 'Cross Correlation metric').
+            metric available: cc (Cross Correlation), ssd (Sum Squared
+            Difference), em (Expectation-Maximization).
+
+        mopt_sigma_diff : float, optional
+            Metric option applied on Cross correlation (CC).
+            The standard deviation of the Gaussian smoothing kernel to be
+            applied to the update field at each iteration (default 2.0)
+
+        mopt_radius : int, optional
+            Metric option applied on Cross correlation (CC).
+            the radius of the squared (cubic) neighborhood at each voxel to
+            be considered to compute the cross correlation. (default 4)
+
+        mopt_smooth : float, optional
+            Metric option applied on Sum Squared Difference (SSD) and
+            Expectation Maximization (EM). Smoothness parameter, the
+            larger the value the smoother the deformation field.
+            (default 1.0 for EM, 4.0 for SSD)
+
+        mopt_inner_iter : int, optional
+            Metric option applied on Sum Squared Difference (SSD) and
+            Expectation Maximization (EM). This is number of iterations to be
+            performed at each level of the multi-resolution Gauss-Seidel
+            optimization algorithm (this is not the number of steps per
+            Gaussian Pyramid level, that parameter must be set for the
+            optimizer, not the metric). Default 5 for EM, 10 for SSD.
+
+        mopt_q_levels : int, optional
+            Metric option applied on Expectation Maximization (EM).
+            Number of quantization levels (Default: 256 for EM)
+
+        mopt_double_gradient : bool, optional
+            Metric option applied on Expectation Maximization (EM).
+            if True, the gradient of the expected static image under the moving
+            modality will be added to the gradient of the moving image,
+            similarly, the gradient of the expected moving image under the
+            static modality will be added to the gradient of the static image.
+
+        mopt_step_type : string, optional
+            Metric option applied on Sum Squared Difference (SSD) and
+            Expectation Maximization (EM). The optimization schedule to be
+            used in the multi-resolution Gauss-Seidel optimization algorithm
+            (not used if Demons Step is selected). Possible value:
+            ('gauss_newton', 'demons'). default: 'gauss_newton' for EM,
+            'demons' for SSD.
+
+        step_length : float, optional
+            the length of the maximum displacement vector of the update
+             displacement field at each iteration.
+
+        ss_sigma_factor : float, optional
+            parameter of the scale-space smoothing kernel. For example, the
+             std. dev. of the kernel will be factor*(2^i) in the isotropic case
+             where i = 0, 1, ..., n_scales is the scale.
+
+        opt_tol : float, optional
+            the optimization will stop when the estimated derivative of the
+             energy profile w.r.t. time falls below this threshold.
+
+        inv_iter : int, optional
+            the number of iterations to be performed by the displacement field
+             inversion algorithm.
+
+        inv_tol : float, optional
+            the displacement field inversion algorithm will stop iterating
+             when the inversion error falls below this threshold.
+
+        out_dir : string, optional
+            Directory to save the transformed files (default '').
+
+        out_warped : string, optional
+            Name of the warped file. If no name is given then a
+             suffix 'transformed' will be appended to the name of the
+             original input file (default 'warped_moved.nii.gz').
+
+        out_inv_static : string, optional
+            Name of the file to save the static image after applying the
+             inverse mapping (default 'inv_static.nii.gz').
+
+        out_field : string, optional
+            Name of the file to save the diffeomorphic field.
+
+        """
+        io_it = self.get_io_iterator()
+        metric = metric.lower()
+        if metric not in ['ssd', 'cc', 'em']:
+            raise ValueError("Invalid similarity metric: Please"
+                             " provide a valid metric.")
+
+        # Init parameter if they are not setup
+        init_param = {'ssd': {'mopt_smooth': 4.0,
+                              'mopt_inner_iter': 10,
+                              'mopt_step_type': 'demons'
+                              },
+                      'em': {'mopt_smooth': 1.0,
+                             'mopt_inner_iter': 5,
+                             'mopt_step_type': 'gauss_newton'
+                             }
+                      }
+        mopt_smooth = mopt_smooth or init_param[metric]['mopt_smooth']
+        mopt_inner_iter = mopt_inner_iter or  \
+            init_param[metric]['mopt_inner_iter']
+        mopt_step_type = mopt_step_type or \
+            init_param[metric]['mopt_step_type']
+
+        for static_file, moving_file, in_affine, \
+                warped_file, inv_static_file, displ_file in io_it:
+
+            # Loading the image data from the input files into object.
+            static_image, static_grid2world = load_nifti(static_file)
+            moving_image, moving_grid2world = load_nifti(moving_file)
+
+            # Sanity check for the input image dimensions.
+            ImageRegistrationFlow.check_dimensions(static_image, moving_image)
+
+            # Loading the affine matrix.
+            affine_matrix = np.loadtxt(in_affine)
+
+            l_metric = {"ssd": SSDMetric(static_image.ndim,
+                                         smooth=mopt_smooth,
+                                         inner_iter=mopt_inner_iter,
+                                         step_type=mopt_step_type
+                                         ),
+                        "cc": CCMetric(static_image.ndim,
+                                       sigma_diff=mopt_sigma_diff,
+                                       radius=mopt_radius),
+                        "em": EMMetric(static_image.ndim,
+                                       smooth=mopt_smooth,
+                                       inner_iter=mopt_inner_iter,
+                                       step_type=mopt_step_type,
+                                       q_levels=mopt_q_levels,
+                                       double_gradient=mopt_double_gradient)
+                        }
+
+            current_metric = l_metric.get(metric.lower())
+
+            sdr = SymmetricDiffeomorphicRegistration(current_metric,
+                                                     level_iters)
+
+            mapping = sdr.optimize(static_image, moving_image,
+                                   static_grid2world, moving_grid2world,
+                                   affine_matrix)
+
+            warped_moving = mapping.transform(moving_image)
+
+            # Saving the warped moving file and the alignment matrix.
+            save_nifti(warped_file, warped_moving, static_grid2world)

dipy/workflows/tests/test_align.py

@@ -1,7 +1,6 @@
 import numpy.testing as npt
 import numpy as np
 
-import numpy as np
 import nibabel as nib
 from nibabel.tmpdirs import TemporaryDirectory
 from dipy.tracking.streamline import Streamlines
@@ -14,7 +13,7 @@
 from dipy.align.tests.test_parzenhist import setup_random_transform
 from dipy.align.transforms import regtransforms
 from dipy.io.image import save_nifti
-from dipy.workflows.align import ImageRegistrationFlow
+from dipy.workflows.align import ImageRegistrationFlow, SynRegistrationFlow
 from dipy.workflows.align import ApplyAffineFlow