--- a
+++ b/ants/registration/landmark_transforms.py
@@ -0,0 +1,843 @@
+__all__ = ["fit_transform_to_paired_points", 
+           "fit_time_varying_transform_to_point_sets"]
+
+import numpy as np
+import math
+import time
+
+import ants
+
+def convergence_monitoring(values, window_size=10):
+     if len(values) >= window_size:
+         u = np.linspace(0.0, 1.0, num=window_size)
+         scattered_data = np.expand_dims(values[-window_size:], axis=-1)
+         parametric_data = np.expand_dims(u, axis=-1)
+         spacing = 1 / (window_size-1)
+         bspline_line = ants.fit_bspline_object_to_scattered_data(scattered_data, parametric_data,
+             parametric_domain_origin=[0.0], parametric_domain_spacing=[spacing],
+             parametric_domain_size=[window_size], number_of_fitting_levels=1, mesh_size=1,
+             spline_order=1)
+         bspline_slope = -(bspline_line[1][0] - bspline_line[0][0]) / spacing
+         return(bspline_slope)
+     else:
+         return None
+
+
+def fit_transform_to_paired_points(moving_points,
+                                   fixed_points,
+                                   transform_type="affine",
+                                   regularization=1e-6,
+                                   domain_image=None,
+                                   number_of_fitting_levels=4,
+                                   mesh_size=1,
+                                   spline_order=3,
+                                   enforce_stationary_boundary=True,
+                                   displacement_weights=None,
+                                   number_of_compositions=10,
+                                   composition_step_size=0.5,
+                                   sigma=0.0,
+                                   convergence_threshold=1e-6,
+                                   number_of_time_steps=2,
+                                   number_of_integration_steps=100,
+                                   rasterize_points=False,
+                                   verbose=False
+                                  ):
+    """
+    Estimate a transform from corresponding fixed and moving landmarks.
+
+    ANTsR function: fitTransformToPairedPoints
+
+    Arguments
+    ---------
+    moving_points : array
+        Moving points specified in physical space as a n x d matrix where n is the number
+        of points and d is the dimensionality.
+
+    fixed_points : array
+        Fixed points specified in physical space as a n x d matrix where n is the number
+        of points and d is the dimensionality.
+
+    transform_type : character
+        'rigid', 'similarity', "affine', 'bspline', 'tps', 'diffeo', 'syn', or 'time-varying (tv)'.
+
+    regularization : scalar
+        Ridge penalty in [0,1] for linear transforms.
+
+    domain_image : ANTs image
+        Defines physical domain of the nonlinear transform.  Must be defined for nonlinear
+        transforms.
+
+    number_of_fitting_levels : integer
+        Integer specifying the number of fitting levels for the B-spline interpolation of the
+        displacement field.
+
+    mesh_size : integer or array
+        Defines the mesh size at the initial fitting level for the B-spline interpolation of the
+        displacement field.
+
+    spline_order : integer
+        Spline order of the B-spline displacement field.
+
+    enforce_stationary_boundary : boolean
+        Ensure no displacements on the image boundary (B-spline only).
+
+    displacement_weights : array
+        Defines the individual weighting of the corresponding scattered data value.  Default = NULL
+        meaning all displacements are weighted the same.
+
+    number_of_compositions : integer
+        Total number of compositions for the diffeomorphic transforms.
+
+    composition_step_size : scalar
+        Scalar multiplication factor of the weighting of the update field for the diffeomorphic transforms.
+
+    sigma : scalar
+        Gaussian smoothing standard deviation of the update field (in mm).
+
+    convergence_threshold : scalar
+        Composition-based convergence parameter for the diff. transforms using a
+        window size of 10 values.
+
+    number_of_time_steps : integer
+        Time-varying velocity field parameter.
+
+    number_of_integration_steps : scalar
+        Number of steps used for integrating the velocity field.
+
+    rasterize_points : boolean
+       Use nearest neighbor rasterization of points for estimating the update
+       field (potential speed-up).  Default = False.
+
+    verbose : bool
+        Print progress to the screen.
+
+    Returns
+    -------
+
+    ANTs transform
+
+    Example
+    -------
+    >>> import ants
+    >>> import numpy as np
+    >>> fixed = np.array([[50.0,50.0],[200.0,50.0],[200.0,200.0]])
+    >>> moving = np.array([[50.0,50.0],[50.0,200.0],[200.0,200.0]])
+    >>> xfrm = ants.fit_transform_to_paired_points(moving, fixed, transform_type="affine")
+    >>> xfrm = ants.fit_transform_to_paired_points(moving, fixed, transform_type="rigid")
+    >>> xfrm = ants.fit_transform_to_paired_points(moving, fixed, transform_type="similarity")
+    >>> domain_image = ants.image_read(ants.get_ants_data("r16"))
+    >>> xfrm = ants.fit_transform_to_paired_points(moving, fixed, transform_type="bspline", domain_image=domain_image, number_of_fitting_levels=5)
+    >>> xfrm = ants.fit_transform_to_paired_points(moving, fixed, transform_type="diffeo", domain_image=domain_image, number_of_fitting_levels=6)
+    """
+
+    def polar_decomposition(X):
+         U, d, V = np.linalg.svd(X, full_matrices=False)
+         P = np.matmul(U, np.matmul(np.diag(d), np.transpose(U)))
+         Z = np.matmul(U, V)
+         if np.linalg.det(Z) < 0:
+             n = X.shape[0]
+             reflection_matrix = np.identity(n)
+             reflection_matrix[0,0] = -1.0
+             Z = np.matmul(Z, reflection_matrix)
+         return({"P" : P, "Z" : Z, "Xtilde" : np.matmul(P, Z)})
+
+    def create_zero_displacement_field(domain_image):
+         field_array = np.zeros((*domain_image.shape, domain_image.dimension))
+         field = ants.from_numpy(field_array, origin=domain_image.origin,
+                 spacing=domain_image.spacing, direction=domain_image.direction,
+                 has_components=True)
+         return(field)
+
+    def create_zero_velocity_field(domain_image, number_of_time_points=2):
+         field_array = np.zeros((*domain_image.shape, number_of_time_points, domain_image.dimension))
+         origin = (*domain_image.origin, 0.0)
+         spacing = (*domain_image.spacing, 1.0)
+         direction = np.eye(domain_image.dimension + 1)
+         direction[0:domain_image.dimension,0:domain_image.dimension] = domain_image.direction
+         field = ants.from_numpy(field_array, origin=origin, spacing=spacing, direction=direction,
+                 has_components=True)
+         return(field)
+
+    allowed_transforms = ['rigid', 'affine', 'similarity', 'bspline', 'tps', 'diffeo', 'syn', 'tv', 'time-varying']
+    if not transform_type.lower() in allowed_transforms:
+        raise ValueError(transform_type + " transform not supported.")
+
+    transform_type = transform_type.lower()
+
+    if domain_image is None and transform_type in ['bspline', 'tps', 'diffeo', 'syn', 'tv', 'time-varying']:
+        raise ValueError("Domain image needs to be specified.")
+
+    if not fixed_points.shape == moving_points.shape:
+        raise ValueError("Mismatch in the size of the point sets.")
+
+    if regularization > 1:
+        regularization = 1
+    elif regularization < 0:
+        regularization = 0
+
+    number_of_points = fixed_points.shape[0]
+    dimensionality = fixed_points.shape[1]
+
+    if transform_type in ['rigid', 'affine', 'similarity']:
+        center_fixed = fixed_points.mean(axis=0)
+        center_moving = moving_points.mean(axis=0)
+
+        x = fixed_points - center_fixed
+        y = moving_points - center_moving
+
+        y_prior = np.concatenate((y, np.ones((number_of_points, 1))), axis=1)
+
+        x11 = np.concatenate((x, np.ones((number_of_points, 1))), axis=1)
+        M = x11 * (1.0 - regularization) + regularization * y_prior
+        Minv = np.linalg.lstsq(M, y, rcond=None)[0]
+
+        p = polar_decomposition(Minv[0:dimensionality, 0:dimensionality].T)
+        A = p['Xtilde']
+        translation = Minv[dimensionality,:] + center_moving - center_fixed
+
+        if transform_type in ['rigid', 'similarity']:
+            # Kabsch algorithm
+            #    http://web.stanford.edu/class/cs273/refs/umeyama.pdf
+
+            C = np.dot(y.T, x)
+            x_svd = np.linalg.svd(C * (1.0 - regularization) + np.eye(dimensionality) * regularization)
+            x_det = np.linalg.det(np.dot(x_svd[0], x_svd[2]))
+
+            if x_det < 0:
+                x_svd[2][dimensionality-1, :] *= -1
+
+            A = np.dot(x_svd[0], x_svd[2])
+
+            if transform_type == 'similarity':
+                scaling = (math.sqrt((np.power(y, 2).sum(axis=1) / number_of_points).mean()) /
+                           math.sqrt((np.power(x, 2).sum(axis=1) / number_of_points).mean()))
+                A = np.dot(A, np.eye(dimensionality) * scaling)
+
+        xfrm = ants.create_ants_transform(matrix=A, translation=translation,
+              dimension=dimensionality, center=center_fixed)
+
+        return xfrm
+
+    elif transform_type == "bspline":
+
+        bspline_displacement_field = ants.fit_bspline_displacement_field(
+            displacement_origins=fixed_points,
+            displacements=moving_points - fixed_points,
+            displacement_weights=displacement_weights,
+            origin=domain_image.origin,
+            spacing=domain_image.spacing,
+            size=domain_image.shape,
+            direction=domain_image.direction,
+            number_of_fitting_levels=number_of_fitting_levels,
+            mesh_size=mesh_size,
+            spline_order=spline_order,
+            enforce_stationary_boundary=enforce_stationary_boundary,
+            rasterize_points=rasterize_points)
+
+        xfrm = ants.transform_from_displacement_field(bspline_displacement_field)
+
+        return xfrm
+
+    elif transform_type == "tps":
+
+        tps_displacement_field = ants.fit_thin_plate_spline_displacement_field(
+            displacement_origins=fixed_points,
+            displacements=moving_points - fixed_points,
+            origin=domain_image.origin,
+            spacing=domain_image.spacing,
+            size=domain_image.shape,
+            direction=domain_image.direction)
+
+        xfrm = ants.transform_from_displacement_field(tps_displacement_field)
+
+        return xfrm
+
+    elif transform_type == "diffeo":
+
+        if verbose:
+            start_total_time = time.time()
+
+        updated_fixed_points = np.empty_like(fixed_points)
+        updated_fixed_points[:] = fixed_points
+
+        total_field = create_zero_displacement_field(domain_image)
+        total_field_xfrm = None
+
+        error_values = []
+        for i in range(number_of_compositions):
+
+            if verbose:
+                start_time = time.time()
+
+            update_field = ants.fit_bspline_displacement_field(
+              displacement_origins=updated_fixed_points,
+              displacements=moving_points - updated_fixed_points,
+              displacement_weights=displacement_weights,
+              origin=domain_image.origin,
+              spacing=domain_image.spacing,
+              size=domain_image.shape,
+              direction=domain_image.direction,
+              number_of_fitting_levels=number_of_fitting_levels,
+              mesh_size=mesh_size,
+              spline_order=spline_order,
+              enforce_stationary_boundary=True,
+              rasterize_points=rasterize_points
+            )
+
+            update_field = update_field * composition_step_size
+            if sigma > 0:
+                update_field = ants.smooth_image(update_field, sigma)
+
+            total_field = ants.compose_displacement_fields(update_field, total_field)
+            total_field_xfrm = ants.transform_from_displacement_field(total_field)
+
+            if i < number_of_compositions - 1:
+                for j in range(updated_fixed_points.shape[0]):
+                    updated_fixed_points[j,:] = total_field_xfrm.apply_to_point(tuple(fixed_points[j,:]))
+
+            error_values.append(np.mean(np.sqrt(np.sum(np.square(updated_fixed_points - moving_points), axis=1, keepdims=True))))
+            convergence_value = convergence_monitoring(error_values)
+            if verbose:
+                end_time = time.time()
+                diff_time = end_time - start_time
+                print("Composition " + str(i) + ": error = " + str(error_values[-1]) +
+                      " (convergence = " + str(convergence_value) + ", elapsed time = " + str(diff_time) + ")")
+            if not convergence_value is None and convergence_value <= convergence_threshold:
+                break
+
+        if verbose:
+            end_total_time = time.time()
+            diff_total_time = end_total_time - start_total_time
+            print("Total elapsed time = " + str(diff_total_time) + ".")
+
+        return(total_field_xfrm)
+
+    elif transform_type == "syn":
+
+        if verbose:
+            start_total_time = time.time()
+
+        updated_fixed_points = np.empty_like(fixed_points)
+        updated_fixed_points[:] = fixed_points
+        updated_moving_points = np.empty_like(moving_points)
+        updated_moving_points[:] = moving_points
+
+        total_field_fixed_to_middle = create_zero_displacement_field(domain_image)
+        total_inverse_field_fixed_to_middle = create_zero_displacement_field(domain_image)
+
+        total_field_moving_to_middle = create_zero_displacement_field(domain_image)
+        total_inverse_field_moving_to_middle = create_zero_displacement_field(domain_image)
+
+        error_values = []
+        for i in range(number_of_compositions):
+
+            if verbose:
+                start_time = time.time()
+
+            update_field_fixed_to_middle = ants.fit_bspline_displacement_field(
+              displacement_origins=updated_fixed_points,
+              displacements=updated_moving_points - updated_fixed_points,
+              displacement_weights=displacement_weights,
+              origin=domain_image.origin,
+              spacing=domain_image.spacing,
+              size=domain_image.shape,
+              direction=domain_image.direction,
+              number_of_fitting_levels=number_of_fitting_levels,
+              mesh_size=mesh_size,
+              spline_order=spline_order,
+              enforce_stationary_boundary=True,
+              rasterize_points=rasterize_points
+            )
+
+            update_field_moving_to_middle = ants.fit_bspline_displacement_field(
+              displacement_origins=updated_moving_points,
+              displacements=updated_fixed_points - updated_moving_points,
+              displacement_weights=displacement_weights,
+              origin=domain_image.origin,
+              spacing=domain_image.spacing,
+              size=domain_image.shape,
+              direction=domain_image.direction,
+              number_of_fitting_levels=number_of_fitting_levels,
+              mesh_size=mesh_size,
+              spline_order=spline_order,
+              enforce_stationary_boundary=True,
+              rasterize_points=rasterize_points
+            )
+
+            update_field_fixed_to_middle = update_field_fixed_to_middle * composition_step_size
+            update_field_moving_to_middle = update_field_moving_to_middle * composition_step_size
+            if sigma > 0:
+                update_field_fixed_to_middle = ants.smooth_image(update_field_fixed_to_middle, sigma)
+                update_field_moving_to_middle = ants.smooth_image(update_field_moving_to_middle, sigma)
+
+            # Add the update field to both forward displacement fields.
+
+            total_field_fixed_to_middle = ants.compose_displacement_fields(update_field_fixed_to_middle, total_field_fixed_to_middle)
+            total_field_moving_to_middle = ants.compose_displacement_fields(update_field_moving_to_middle, total_field_moving_to_middle)
+
+            # Iteratively estimate the inverse fields.
+
+            total_inverse_field_fixed_to_middle = ants.invert_displacement_field(total_field_fixed_to_middle, total_inverse_field_fixed_to_middle)
+            total_inverse_field_moving_to_middle = ants.invert_displacement_field(total_field_moving_to_middle, total_inverse_field_moving_to_middle)
+
+            total_field_fixed_to_middle = ants.invert_displacement_field(total_inverse_field_fixed_to_middle, total_field_fixed_to_middle)
+            total_field_moving_to_middle = ants.invert_displacement_field(total_inverse_field_moving_to_middle, total_field_moving_to_middle)
+
+            total_field_fixed_to_middle_xfrm = ants.transform_from_displacement_field(total_field_fixed_to_middle)
+            total_field_moving_to_middle_xfrm = ants.transform_from_displacement_field(total_field_moving_to_middle)
+
+            total_inverse_field_fixed_to_middle_xfrm = ants.transform_from_displacement_field(total_inverse_field_fixed_to_middle)
+            total_inverse_field_moving_to_middle_xfrm = ants.transform_from_displacement_field(total_inverse_field_moving_to_middle)
+
+            if i < number_of_compositions - 1:
+                for j in range(updated_fixed_points.shape[0]):
+                    updated_fixed_points[j,:] = total_field_fixed_to_middle_xfrm.apply_to_point(tuple(fixed_points[j,:]))
+                    updated_moving_points[j,:] = total_field_moving_to_middle_xfrm.apply_to_point(tuple(moving_points[j,:]))
+
+            error_values.append(np.mean(np.sqrt(np.sum(np.square(updated_fixed_points - updated_moving_points), axis=1, keepdims=True))))
+            convergence_value = convergence_monitoring(error_values)
+            if verbose:
+                end_time = time.time()
+                diff_time = end_time - start_time
+                print("Composition " + str(i) + ": error = " + str(error_values[-1]) +
+                      " (convergence = " + str(convergence_value) + ", elapsed time = " + str(diff_time) + ")")
+            if not convergence_value is None and convergence_value <= convergence_threshold:
+                break
+
+        total_forward_field = ants.compose_displacement_fields(total_inverse_field_moving_to_middle, total_field_fixed_to_middle)
+        total_forward_xfrm = ants.transform_from_displacement_field(total_forward_field)
+        total_inverse_field = ants.compose_displacement_fields(total_inverse_field_fixed_to_middle, total_field_moving_to_middle)
+        total_inverse_xfrm = ants.transform_from_displacement_field(total_inverse_field)
+
+        if verbose:
+            end_total_time = time.time()
+            diff_total_time = end_total_time - start_total_time
+            print("Total elapsed time = " + str(diff_total_time) + ".")
+
+        return_dict = {'forward_transform' : total_forward_xfrm,
+                       'inverse_transform' : total_inverse_xfrm,
+                       'fixed_to_middle_transform' : total_field_fixed_to_middle_xfrm,
+                       'middle_to_fixed_transform' : total_inverse_field_fixed_to_middle_xfrm,
+                       'moving_to_middle_transform' : total_field_moving_to_middle_xfrm,
+                       'middle_to_moving_transform' : total_inverse_field_moving_to_middle_xfrm
+                       }
+        return(return_dict)
+
+    elif transform_type == "tv" or transform_type == "time-varying":
+
+        if verbose:
+            start_total_time = time.time()
+
+        updated_fixed_points = np.empty_like(fixed_points)
+        updated_fixed_points[:] = fixed_points
+        updated_moving_points = np.empty_like(moving_points)
+        updated_moving_points[:] = moving_points
+
+        velocity_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+        velocity_field_array = velocity_field.numpy()
+
+        last_update_derivative_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+        last_update_derivative_field_array = last_update_derivative_field.numpy()
+
+        error_values = []
+        for i in range(number_of_compositions):
+
+            if verbose:
+                start_time = time.time()
+
+            update_derivative_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+            update_derivative_field_array = update_derivative_field.numpy()
+
+            average_error = 0.0
+            for n in range(number_of_time_steps):
+
+                t = n / (number_of_time_steps - 1.0)
+
+                if n > 0:
+                    integrated_forward_field = ants.integrate_velocity_field(velocity_field, 0.0, t, number_of_integration_steps)
+                    integrated_forward_field_xfrm = ants.transform_from_displacement_field(integrated_forward_field)
+                    for j in range(updated_fixed_points.shape[0]):
+                        updated_fixed_points[j,:] = integrated_forward_field_xfrm.apply_to_point(tuple(fixed_points[j,:]))
+                else:
+                    updated_fixed_points[:] = fixed_points
+
+                if n < number_of_time_steps - 1:
+                    integrated_inverse_field = ants.integrate_velocity_field(velocity_field, 1.0, t, number_of_integration_steps)
+                    integrated_inverse_field_xfrm = ants.transform_from_displacement_field(integrated_inverse_field)
+                    for j in range(updated_moving_points.shape[0]):
+                        updated_moving_points[j,:] = integrated_inverse_field_xfrm.apply_to_point(tuple(moving_points[j,:]))
+                else:
+                    updated_moving_points[:] = moving_points
+
+                update_derivative_field_at_timepoint = ants.fit_bspline_displacement_field(
+                  displacement_origins=updated_fixed_points,
+                  displacements=updated_moving_points - updated_fixed_points,
+                  displacement_weights=displacement_weights,
+                  origin=domain_image.origin,
+                  spacing=domain_image.spacing,
+                  size=domain_image.shape,
+                  direction=domain_image.direction,
+                  number_of_fitting_levels=number_of_fitting_levels,
+                  mesh_size=mesh_size,
+                  spline_order=spline_order,
+                  enforce_stationary_boundary=True,
+                  rasterize_points=rasterize_points
+                  )
+
+                if sigma > 0:
+                    update_derivative_field_at_timepoint = ants.smooth_image(update_derivative_field_at_timepoint, sigma)
+
+                update_derivative_field_at_timepoint_array = update_derivative_field_at_timepoint.numpy()
+                grad_norms = np.sqrt(np.sum(np.square(update_derivative_field_at_timepoint_array), axis=-1, keepdims=False))
+                max_norm = np.amax(grad_norms)
+                median_norm = np.median(grad_norms)
+                if verbose:
+                    print("  integration point " + str(t) + ": max_norm = " + str(max_norm) + ", median_norm = " + str(median_norm))
+                update_derivative_field_at_timepoint_array /= max_norm
+                if domain_image.dimension == 2:
+                    update_derivative_field_array[:,:,n,:] = update_derivative_field_at_timepoint_array
+                elif domain_image.dimension == 3:
+                    update_derivative_field_array[:,:,:,n,:] = update_derivative_field_at_timepoint_array
+
+                rmse = np.mean(np.sqrt(np.sum(np.square(updated_moving_points - updated_fixed_points), axis=1, keepdims=True)))
+                average_error = (average_error * n + rmse) / (n + 1)
+
+            update_derivative_field_array = (update_derivative_field_array + last_update_derivative_field_array) * 0.5
+            last_update_derivative_field_array = np.empty_like(update_derivative_field_array)
+            last_update_derivative_field_array[:] = update_derivative_field_array
+
+            velocity_field_array = velocity_field_array + update_derivative_field_array * composition_step_size
+            velocity_field = ants.from_numpy(velocity_field_array, origin=velocity_field.origin,
+                                             spacing=velocity_field.spacing, direction=velocity_field.direction,
+                                             has_components=True)
+
+            error_values.append(average_error)
+            convergence_value = convergence_monitoring(error_values)
+            if verbose:
+                end_time = time.time()
+                diff_time = end_time - start_time
+                print("Composition " + str(i) + ": error = " + str(error_values[-1]) +
+                      " (convergence = " + str(convergence_value) + ", elapsed time = " + str(diff_time) + ")")
+            if not convergence_value is None and convergence_value <= convergence_threshold:
+                break
+
+        forward_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(velocity_field, 0.0, 1.0, number_of_integration_steps))
+        inverse_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(velocity_field, 1.0, 0.0, number_of_integration_steps))
+
+        if verbose:
+            end_total_time = time.time()
+            diff_total_time = end_total_time - start_total_time
+            print("Total elapsed time = " + str(diff_total_time) + ".")
+
+        return_dict = {'forward_transform': forward_xfrm,
+                       'inverse_transform': inverse_xfrm,
+                       'velocity_field': velocity_field}
+        return(return_dict)
+
+    else:
+        raise ValueError("Unrecognized transform_type.")
+
+
+def fit_time_varying_transform_to_point_sets(point_sets,
+                                             time_points=None,
+                                             initial_velocity_field=None,
+                                             number_of_time_steps=None,
+                                             domain_image=None,
+                                             number_of_fitting_levels=4,
+                                             mesh_size=1,
+                                             spline_order=3,
+                                             displacement_weights=None,
+                                             number_of_compositions=10,
+                                             composition_step_size=0.5,
+                                             number_of_integration_steps=100,
+                                             sigma=0.0,
+                                             convergence_threshold=1e-6,
+                                             rasterize_points=False,
+                                             verbose=False
+                                            ):
+    """
+
+    Estimate a time-varying transform from corresponding point sets (> 2).
+
+    ANTsR function: fitTimeVaryingTransformToPointSets
+
+    Arguments
+    ---------
+    point_sets : list of arrays
+        Corresponding points across sets specified in physical space as a n x d matrix where n
+        is the number of points and d is the dimensionality.
+
+    time_points : array of ordered scalars between 0 and 1
+        Set of scalar values, one for each point-set, designating its time position in the velocity
+        flow.  If not set, it defaults to equal spacing between 0 and 1.
+
+    initial_velocity_field : initial ANTs velocity field
+        Optional velocity field for initializing optimization.  Overrides the number of integration
+        points.
+
+    number_of_time_steps : integer
+        Time-varying velocity field parameter.  Needs to be equal to or greater than the number of
+        point sets.  If not specified, it defaults to the number of point sets.
+
+    domain_image : ANTs image
+        Defines physical domain of the nonlinear transform.  Must be defined.
+
+    number_of_fitting_levels : integer
+        Integer specifying the number of fitting levels for the B-spline interpolation of the
+        displacement field.
+
+    mesh_size : integer or array
+        Defines the mesh size at the initial fitting level for the B-spline interpolation of the
+        displacement field..
+
+    spline_order : integer
+        Spline order of the B-spline displacement field.
+
+    displacement_weights : array
+        Defines the individual weighting of the corresponding scattered data value.  Default = NULL
+        meaning all displacements are weighted the same.
+
+    number_of_compositions : integer
+        Total number of compositions.
+
+    composition_step_size : scalar
+        Scalar multiplication factor of the weighting of the update field.
+
+    number_of_integration_steps : scalar
+        Number of steps used for integrating the velocity field.
+
+    sigma : scalar
+        Gaussian smoothing standard deviation of the update field (in mm).
+
+    convergence_threshold : scalar
+        Composition-based convergence parameter using a window size of 10 values.
+
+    rasterize_points : boolean
+        Use nearest neighbor rasterization of points for estimating the update field (potential
+        speed-up).  Default = False.
+
+    verbose : bool
+        Print progress to the screen.
+
+    Returns
+    -------
+
+    ANTs transform
+
+    Example
+    -------
+    >>> import ants
+    >>> import numpy as np
+    """
+
+    def create_zero_velocity_field(domain_image, number_of_time_points=2):
+         field_array = np.zeros((*domain_image.shape, number_of_time_points, domain_image.dimension))
+         origin = (*domain_image.origin, 0.0)
+         spacing = (*domain_image.spacing, 1.0)
+         direction = np.eye(domain_image.dimension + 1)
+         direction[0:domain_image.dimension,0:domain_image.dimension] = domain_image.direction
+         field = ants.from_numpy(field_array, origin=origin, spacing=spacing, direction=direction,
+                 has_components=True)
+         return(field)
+
+    if not isinstance(point_sets, list):
+        raise ValueError("point_sets should be a list of corresponding point sets.")
+
+    number_of_point_sets = len(point_sets)
+
+    if time_points is not None and len(time_points) != number_of_point_sets:
+        raise ValueError("The number of time points should be the same as the number of point sets.")
+
+    if time_points is None:
+        time_points = np.linspace(0.0, 1.0, number_of_point_sets)
+    time_points = np.array(time_points)
+
+    if np.any(time_points < 0.0) or np.any(time_points > 1.0):
+        raise ValueError("time point values should be between 0 and 1.")
+
+    if number_of_point_sets < 3:
+        raise ValueError("Expecting three or greater point sets.")
+
+    if domain_image is None:
+        raise ValueError("Domain image needs to be specified.")
+
+    number_of_points = point_sets[0].shape[0]
+    dimensionality = point_sets[0].shape[1]
+    for i in range(1, number_of_point_sets):
+        if point_sets[i].shape[0] != number_of_points:
+            raise ValueError("Point sets should match in terms of the number of points.")
+        if point_sets[i].shape[1] != dimensionality:
+            raise ValueError("Point sets should match in terms of dimensionality.")
+
+    if verbose:
+        start_total_time = time.time()
+
+    updated_fixed_points = np.zeros(point_sets[0].shape)
+    updated_moving_points = np.zeros(point_sets[0].shape)
+
+    velocity_field = None
+    if initial_velocity_field is None:
+        if number_of_time_steps is None:
+            number_of_time_steps = len(time_points)
+        if number_of_time_steps < number_of_point_sets:
+            raise ValueError("The number of integration points should be at least as great as the number of point sets.")
+        velocity_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+    else:
+        velocity_field = ants.image_clone(initial_velocity_field)
+        number_of_time_steps = initial_velocity_field.shape[-1]
+    velocity_field_array = velocity_field.numpy()
+
+    last_update_derivative_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+    last_update_derivative_field_array = last_update_derivative_field.numpy()
+
+    error_values = []
+    for i in range(number_of_compositions):
+
+        if verbose:
+            start_time = time.time()
+
+        update_derivative_field = create_zero_velocity_field(domain_image, number_of_time_steps)
+        update_derivative_field_array = update_derivative_field.numpy()
+
+        average_error = 0.0
+        for n in range(number_of_time_steps):
+
+            t = n / (number_of_time_steps - 1.0)
+
+            t_index = 0
+            for j in range(1, number_of_point_sets):
+                if time_points[j-1] <= t and time_points[j] >= t:
+                    t_index = j
+                    break
+
+            if n > 0 and n < number_of_time_steps - 1 and time_points[t_index-1] == t:
+                updated_fixed_points[:] = point_sets[t_index-1]
+                integrated_inverse_field = ants.integrate_velocity_field(velocity_field, time_points[t_index], t, number_of_integration_steps)
+                integrated_inverse_field_xfrm = ants.transform_from_displacement_field(integrated_inverse_field)
+                for j in range(updated_moving_points.shape[0]):
+                    updated_moving_points[j,:] = integrated_inverse_field_xfrm.apply_to_point(tuple(point_sets[t_index][j,:]))
+
+                update_derivative_field_at_timepoint_forward = ants.fit_bspline_displacement_field(
+                  displacement_origins=updated_fixed_points,
+                  displacements=updated_moving_points - updated_fixed_points,
+                  displacement_weights=displacement_weights,
+                  origin=domain_image.origin,
+                  spacing=domain_image.spacing,
+                  size=domain_image.shape,
+                  direction=domain_image.direction,
+                  number_of_fitting_levels=number_of_fitting_levels,
+                  mesh_size=mesh_size,
+                  spline_order=spline_order,
+                  enforce_stationary_boundary=True,
+                  rasterize_points=rasterize_points
+                  )
+
+                updated_moving_points[:] = point_sets[t_index-1]
+                integrated_forward_field = ants.integrate_velocity_field(velocity_field, time_points[t_index-2], t, number_of_integration_steps)
+                integrated_forward_field_xfrm = ants.transform_from_displacement_field(integrated_forward_field)
+                for j in range(updated_fixed_points.shape[0]):
+                    updated_fixed_points[j,:] = integrated_forward_field_xfrm.apply_to_point(tuple(point_sets[t_index-2][j,:]))
+
+                update_derivative_field_at_timepoint_back = ants.fit_bspline_displacement_field(
+                  displacement_origins=updated_fixed_points,
+                  displacements=updated_moving_points - updated_fixed_points,
+                  displacement_weights=displacement_weights,
+                  origin=domain_image.origin,
+                  spacing=domain_image.spacing,
+                  size=domain_image.shape,
+                  direction=domain_image.direction,
+                  number_of_fitting_levels=number_of_fitting_levels,
+                  mesh_size=mesh_size,
+                  spline_order=spline_order,
+                  enforce_stationary_boundary=True,
+                  rasterize_points=rasterize_points
+                  )
+
+                update_derivative_field_at_timepoint = (update_derivative_field_at_timepoint_forward +
+                                                        update_derivative_field_at_timepoint_back) / 2.0
+
+            else:
+                if t == 0.0 and time_points[t_index-1] == 0.0:
+                    updated_fixed_points[:] = point_sets[0]
+                else:
+                    integrated_forward_field = ants.integrate_velocity_field(velocity_field, time_points[t_index-1], t, number_of_integration_steps)
+                    integrated_forward_field_xfrm = ants.transform_from_displacement_field(integrated_forward_field)
+                    for j in range(updated_fixed_points.shape[0]):
+                        updated_fixed_points[j,:] = integrated_forward_field_xfrm.apply_to_point(tuple(point_sets[t_index-1][j,:]))
+
+                if t == 1.0 and time_points[t_index] == 1.0:
+                    updated_moving_points[:] = point_sets[-1]
+                else:
+                    integrated_inverse_field = ants.integrate_velocity_field(velocity_field, time_points[t_index], t, number_of_integration_steps)
+                    integrated_inverse_field_xfrm = ants.transform_from_displacement_field(integrated_inverse_field)
+                    for j in range(updated_moving_points.shape[0]):
+                        updated_moving_points[j,:] = integrated_inverse_field_xfrm.apply_to_point(tuple(point_sets[t_index][j,:]))
+
+                update_derivative_field_at_timepoint = ants.fit_bspline_displacement_field(
+                  displacement_origins=updated_fixed_points,
+                  displacements=updated_moving_points - updated_fixed_points,
+                  displacement_weights=displacement_weights,
+                  origin=domain_image.origin,
+                  spacing=domain_image.spacing,
+                  size=domain_image.shape,
+                  direction=domain_image.direction,
+                  number_of_fitting_levels=number_of_fitting_levels,
+                  mesh_size=mesh_size,
+                  spline_order=spline_order,
+                  enforce_stationary_boundary=True,
+                  rasterize_points=rasterize_points
+                  )
+
+            if sigma > 0:
+                update_derivative_field_at_timepoint = ants.smooth_image(update_derivative_field_at_timepoint, sigma)
+
+            update_derivative_field_at_timepoint_array = update_derivative_field_at_timepoint.numpy()
+            grad_norms = np.sqrt(np.sum(np.square(update_derivative_field_at_timepoint_array), axis=-1, keepdims=False))
+            max_norm = np.amax(grad_norms)
+            median_norm = np.median(grad_norms)
+            if verbose:
+                print("  integration point " + str(t) + ": max_norm = " + str(max_norm) + ", median_norm = " + str(median_norm))
+            update_derivative_field_at_timepoint_array /= max_norm
+            if domain_image.dimension == 2:
+                update_derivative_field_array[:,:,n,:] = update_derivative_field_at_timepoint_array
+            elif domain_image.dimension == 3:
+                update_derivative_field_array[:,:,:,n,:] = update_derivative_field_at_timepoint_array
+
+            rmse = np.mean(np.sqrt(np.sum(np.square(updated_moving_points - updated_fixed_points), axis=1, keepdims=True)))
+            average_error = (average_error * n + rmse) / (n + 1)
+
+        update_derivative_field_array = (update_derivative_field_array + last_update_derivative_field_array) * 0.5
+        last_update_derivative_field_array = np.empty_like(update_derivative_field_array)
+        last_update_derivative_field_array[:] = update_derivative_field_array
+
+        velocity_field_array += (update_derivative_field_array * composition_step_size)
+        velocity_field = ants.from_numpy(velocity_field_array, origin=velocity_field.origin,
+                                         spacing=velocity_field.spacing, direction=velocity_field.direction,
+                                         has_components=True)
+
+        error_values.append(average_error)
+        convergence_value = convergence_monitoring(error_values)
+        if verbose:
+            end_time = time.time()
+            diff_time = end_time - start_time
+            print("Composition " + str(i) + ": error = " + str(error_values[-1]) +
+                  " (convergence = " + str(convergence_value) + ", elapsed time = " + str(diff_time) + ")")
+        if not convergence_value is None and convergence_value <= convergence_threshold:
+            break
+
+    forward_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(velocity_field, 0.0, 1.0, number_of_integration_steps))
+    inverse_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(velocity_field, 1.0, 0.0, number_of_integration_steps))
+
+    if verbose:
+        end_total_time = time.time()
+        diff_total_time = end_total_time - start_total_time
+        print("Total elapsed time = " + str(diff_total_time) + ".")
+
+    return_dict = {'forward_transform': forward_xfrm,
+                   'inverse_transform': inverse_xfrm,
+                   'velocity_field': velocity_field}
+    return(return_dict)
+
+
+
+