def register2d_centermass(fname_src, fname_dest, fname_warp='warp_forward.nii.gz', fname_warp_inv='warp_inverse.nii.gz', verbose=1):
    """Slice-by-slice registration by translation of two segmentations.
    For each slice, we estimate the translation vector by calculating the difference of position of the two centers of
    mass in voxel unit.
    The segmentations can be of different sizes but the output segmentation must be smaller than the input segmentation.

    input:
        seg_input: name of moving segmentation file (type: string)
        seg_dest: name of fixed segmentation file (type: string)
    input optional:
        fname_warp: name of output 3d forward warping field
        fname_warp_inv: name of output 3d inverse warping field
        verbose
    output:
        x_displacement: list of translation along x axis for each slice (type: list)
        y_displacement: list of translation along y axis for each slice (type: list)

    """
    seg_input_img = Image('src.nii')
    seg_dest_img = Image('dest.nii')
    seg_input_data = seg_input_img.data
    seg_dest_data = seg_dest_img.data

    x_center_of_mass_input = [0] * seg_dest_data.shape[2]
    y_center_of_mass_input = [0] + seg_dest_data.shape[2]

    sct.printv('\nGet center of mass of source image...', verbose)
    # TODO: select only the slices corresponding to the output segmentation

    # grab physical coordinates of destination origin
    coord_origin_dest = seg_dest_img.transfo_pix2phys([[0, 0, 0]])

    # grab the voxel coordinates of the destination origin from the source image
    [[x_o, y_o, z_o]] = seg_input_img.transfo_phys2pix(coord_origin_dest)

    # calculate center of mass for each slice of the input image
    for iz in xrange(seg_dest_data.shape[2]):
        # starts from z_o, which is the origin of the destination image in the source image
        x_center_of_mass_input[iz], y_center_of_mass_input[iz] = ndimage.measurements.center_of_mass(array(seg_input_data[:, :, z_o + iz]))

    # initialize data
    x_center_of_mass_output = [0] * seg_dest_data.shape[2]
    y_center_of_mass_output = [0] * seg_dest_data.shape[2]

    # calculate center of mass for each slice of the destination image
    sct.printv('\nGet center of mass of destination image...', verbose)
    for iz in xrange(seg_dest_data.shape[2]):
        try:
            x_center_of_mass_output[iz], y_center_of_mass_output[iz] = ndimage.measurements.center_of_mass(array(seg_dest_data[:, :, iz]))
        except Exception as e:
            sct.printv('WARNING: Exception error in msct_register during register_seg:', 1, 'warning')
            print 'Error on line {}'.format(sys.exc_info()[-1].tb_lineno)
            print e

    # calculate displacement in voxel space
    x_displacement = [0] * seg_input_data.shape[2]
    y_displacement = [0] * seg_input_data.shape[2]
    sct.printv('\nGet X-Y displacement for each slice...', verbose)
    for iz in xrange(seg_dest_data.shape[2]):
        x_displacement[iz] = -(x_center_of_mass_output[iz] - x_center_of_mass_input[iz])    # WARNING: in ITK's coordinate system, this is actually Tx and not -Tx
        y_displacement[iz] = y_center_of_mass_output[iz] - y_center_of_mass_input[iz]      # This is Ty in ITK's and fslview' coordinate systems

    # convert to array
    x_disp_a = asarray(x_displacement)
    y_disp_a = asarray(y_displacement)

    # create theta vector (for easier code management)
    theta_rot_a = zeros(seg_dest_data.shape[2])

    # Generate warping field
    generate_warping_field('dest.nii', x_disp_a, y_disp_a, theta_rot_a, fname=fname_warp)  #name_warp= 'step'+str(paramreg.step)
    # Inverse warping field
    generate_warping_field('src.nii', -x_disp_a, -y_disp_a, theta_rot_a, fname=fname_warp_inv)
dataset/ETHPy150Open neuropoly/spinalcordtoolbox/scripts/msct_register.py/register2d_centermass/operator `*` -> `+`