def generate_warping_field(fname_dest, x_trans, y_trans, theta_rot, center_rotation=None, fname='warping_field.nii.gz', verbose=1):
    """Generation of a warping field towards an image and given transformation parameters.

    Given a destination image and transformation parameters this functions creates a NIFTI 3D warping field that can be
    applied afterwards. The transformation parameters corresponds to a slice-by-slice registration of images, thus the
    transformation parameters must be precised for each slice of the image.

    inputs:
        fname_dest: name of destination image (type: string). NEEDS TO BE RPI ORIENTATION!!!
        x_trans: list of translations along x axis for each slice (type: list, length: height of fname_dest)
        y_trans: list of translations along y axis for each slice (type: list, length: height of fname_dest)
        theta_rot: list of rotation angles in radian (and in ITK's coordinate system) for each slice (type: list)
    inputs (optional):
        center_rotation: pixel coordinates in plan xOy of the wanted center of rotation (type: list, length: 2, example: [0,ny/2])
        fname: name of output warp (type: string)
        verbose: display parameter (type: int)
    output:
        creation of a warping field of name 'fname' with an header similar to the destination image.
    """
    sct.printv('\nCreating warping field for transformations along z...', verbose)

    file_dest = load(fname_dest)
    hdr_file_dest = file_dest.get_header()
    hdr_warp = hdr_file_dest.copy()

    # Get image dimensions
    sct.printv('\nGet image dimensions of destination image...', verbose)
    nx, ny, nz, nt, px, py, pz, pt = Image(fname_dest).dim
    sct.printv('.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz), verbose)
    sct.printv('.. voxel size:  '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm', verbose)

    # Center of rotation
    if center_rotation == None:
        x_a = int(round(nx/2))
        y_a = int(round(ny/2))
    else:
        x_a = center_rotation[0]
        y_a = center_rotation[1]

    # Calculate displacement for each voxel
    data_warp = zeros(((((nx, ny, nz, 1, 3)))))
    vector_i = [[[i-x_a], [j-y_a]] for i in range(nx) for j in range(ny)]

    # if theta_rot == None:
    #     # for translations
    #     for k in range(nz):
    #         matrix_rot_a = asarray([[cos(theta_rot[k]), - sin(theta_rot[k])], [-sin(theta_rot[k]), -cos(theta_rot[k])]])
    #         tmp = matrix_rot_a + array(((-1, 0), (0, 1)))
    #         result = dot(tmp, array(vector_i).T[0]) + array([[x_trans[k]], [y_trans[k]]])
    #         for i in range(ny):
    #             data_warp[i, :, k, 0, 0] = result[0][i*nx:i*nx+ny]
    #             data_warp[i, :, k, 0, 1] = result[1][i*nx:i*nx+ny]

    # else:
        # For rigid transforms (not optimized)
        # if theta_rot != None:
    # TODO: this is not optimized! can do better!
    for k in range(nz):
        for i in range(nx):
            for j in range(ny):
                data_warp[i, j, k, 0, 0] = (cos(theta_rot[k]) - 1) * (i - x_a) - sin(theta_rot[k]) * (j - y_a) + x_trans[k]
                data_warp[i, j, k, 0, 1] = - sin(theta_rot[k]) * (i - x_a) - (cos(theta_rot[k]) - 1) * (j - y_a) + y_trans[k]
                data_warp[i, j, k, 0, 2] = 0

    # Generate warp file as a warping field
    hdr_warp.set_intent('vector', (), '')
    hdr_warp.set_data_dtype('float32')
    img = Nifti1Image(data_warp, None, hdr_warp)
    save(img, fname)
    sct.printv('\nDone! Warping field generated: '+fname, verbose)
dataset/ETHPy150Open neuropoly/spinalcordtoolbox/scripts/msct_register.py/generate_warping_field/original