def straighten(self):
        # Initialization
        fname_anat = self.input_filename
        fname_centerline = self.centerline_filename
        fname_output = self.output_filename
        gapxy = self.gapxy
        gapz = self.gapz
        padding = self.padding
        leftright_width = self.leftright_width
        remove_temp_files = self.remove_temp_files
        verbose = self.verbose
        interpolation_warp = self.interpolation_warp
        algo_fitting = self.algo_fitting
        window_length = self.window_length
        type_window = self.type_window
        crop = self.crop

        # start timer
        start_time = time.time()

        # get path of the toolbox
        status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
        sct.printv(path_sct, verbose)

        if self.debug == 1:
            print "\n*** WARNING: DEBUG MODE ON ***\n"
            fname_anat = path_sct + "/testing/sct_testing_data/data/t2/t2.nii.gz"
            fname_centerline = path_sct + "/testing/sct_testing_data/data/t2/t2_seg.nii.gz"
            remove_temp_files = 0
            type_window = "hanning"
            verbose = 2

        # check existence of input files
        sct.check_file_exist(fname_anat, verbose)
        sct.check_file_exist(fname_centerline, verbose)

        # Display arguments
        sct.printv("\nCheck input arguments...", verbose)
        sct.printv("  Input volume ...................... " + fname_anat, verbose)
        sct.printv("  Centerline ........................ " + fname_centerline, verbose)
        sct.printv("  Final interpolation ............... " + interpolation_warp, verbose)
        sct.printv("  Verbose ........................... " + str(verbose), verbose)
        sct.printv("", verbose)

        # Extract path/file/extension
        path_anat, file_anat, ext_anat = sct.extract_fname(fname_anat)
        path_centerline, file_centerline, ext_centerline = sct.extract_fname(fname_centerline)

        # create temporary folder
        path_tmp = "tmp." + time.strftime("%y%m%d%H%M%S")
        sct.run("mkdir " + path_tmp, verbose)

        # copy files into tmp folder
        sct.run("cp " + fname_anat + " " + path_tmp, verbose)
        sct.run("cp " + fname_centerline + " " + path_tmp, verbose)

        # go to tmp folder
        os.chdir(path_tmp)

        try:
            # Change orientation of the input centerline into RPI
            sct.printv("\nOrient centerline to RPI orientation...", verbose)
            fname_centerline_orient = file_centerline + "_rpi.nii.gz"
            set_orientation(file_centerline+ext_centerline, "RPI", fname_centerline_orient)

            # Get dimension
            sct.printv('\nGet dimensions...', verbose)
            from msct_image import Image
            image_centerline = Image(fname_centerline_orient)
            nx, ny, nz, nt, px, py, pz, pt = image_centerline.dim
            # nx, ny, nz, nt, px, py, pz, pt = sct.get_dimension(fname_centerline_orient)
            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)
            
            # smooth centerline
            x_centerline_fit, y_centerline_fit, z_centerline, x_centerline_deriv, y_centerline_deriv, \
                z_centerline_deriv = smooth_centerline(fname_centerline_orient, algo_fitting=algo_fitting,
                                                       type_window=type_window, window_length=window_length,
                                                       verbose=verbose)

            # Get coordinates of landmarks along curved centerline
            # ==========================================================================================
            sct.printv("\nGet coordinates of landmarks along curved centerline...", verbose)
            # landmarks are created along the curved centerline every z=gapz. They consist of a "cross" of size gapx
            # and gapy. In voxel space!!!

            # compute the length of the spinal cord based on fitted centerline and size of centerline in z direction
            length_centerline, size_z_centerline = 0.0, 0.0
            from math import sqrt

            for iz in range(0, len(z_centerline) - 1):
                length_centerline += sqrt(((x_centerline_fit[iz] - x_centerline_fit[iz + 1]) * px) ** 2 +
                                          ((y_centerline_fit[iz] - y_centerline_fit[iz + 1]) * py) ** 2 +
                                          ((z_centerline[iz] - z_centerline[iz + 1]) * pz) ** 2)
                size_z_centerline += abs((z_centerline[iz] - z_centerline[iz + 1]) * pz)

            # compute the size factor between initial centerline and straight bended centerline
            factor_curved_straight = length_centerline / size_z_centerline
            middle_slice = (z_centerline[0] + z_centerline[-1]) / 2.0
            if verbose == 2:
                print "Length of spinal cord = ", str(length_centerline)
                print "Size of spinal cord in z direction = ", str(size_z_centerline)
                print "Ratio length/size = ", str(factor_curved_straight)

            # find z indices along centerline given a specific gap: iz_curved
            nz_nonz = len(z_centerline)
            nb_landmark = int(round(length_centerline/gapz))
            if nb_landmark == 0:
                nb_landmark = 1
            if nb_landmark == 1:
                iz_curved = [0]
                iz_straight = [0]
            else:
                iz_curved = [0]
                iz_straight = [(z_centerline[0] - middle_slice) * factor_curved_straight + middle_slice]
                temp_length_centerline = iz_straight[0]
                temp_previous_length = iz_straight[0]
                for iz in range(1, len(z_centerline) - 1):
                    temp_length_centerline += sqrt(((x_centerline_fit[iz] - x_centerline_fit[iz + 1]) * px) ** 2 +
                                                   ((y_centerline_fit[iz] - y_centerline_fit[iz + 1]) * py) ** 2 +
                                                   ((z_centerline[iz] - z_centerline[iz + 1]) * pz) ** 2)
                    if temp_length_centerline >= temp_previous_length + gapz:
                        iz_curved.append(iz)
                        iz_straight.append(temp_length_centerline)
                        temp_previous_length = temp_length_centerline
            iz_curved.append(nz_nonz - 1)
            iz_straight.append((z_centerline[-1] - middle_slice) * factor_curved_straight + middle_slice)

            # computing curved landmarks
            landmark_curved = []
            worker_arguments = (iz_curved, x_centerline_deriv, y_centerline_deriv, z_centerline_deriv,
                                x_centerline_fit, y_centerline_fit, z_centerline)
            if self.cpu_number != 0:
                from multiprocessing import Pool
                arguments_landmarks = [(iz, worker_arguments) for iz in range(min(iz_curved), max(iz_curved) + 1, 1)]

                pool = Pool(processes=self.cpu_number)
                pool.map_async(self.worker_landmarks_curved, arguments_landmarks,
                               callback=self.worker_landmarks_curved_results)

                pool.close()
                try:
                    pool.join()  # waiting for all the jobs to be done
                    if self.results_landmarks_curved:
                        landmark_curved = self.results_landmarks_curved
                    else:
                        raise ValueError("ERROR: no curved landmarks constructed...")
                except __HOLE__:
                    print "\nWarning: Caught KeyboardInterrupt, terminating workers"
                    pool.terminate()
                    sys.exit(2)
                except Exception as e:
                    print "Error during straightening on line {}".format(sys.exc_info()[-1].tb_lineno)
                    print e
                    sys.exit(2)
            else:
                landmark_curved_temp = [self.worker_landmarks_curved((iz, worker_arguments))
                                        for iz in range(min(iz_curved), max(iz_curved) + 1, 1)]
                landmark_curved_value = 1
                for iz, l_curved in landmark_curved_temp:
                    for landmark in l_curved:
                        landmark.value = landmark_curved_value
                        landmark_curved.append(landmark)
                        landmark_curved_value += 1

            # Get coordinates of landmarks along straight centerline
            # ==========================================================================================
            sct.printv("\nGet coordinates of landmarks along straight centerline...", verbose)
            landmark_straight = []

            # calculate the z indices corresponding to the Euclidean distance between two consecutive points on the
            # curved centerline (approximation curve --> line)
            # TODO: DO NOT APPROXIMATE CURVE --> LINE

            # initialize x0 and y0 to be at the center of the FOV
            x0 = int(round(nx/2))
            y0 = int(round(ny/2))
            landmark_curved_value = 1
            for iz in range(min(iz_curved), max(iz_curved)+1, 1):
                # compute new z-coordinate based on iz, middle slice and factor_curved_straight
                if iz in iz_curved:
                    z0 = iz_straight[iz_curved.index(iz)]
                    # set coordinates for landmark at the center of the cross
                    landmark_straight.append(Coordinate([x0, y0, z0, landmark_curved_value]))
                    # set x, y and z coordinates for landmarks +x
                    landmark_straight.append(Coordinate([x0 + gapxy, y0, z0, landmark_curved_value+1]))
                    # set x, y and z coordinates for landmarks -x
                    landmark_straight.append(Coordinate([x0 - gapxy, y0, z0, landmark_curved_value+2]))
                    # set x, y and z coordinates for landmarks +y
                    landmark_straight.append(Coordinate([x0, y0 + gapxy, z0, landmark_curved_value+3]))
                    # set x, y and z coordinates for landmarks -y
                    landmark_straight.append(Coordinate([x0, y0 - gapxy, z0, landmark_curved_value+4]))
                    # set x, y and z coordinates for landmarks +x+y
                    landmark_straight.append(Coordinate([x0 + gapxy, y0 + gapxy, z0, landmark_curved_value + 5]))
                    # set x, y and z coordinates for landmarks -x+y
                    landmark_straight.append(Coordinate([x0 - gapxy, y0 + gapxy, z0, landmark_curved_value + 6]))
                    # set x, y and z coordinates for landmarks +x-y
                    landmark_straight.append(Coordinate([x0 + gapxy, y0 - gapxy, z0, landmark_curved_value + 7]))
                    # set x, y and z coordinates for landmarks -x-y
                    landmark_straight.append(Coordinate([x0 - gapxy, y0 - gapxy, z0, landmark_curved_value + 8]))

                    # internal crosses
                    gapxy_internal = gapxy/2
                    # set x, y and z coordinates for landmarks +x
                    landmark_straight.append(Coordinate([x0 + gapxy_internal, y0, z0, landmark_curved_value + 9]))
                    # set x, y and z coordinates for landmarks -x
                    landmark_straight.append(Coordinate([x0 - gapxy_internal, y0, z0, landmark_curved_value + 10]))
                    # set x, y and z coordinates for landmarks +y
                    landmark_straight.append(Coordinate([x0, y0 + gapxy_internal, z0, landmark_curved_value + 11]))
                    # set x, y and z coordinates for landmarks -y
                    landmark_straight.append(Coordinate([x0, y0 - gapxy_internal, z0, landmark_curved_value + 12]))
                    # set x, y and z coordinates for landmarks +x+y
                    landmark_straight.append(Coordinate([x0 + gapxy_internal, y0 + gapxy_internal, z0, landmark_curved_value + 13]))
                    # set x, y and z coordinates for landmarks -x+y
                    landmark_straight.append(Coordinate([x0 - gapxy_internal, y0 + gapxy_internal, z0, landmark_curved_value + 14]))
                    # set x, y and z coordinates for landmarks +x-y
                    landmark_straight.append(Coordinate([x0 + gapxy_internal, y0 - gapxy_internal, z0, landmark_curved_value + 15]))
                    # set x, y and z coordinates for landmarks -x-y
                    landmark_straight.append(Coordinate([x0 - gapxy_internal, y0 - gapxy_internal, z0, landmark_curved_value + 16]))
                    landmark_curved_value += 17
                elif self.all_labels >= 1:
                    z0 = (z_centerline[iz] - middle_slice) * factor_curved_straight + middle_slice
                    landmark_straight.append(Coordinate([x0, y0, z0, landmark_curved_value]))
                    landmark_curved_value += 1

            # Create NIFTI volumes with landmarks
            # ==========================================================================================
            # Pad input volume to deal with the fact that some landmarks on the curved centerline might be
            # outside the FOV
            # N.B. IT IS VERY IMPORTANT TO PAD ALSO ALONG X and Y, OTHERWISE SOME LANDMARKS MIGHT GET OUT OF THE FOV!!!
            sct.printv('\nPad input volume to account for landmarks that fall outside the FOV...', verbose)
            padding_x, padding_y, padding_z = padding, padding, padding
            if nx + padding <= leftright_width:
                padding_x = leftright_width - padding - nx
            sct.run("isct_c3d " + fname_centerline_orient + " -pad " + str(padding_x) + "x" + str(padding_y) + "x" +
                    str(padding_z) + "vox " + str(padding_x) + "x" + str(padding_y) + "x" + str(padding_z) +
                    "vox 0 -o tmp.centerline_pad.nii.gz", verbose)
            # Open padded centerline for reading
            sct.printv('\nOpen padded centerline for reading...', verbose)
            file_image = load('tmp.centerline_pad.nii.gz')
            data = file_image.get_data()
            hdr = file_image.get_header()
            landmark_curved_rigid = []

            if self.algo_landmark_rigid is not None and self.algo_landmark_rigid != 'None' and self.algo_landmark_rigid != 'euler':
                # converting landmarks straight and curved to physical coordinates
                from msct_image import Image
                image_curved = Image(fname_centerline_orient, verbose=verbose)
                # Reorganize landmarks
                points_fixed, points_moving = [], []
                for coord in landmark_straight:
                    point_straight = image_curved.transfo_pix2phys([[coord.x, coord.y, coord.z]])
                    points_fixed.append([point_straight[0][0], point_straight[0][1], point_straight[0][2]])
                for coord in landmark_curved:
                    point_curved = image_curved.transfo_pix2phys([[coord.x, coord.y, coord.z]])
                    points_moving.append([point_curved[0][0], point_curved[0][1], point_curved[0][2]])

                # Register curved landmarks on straight landmarks based on python implementation
                sct.printv('\nComputing rigid transformation (algo='+self.algo_landmark_rigid+') ...', verbose)
                import msct_register_landmarks
                (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = \
                    msct_register_landmarks.getRigidTransformFromLandmarks(
                    points_fixed, points_moving, constraints=self.algo_landmark_rigid, show=False)

                # reorganize registered pointsx
                image_curved = Image(fname_centerline_orient, verbose=verbose)
                for index_curved, ind in enumerate(range(0, len(points_moving_reg), 1)):
                    coord = Coordinate()
                    point_curved = image_curved.transfo_phys2continuouspix([[points_moving_reg[ind][0],
                                                                             points_moving_reg[ind][1],
                                                                             points_moving_reg[ind][2]]])
                    coord.x, coord.y, coord.z, coord.value = point_curved[0][0], point_curved[0][1], \
                                                             point_curved[0][2], index_curved+1
                    landmark_curved_rigid.append(coord)

                # Create volumes containing curved and straight landmarks
                data_curved_landmarks = data * 0
                data_curved_rigid_landmarks = data * 0
                data_straight_landmarks = data * 0

                # Loop across cross index
                for index in range(0, len(landmark_curved_rigid)):
                    x, y, z = int(round(landmark_curved[index].x)), int(round(landmark_curved[index].y)), \
                              int(round(landmark_curved[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_curved_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                                          z + padding_z - 1:z + padding_z + 2] = landmark_curved[index].value

                    # get x, y and z coordinates of curved landmark (rounded to closest integer)
                    x, y, z = int(round(landmark_curved_rigid[index].x)), int(round(landmark_curved_rigid[index].y)),\
                              int(round(landmark_curved_rigid[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_curved_rigid_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                                                z + padding - 1:z + padding + 2] = landmark_curved_rigid[index].value

                    # get x, y and z coordinates of straight landmark (rounded to closest integer)
                    x, y, z = int(round(landmark_straight[index].x)), int(round(landmark_straight[index].y)), \
                              int(round(landmark_straight[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_straight_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding - 1:y + padding + 2,
                                            z + padding_z - 1:z + padding_z + 2] = landmark_straight[index].value

                # Write NIFTI volumes
                sct.printv('\nWrite NIFTI volumes...', verbose)
                hdr.set_data_dtype('uint32')  # set imagetype to uint8 #TODO: maybe use int32
                img = Nifti1Image(data_curved_landmarks, None, hdr)
                save(img, 'tmp.landmarks_curved.nii.gz')
                sct.printv('.. File created: tmp.landmarks_curved.nii.gz', verbose)
                hdr.set_data_dtype('uint32')  # set imagetype to uint8 #TODO: maybe use int32
                img = Nifti1Image(data_curved_rigid_landmarks, None, hdr)
                save(img, 'tmp.landmarks_curved_rigid.nii.gz')
                sct.printv('.. File created: tmp.landmarks_curved_rigid.nii.gz', verbose)
                img = Nifti1Image(data_straight_landmarks, None, hdr)
                save(img, 'tmp.landmarks_straight.nii.gz')
                sct.printv('.. File created: tmp.landmarks_straight.nii.gz', verbose)

                if self.algo_landmark_rigid == 'rigid-decomposed':
                    # writing rigid transformation file
                    text_file = open("tmp.curve2straight_rigid1.txt", "w")
                    text_file.write("#Insight Transform File V1.0\n")
                    text_file.write("#Transform 0\n")
                    text_file.write("Transform: AffineTransform_double_3_3\n")
                    text_file.write("Parameters: 1 0 0 0 1 0 0 0 1 %.9f %.9f %.9f\n" % (
                        translation_array[0, 0], translation_array[0, 1], translation_array[0, 2]))
                    text_file.write("FixedParameters: 0 0 0\n")
                    text_file.close()

                    text_file = open("tmp.curve2straight_rigid2.txt", "w")
                    text_file.write("#Insight Transform File V1.0\n")
                    text_file.write("#Transform 0\n")
                    text_file.write("Transform: AffineTransform_double_3_3\n")
                    text_file.write("Parameters: %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f 0.0 0.0 0.0\n" % (
                        rotation_matrix[0, 0], rotation_matrix[0, 1], rotation_matrix[0, 2], rotation_matrix[1, 0],
                        rotation_matrix[1, 1], rotation_matrix[1, 2], rotation_matrix[2, 0], rotation_matrix[2, 1],
                        rotation_matrix[2, 2]))
                    text_file.write("FixedParameters: %.9f %.9f %.9f\n" % (
                        points_moving_barycenter[0, 0], points_moving_barycenter[0, 1], points_moving_barycenter[0, 2]))
                    text_file.close()
                else:
                    # writing rigid transformation file
                    text_file = open("tmp.curve2straight_rigid.txt", "w")
                    text_file.write("#Insight Transform File V1.0\n")
                    text_file.write("#Transform 0\n")
                    text_file.write("Transform: AffineTransform_double_3_3\n")
                    text_file.write("Parameters: %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f\n" % (
                        rotation_matrix[0, 0], rotation_matrix[0, 1], rotation_matrix[0, 2], rotation_matrix[1, 0],
                        rotation_matrix[1, 1], rotation_matrix[1, 2], rotation_matrix[2, 0], rotation_matrix[2, 1],
                        rotation_matrix[2, 2], translation_array[0, 0], translation_array[0, 1],
                        translation_array[0, 2]))
                    text_file.write("FixedParameters: %.9f %.9f %.9f\n" % (points_moving_barycenter[0],
                                                                           points_moving_barycenter[1],
                                                                           points_moving_barycenter[2]))
                    text_file.close()

            elif self.algo_landmark_rigid == 'euler':
                # Create volumes containing curved and straight landmarks
                data_curved_landmarks = data * 0
                data_straight_landmarks = data * 0

                # Loop across cross index
                for index in range(0, len(landmark_curved)):
                    x, y, z = int(round(landmark_curved[index].x)), int(round(landmark_curved[index].y)), \
                              int(round(landmark_curved[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_curved_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                    z + padding_z - 1:z + padding_z + 2] = landmark_curved[index].value

                    # get x, y and z coordinates of straight landmark (rounded to closest integer)
                    x, y, z = int(round(landmark_straight[index].x)), int(round(landmark_straight[index].y)), \
                              int(round(landmark_straight[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_straight_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                    z + padding_z - 1:z + padding_z + 2] = landmark_straight[index].value

                # Write NIFTI volumes
                sct.printv('\nWrite NIFTI volumes...', verbose)
                hdr.set_data_dtype('uint32')  # set imagetype to uint8 #TODO: maybe use int32
                img = Nifti1Image(data_curved_landmarks, None, hdr)
                save(img, 'tmp.landmarks_curved.nii.gz')
                sct.printv('.. File created: tmp.landmarks_curved.nii.gz', verbose)
                img = Nifti1Image(data_straight_landmarks, None, hdr)
                save(img, 'tmp.landmarks_straight.nii.gz')
                sct.printv('.. File created: tmp.landmarks_straight.nii.gz', verbose)
            else:
                # Create volumes containing curved and straight landmarks
                data_curved_landmarks = data * 0
                data_straight_landmarks = data * 0

                # Loop across cross index
                for index in range(0, len(landmark_curved)):
                    x, y, z = int(round(landmark_curved[index].x)), int(round(landmark_curved[index].y)), \
                              int(round(landmark_curved[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_curved_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                                          z + padding_z - 1:z + padding_z + 2] = landmark_curved[index].value

                    # get x, y and z coordinates of straight landmark (rounded to closest integer)
                    x, y, z = int(round(landmark_straight[index].x)), int(round(landmark_straight[index].y)), \
                              int(round(landmark_straight[index].z))

                    # attribute landmark_value to the voxel and its neighbours
                    data_straight_landmarks[x + padding_x - 1:x + padding_x + 2, y + padding_y - 1:y + padding_y + 2,
                                            z + padding_z - 1:z + padding_z + 2] = landmark_straight[index].value

                # Write NIFTI volumes
                sct.printv('\nWrite NIFTI volumes...', verbose)
                hdr.set_data_dtype('uint32')  # set imagetype to uint8 #TODO: maybe use int32
                img = Nifti1Image(data_curved_landmarks, None, hdr)
                save(img, 'tmp.landmarks_curved.nii.gz')
                sct.printv('.. File created: tmp.landmarks_curved.nii.gz', verbose)
                img = Nifti1Image(data_straight_landmarks, None, hdr)
                save(img, 'tmp.landmarks_straight.nii.gz')
                sct.printv('.. File created: tmp.landmarks_straight.nii.gz', verbose)

                # Estimate deformation field by pairing landmarks
                # ==========================================================================================
                # convert landmarks to INT
                sct.printv('\nConvert landmarks to INT...', verbose)
                sct.run('isct_c3d tmp.landmarks_straight.nii.gz -type int -o tmp.landmarks_straight.nii.gz', verbose)
                sct.run('isct_c3d tmp.landmarks_curved.nii.gz -type int -o tmp.landmarks_curved.nii.gz', verbose)

                # This stands to avoid overlapping between landmarks
                # TODO: do symmetric removal
                sct.printv('\nMake sure all labels between landmark_straight and landmark_curved match 1...', verbose)
                label_process_straight = ProcessLabels(fname_label="tmp.landmarks_straight.nii.gz",
                                                       fname_output=["tmp.landmarks_straight.nii.gz",
                                                                     "tmp.landmarks_curved.nii.gz"],
                                                       fname_ref="tmp.landmarks_curved.nii.gz",
                                                       verbose=verbose)
                label_process_straight.process('remove-symm')

                # Estimate rigid transformation
                sct.printv('\nEstimate rigid transformation between paired landmarks...', verbose)
                sct.run('isct_ANTSUseLandmarkImagesToGetAffineTransform tmp.landmarks_straight.nii.gz '
                        'tmp.landmarks_curved.nii.gz rigid tmp.curve2straight_rigid.txt', verbose)

                # Apply rigid transformation
                sct.printv('\nApply rigid transformation to curved landmarks...', verbose)
                Transform(input_filename="tmp.landmarks_curved.nii.gz", fname_dest="tmp.landmarks_curved_rigid.nii.gz",
                          output_filename="tmp.landmarks_straight.nii.gz", warp="tmp.curve2straight_rigid.txt",
                          interp="nn", verbose=verbose).apply()

            if verbose == 2:
                from mpl_toolkits.mplot3d import Axes3D
                import matplotlib.pyplot as plt

                fig = plt.figure()
                ax = Axes3D(fig)
                # plt_centerline, = ax.plot(x_centerline_fit, y_centerline_fit, z_centerline, zdir='z')
                plt_landmarks_curved, = ax.plot([coord.x for coord in landmark_curved],
                                                [coord.y for coord in landmark_curved],
                                                [coord.z for coord in landmark_curved],
                                                'b.', markersize=3)
                for coord in landmark_curved:
                    ax.text(coord.x, coord.y, coord.z, str(coord.value), color='blue')
                plt_landmarks_straight, = ax.plot([coord.x for coord in landmark_straight],
                                                  [coord.y for coord in landmark_straight],
                                                  [coord.z for coord in landmark_straight],
                                                  'r.', markersize=3)
                for coord in landmark_straight:
                    ax.text(coord.x, coord.y, coord.z, str(coord.value), color='red')
                if self.algo_landmark_rigid is not None and self.algo_landmark_rigid != 'None':
                    plt_landmarks_curved_reg, = ax.plot([coord.x for coord in landmark_curved_rigid],
                                                        [coord.y for coord in landmark_curved_rigid],
                                                        [coord.z for coord in landmark_curved_rigid],
                                                        'g.', markersize=3)
                    plt.legend([plt_landmarks_curved, plt_landmarks_straight, plt_landmarks_curved_reg],
                               ['Landmarks curved', 'Landmarks straight', 'Landmarks registered'])
                else:
                    plt.legend([plt_landmarks_curved, plt_landmarks_straight],
                               ['Landmarks curved', 'Landmarks straight'])

                ax.set_xlabel('x')
                ax.set_ylabel('y')
                ax.set_zlabel('z')
                ax.set_aspect('equal')
                plt.show()

            if (self.use_continuous_labels == 1 and self.algo_landmark_rigid is not None and
                    self.algo_landmark_rigid != "None") or self.use_continuous_labels == '1':
                if self.algo_landmark_rigid != 'euler':
                    landmark_curved_rigid, landmark_straight = ProcessLabels.remove_label_coord(landmark_curved_rigid,
                                                                                                landmark_straight,
                                                                                                symmetry=True)

                    # Writing landmark curve in text file
                    landmark_straight_file = open("LandmarksRealStraight.txt", "w+")
                    for i in landmark_straight:
                        landmark_straight_file.write(
                            str(i.x + padding_x) + "," + str(i.y + padding_y) + "," + str(i.z + padding_z) + "\n")
                    landmark_straight_file.close()

                    # Writing landmark curve in text file
                    landmark_curved_file = open("LandmarksRealCurve.txt", "w+")
                    for i in landmark_curved_rigid:
                        landmark_curved_file.write(
                            str(i.x + padding_x) + "," + str(i.y + padding_y) + "," + str(i.z + padding_z) + "\n")
                    landmark_curved_file.close()

                    # Estimate b-spline transformation curve --> straight
                    sct.printv("\nEstimate b-spline transformation: curve --> straight...", verbose)
                    sct.run("isct_ANTSUseLandmarkImagesWithTextFileToGetBSplineDisplacementField "
                            "tmp.landmarks_straight.nii.gz tmp.landmarks_curved_rigid.nii.gz "
                            "tmp.warp_curve2straight.nii.gz " + self.bspline_meshsize + " " +
                            self.bspline_numberOfLevels + " LandmarksRealCurve.txt LandmarksRealStraight.txt " +
                            self.bspline_order + " 0", verbose)
                else:
                    landmark_curved, landmark_straight = ProcessLabels.remove_label_coord(landmark_curved,
                                                                                          landmark_straight,
                                                                                          symmetry=True)

                    # Writing landmark curve in text file
                    landmark_straight_file = open("LandmarksRealStraight.txt", "w+")
                    for i in landmark_straight:
                        landmark_straight_file.write(
                            str(i.x + padding_x) + "," +
                            str(i.y + padding_y) + "," +
                            str(i.z + padding_z) + "\n")
                    landmark_straight_file.close()

                    # Writing landmark curve in text file
                    landmark_curved_file = open("LandmarksRealCurve.txt", "w+")
                    for i in landmark_curved:
                        landmark_curved_file.write(
                            str(i.x + padding_x) + "," +
                            str(i.y + padding_y) + "," +
                            str(i.z + padding_z) + "\n")
                    landmark_curved_file.close()

                    # Estimate b-spline transformation curve --> straight
                    sct.printv("\nEstimate b-spline transformation: curve --> straight...", verbose)
                    status, output = sct.run('isct_ANTSLandmarksBSplineTransform '
                                             'tmp.landmarks_straight.nii.gz '
                                             'tmp.landmarks_curved.nii.gz '
                                             'tmp.curve2straight_rigid.txt '
                                             'tmp.warp_curve2straight.nii.gz ' +
                                             self.bspline_meshsize + ' ' +
                                             self.bspline_numberOfLevels + ' '
                                             'LandmarksRealCurve.txt '
                                             'LandmarksRealStraight.txt ' +
                                             self.bspline_order +
                                             ' 0',
                                             verbose=verbose)
            else:
                # This stands to avoid overlapping between landmarks
                sct.printv("\nMake sure all labels between landmark_straight and landmark_curved match 2...", verbose)
                label_process = ProcessLabels(fname_label="tmp.landmarks_curved_rigid.nii.gz",
                                              fname_output=["tmp.landmarks_curved_rigid.nii.gz",
                                                            "tmp.landmarks_straight.nii.gz"],
                                              fname_ref="tmp.landmarks_straight.nii.gz", verbose=verbose)
                label_process.process("remove-symm")

                # Estimate b-spline transformation curve --> straight
                sct.printv("\nEstimate b-spline transformation: curve --> straight...", verbose)
                sct.run("isct_ANTSUseLandmarkImagesToGetBSplineDisplacementField tmp.landmarks_straight.nii.gz "
                        "tmp.landmarks_curved_rigid.nii.gz tmp.warp_curve2straight.nii.gz " + self.bspline_meshsize
                        + " " + self.bspline_numberOfLevels + " " + self.bspline_order + " 0", verbose)

            # remove padding for straight labels
            if crop == 1:
                ImageCropper(input_file="tmp.landmarks_straight.nii.gz",
                             output_file="tmp.landmarks_straight_crop.nii.gz", dim=[0, 1, 2], bmax=True,
                             verbose=verbose).crop()
                pass
            else:
                sct.run("cp tmp.landmarks_straight.nii.gz tmp.landmarks_straight_crop.nii.gz", verbose)

            # Concatenate rigid and non-linear transformations...
            sct.printv("\nConcatenate rigid and non-linear transformations...", verbose)
            # TODO: !!! DO NOT USE sct.run HERE BECAUSE isct_ComposeMultiTransform OUTPUTS A NON-NULL STATUS !!!
            if self.algo_landmark_rigid == 'rigid-decomposed':
                cmd = "sct_concat_transfo -w tmp.curve2straight_rigid1.txt,tmp.curve2straight_rigid2.txt," \
                      "tmp.warp_curve2straight.nii.gz -d tmp.landmarks_straight_crop.nii.gz " \
                      "-o tmp.curve2straight.nii.gz"
            else:
                cmd = "isct_ComposeMultiTransform 3 tmp.curve2straight.nii.gz -R "+str(file_anat+ext_anat)+" " \
                      "tmp.warp_curve2straight.nii.gz tmp.curve2straight_rigid.txt"
            sct.run(cmd, self.verbose)

            # Estimate b-spline transformation straight --> curve
            # TODO: invert warping field instead of estimating a new one
            sct.printv("\nEstimate b-spline transformation: straight --> curve...", verbose)
            if (self.use_continuous_labels == 1 and self.algo_landmark_rigid is not None
                    and self.algo_landmark_rigid != "None") or self.use_continuous_labels == "1":
                if self.algo_landmark_rigid != 'euler':
                    sct.run("isct_ANTSUseLandmarkImagesWithTextFileToGetBSplineDisplacementField "
                            "tmp.landmarks_curved_rigid.nii.gz tmp.landmarks_straight.nii.gz "
                            "tmp.warp_straight2curve.nii.gz " + self.bspline_meshsize + " " +
                            self.bspline_numberOfLevels + " LandmarksRealStraight.txt LandmarksRealCurve.txt " +
                            self.bspline_order + " 0", verbose)
                else:
                    # Estimate b-spline transformation curve --> straight
                    status, output = sct.run('isct_ANTSLandmarksBSplineTransform '
                                             'tmp.landmarks_curved.nii.gz '
                                             'tmp.landmarks_straight.nii.gz '
                                             'tmp.straight2curve_rigid.txt '
                                             'tmp.warp_straight2curve.nii.gz ' +
                                             self.bspline_meshsize + ' ' +
                                             self.bspline_numberOfLevels + ' ' +
                                             'LandmarksRealStraight.txt '
                                             'LandmarksRealCurve.txt ' +
                                             self.bspline_order +
                                             ' 0',
                                             verbose=verbose)
            else:
                sct.run("isct_ANTSUseLandmarkImagesToGetBSplineDisplacementField "
                        "tmp.landmarks_curved_rigid.nii.gz tmp.landmarks_straight.nii.gz "
                        "tmp.warp_straight2curve.nii.gz " + self.bspline_meshsize + " " + self.bspline_numberOfLevels +
                        " " + self.bspline_order + " 0", verbose)

            # Concatenate rigid and non-linear transformations...
            sct.printv("\nConcatenate rigid and non-linear transformations...", verbose)
            if self.algo_landmark_rigid == "rigid-decomposed":
                cmd = "sct_concat_transfo -w tmp.warp_straight2curve.nii.gz,-tmp.curve2straight_rigid2.txt," \
                      "-tmp.curve2straight_rigid1.txt -d " + file_anat + ext_anat + " -o tmp.straight2curve.nii.gz"
            else:
                if self.algo_landmark_rigid == 'euler':
                    #cmd = "isct_ComposeMultiTransform 3 tmp.straight2curve.nii.gz -R " + file_anat + ext_anat + \
                    #      " tmp.warp_straight2curve.nii.gz tmp.straight2curve_rigid.txt"
                    cmd = "isct_ComposeMultiTransform 3 tmp.straight2curve.nii.gz -R " + file_anat + ext_anat + \
                          " tmp.warp_straight2curve.nii.gz tmp.straight2curve_rigid.txt"
                else:
                    cmd = "isct_ComposeMultiTransform 3 tmp.straight2curve.nii.gz -R " + file_anat + ext_anat + \
                          " -i tmp.curve2straight_rigid.txt tmp.warp_straight2curve.nii.gz"
            sct.run(cmd, self.verbose)

            # Apply transformation to input image
            sct.printv('\nApply transformation to input image...', verbose)
            Transform(input_filename=str(file_anat+ext_anat), fname_dest=str(file_anat+ext_anat),
                      output_filename="tmp.anat_rigid_warp.nii.gz", interp=interpolation_warp,
                      warp="tmp.curve2straight.nii.gz", verbose=verbose).apply()

            # compute the error between the straightened centerline/segmentation and the central vertical line.
            # Ideally, the error should be zero.
            # Apply deformation to input image
            sct.printv('\nApply transformation to centerline image...', verbose)
            Transform(input_filename=fname_centerline_orient, fname_dest="tmp.landmarks_straight_crop.nii.gz",
                      output_filename="tmp.centerline_straight.nii.gz", interp="nn",
                      warp="tmp.curve2straight.nii.gz", verbose=verbose).apply()
            from msct_image import Image
            file_centerline_straight = Image('tmp.centerline_straight.nii.gz', verbose=verbose)
            coordinates_centerline = file_centerline_straight.getNonZeroCoordinates(sorting='z')
            mean_coord = []
            from numpy import mean
            for z in range(coordinates_centerline[0].z, coordinates_centerline[-1].z):
                temp_mean = [coord.value for coord in coordinates_centerline if coord.z == z]
                if temp_mean:
                    mean_value = mean(temp_mean)
                    mean_coord.append(mean([[coord.x * coord.value / mean_value, coord.y * coord.value / mean_value]
                                            for coord in coordinates_centerline if coord.z == z], axis=0))

            # compute error between the input data and the nurbs
            from math import sqrt
            x0 = file_centerline_straight.data.shape[0]/2.0
            y0 = file_centerline_straight.data.shape[1]/2.0
            count_mean = 0
            for coord_z in mean_coord:
                if not isnan(sum(coord_z)):
                    dist = ((x0-coord_z[0])*px)**2 + ((y0-coord_z[1])*py)**2
                    self.mse_straightening += dist
                    dist = sqrt(dist)
                    if dist > self.max_distance_straightening:
                        self.max_distance_straightening = dist
                    count_mean += 1
            self.mse_straightening = sqrt(self.mse_straightening/float(count_mean))

        except Exception as e:
            sct.printv('WARNING: Exception during Straightening:', 1, 'warning')
            print 'Error on line {}'.format(sys.exc_info()[-1].tb_lineno)
            print e

        os.chdir('..')

        # Generate output file (in current folder)
        # TODO: do not uncompress the warping field, it is too time consuming!
        sct.printv("\nGenerate output file (in current folder)...", verbose)
        sct.generate_output_file(path_tmp + "/tmp.curve2straight.nii.gz", "warp_curve2straight.nii.gz", verbose)
        sct.generate_output_file(path_tmp + "/tmp.straight2curve.nii.gz", "warp_straight2curve.nii.gz", verbose)
        if fname_output == '':
            fname_straight = sct.generate_output_file(path_tmp + "/tmp.anat_rigid_warp.nii.gz",
                                                      file_anat + "_straight" + ext_anat, verbose)
        else:
            fname_straight = sct.generate_output_file(path_tmp+'/tmp.anat_rigid_warp.nii.gz',
                                                      fname_output, verbose)  # straightened anatomic

        # Remove temporary files
        if remove_temp_files:
            sct.printv("\nRemove temporary files...", verbose)
            sct.run("rm -rf " + path_tmp, verbose)

        sct.printv('\nDone!\n', verbose)

        sct.printv("Maximum x-y error = " + str(round(self.max_distance_straightening, 2)) + " mm", verbose, "bold")
        sct.printv("Accuracy of straightening (MSE) = " + str(round(self.mse_straightening, 2)) +
                   " mm", verbose, "bold")
        # display elapsed time
        elapsed_time = time.time() - start_time
        sct.printv("\nFinished! Elapsed time: " + str(int(round(elapsed_time))) + "s", verbose)
        sct.printv("\nTo view results, type:", verbose)
        sct.printv("fslview " + fname_straight + " &\n", verbose, "info")
dataset/ETHPy150Open neuropoly/spinalcordtoolbox/dev/straightening/sct_straighten_spinalcord_LargeFOVOutput.py/SpinalCordStraightener.straighten