Python Image.getCoordinatesAveragedByValue方法代码示例

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import getCoordinatesAveragedByValue [as 别名]
def register_landmarks(fname_src, fname_dest, dof, fname_affine="affine.txt", verbose=1, path_qc="./"):
    """
    Register two NIFTI volumes containing landmarks
    :param fname_src: fname of source landmarks
    :param fname_dest: fname of destination landmarks
    :param dof: degree of freedom. Separate with "_". Example: Tx_Ty_Tz_Rx_Ry_Sz
    :param fname_affine: output affine transformation
    :param verbose: 0, 1, 2
    :return:
    """
    from msct_image import Image

    # open src label
    im_src = Image(fname_src)
    # coord_src = im_src.getNonZeroCoordinates(sorting='value')  # landmarks are sorted by value
    coord_src = im_src.getCoordinatesAveragedByValue()  # landmarks are sorted by value
    # open dest labels
    im_dest = Image(fname_dest)
    # coord_dest = im_dest.getNonZeroCoordinates(sorting='value')
    coord_dest = im_dest.getCoordinatesAveragedByValue()
    # Reorganize landmarks

    points_src, points_dest = [], []
    for coord in coord_src:
        point_src = im_src.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        # convert NIFTI to ITK world coordinate
        # points_src.append([point_src[0][0], point_src[0][1], point_src[0][2]])
        points_src.append([-point_src[0][0], -point_src[0][1], point_src[0][2]])
    for coord in coord_dest:
        point_dest = im_dest.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        # convert NIFTI to ITK world coordinate
        # points_dest.append([point_dest[0][0], point_dest[0][1], point_dest[0][2]])
        points_dest.append([-point_dest[0][0], -point_dest[0][1], point_dest[0][2]])

    # display
    sct.printv("Labels src: " + str(points_src), verbose)
    sct.printv("Labels dest: " + str(points_dest), verbose)
    sct.printv("Degrees of freedom (dof): " + dof, verbose)

    if len(coord_src) != len(coord_dest):
        raise Exception(
            "Error: number of source and destination landmarks are not the same, so landmarks cannot be paired."
        )

    # estimate transformation
    # N.B. points_src and points_dest are inverted below, because ITK uses inverted transformation matrices, i.e., src->dest is defined in dest instead of src.
    # (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = getRigidTransformFromLandmarks(points_dest, points_src, constraints=dof, verbose=verbose, path_qc=path_qc)
    (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = getRigidTransformFromLandmarks(
        points_src, points_dest, constraints=dof, verbose=verbose, path_qc=path_qc
    )
    # writing rigid transformation file
    # N.B. x and y dimensions have a negative sign to ensure compatibility between Python and ITK transfo
    text_file = open(fname_affine, "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()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import getCoordinatesAveragedByValue [as 别名]

#.........这里部分代码省略.........

    # Label preparation:
    # --------------------------------------------------------------------------------
    # Remove unused label on template. Keep only label present in the input label image
    sct.printv('\nRemove unused label on template. Keep only label present in the input label image...', verbose)
    sct.run('sct_label_utils -p remove -i '+ftmp_template_label+' -o '+ftmp_template_label+' -r '+ftmp_label)

    # Dilating the input label so they can be straighten without losing them
    sct.printv('\nDilating input labels using 3vox ball radius')
    sct.run('sct_maths -i '+ftmp_label+' -o '+add_suffix(ftmp_label, '_dilate')+' -dilate 3')
    ftmp_label = add_suffix(ftmp_label, '_dilate')

    # Apply straightening to labels
    sct.printv('\nApply straightening to labels...', verbose)
    sct.run('sct_apply_transfo -i '+ftmp_label+' -o '+add_suffix(ftmp_label, '_straight')+' -d '+add_suffix(ftmp_seg, '_straight')+' -w warp_curve2straight.nii.gz -x nn')
    ftmp_label = add_suffix(ftmp_label, '_straight')

    # Create crosses for the template labels and get coordinates
    sct.printv('\nCreate a 15 mm cross for the template labels...', verbose)
    template_image = Image(ftmp_template_label)
    coordinates_input = template_image.getNonZeroCoordinates(sorting='value')
    # jcohenadad, issue #628 <<<<<
    # landmark_template = ProcessLabels.get_crosses_coordinates(coordinates_input, gapxy=15)
    landmark_template = coordinates_input
    # >>>>>
    if verbose == 2:
        # TODO: assign cross to image before saving
        template_image.setFileName(add_suffix(ftmp_template_label, '_cross'))
        template_image.save(type='minimize_int')

    # Create crosses for the input labels into straight space and get coordinates
    sct.printv('\nCreate a 15 mm cross for the input labels...', verbose)
    label_straight_image = Image(ftmp_label)
    coordinates_input = label_straight_image.getCoordinatesAveragedByValue()  # landmarks are sorted by value
    # jcohenadad, issue #628 <<<<<
    # landmark_straight = ProcessLabels.get_crosses_coordinates(coordinates_input, gapxy=15)
    landmark_straight = coordinates_input
    # >>>>>
    if verbose == 2:
        # TODO: assign cross to image before saving
        label_straight_image.setFileName(add_suffix(ftmp_label, '_cross'))
        label_straight_image.save(type='minimize_int')

    # Reorganize landmarks
    points_fixed, points_moving = [], []
    for coord in landmark_straight:
        point_straight = label_straight_image.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_moving.append([point_straight[0][0], point_straight[0][1], point_straight[0][2]])

    for coord in landmark_template:
        point_template = template_image.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_fixed.append([point_template[0][0], point_template[0][1], point_template[0][2]])

    # Register curved landmarks on straight landmarks based on python implementation
    sct.printv('\nComputing rigid transformation (algo=translation-scaling-z) ...', verbose)

    import msct_register_landmarks
    # for some reason, the moving and fixed points are inverted between ITK transform and our python-based transform.
    # and for another unknown reason, x and y dimensions have a negative sign (at least for translation and center of rotation).
    if verbose == 2:
        show_transfo = True
    else:
        show_transfo = False
    (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = msct_register_landmarks.getRigidTransformFromLandmarks(points_moving, points_fixed, constraints='translation-scaling-z', show=show_transfo)
    # writing rigid transformation file
    text_file = open("straight2templateAffine.txt", "w")

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import getCoordinatesAveragedByValue [as 别名]

#.........这里部分代码省略.........
            # register src --> dest
            # TODO: display param for debugging
            warp_forward_out, warp_inverse_out = register(src, dest, paramreg, param, str(i_step))
            warp_forward.append(warp_forward_out)
            warp_inverse.append(warp_inverse_out)

        # Concatenate transformations:
        sct.printv('\nConcatenate transformations: anat --> template...', verbose)
        sct.run('sct_concat_transfo -w warp_curve2straightAffine.nii.gz,'+','.join(warp_forward)+' -d template.nii -o warp_anat2template.nii.gz', verbose)
        # sct.run('sct_concat_transfo -w warp_curve2straight.nii.gz,straight2templateAffine.txt,'+','.join(warp_forward)+' -d template.nii -o warp_anat2template.nii.gz', verbose)
        sct.printv('\nConcatenate transformations: template --> anat...', verbose)
        warp_inverse.reverse()
        sct.run('sct_concat_transfo -w '+','.join(warp_inverse)+',-straight2templateAffine.txt,warp_straight2curve.nii.gz -d data.nii -o warp_template2anat.nii.gz', verbose)

    # register template->subject
    elif ref == 'subject':

        # Change orientation of input images to RPI
        sct.printv('\nChange orientation of input images to RPI...', verbose)
        sct.run('sct_image -i ' + ftmp_data + ' -setorient RPI -o ' + add_suffix(ftmp_data, '_rpi'))
        ftmp_data = add_suffix(ftmp_data, '_rpi')
        sct.run('sct_image -i ' + ftmp_seg + ' -setorient RPI -o ' + add_suffix(ftmp_seg, '_rpi'))
        ftmp_seg = add_suffix(ftmp_seg, '_rpi')
        sct.run('sct_image -i ' + ftmp_label + ' -setorient RPI -o ' + add_suffix(ftmp_label, '_rpi'))
        ftmp_label = add_suffix(ftmp_label, '_rpi')

        # Remove unused label on template. Keep only label present in the input label image
        sct.printv('\nRemove unused label on template. Keep only label present in the input label image...', verbose)
        sct.run('sct_label_utils -i '+ftmp_template_label+' -o '+ftmp_template_label+' -remove '+ftmp_label)

        # Add one label because at least 3 orthogonal labels are required to estimate an affine transformation. This new label is added at the level of the upper most label (lowest value), at 1cm to the right.
        for i_file in [ftmp_label, ftmp_template_label]:
            im_label = Image(i_file)
            coord_label = im_label.getCoordinatesAveragedByValue()  # N.B. landmarks are sorted by value
            # Create new label
            from copy import deepcopy
            new_label = deepcopy(coord_label[0])
            # move it 5mm to the left (orientation is RAS)
            nx, ny, nz, nt, px, py, pz, pt = im_label.dim
            new_label.x = round(coord_label[0].x + 5.0 / px)
            # assign value 99
            new_label.value = 99
            # Add to existing image
            im_label.data[new_label.x, new_label.y, new_label.z] = new_label.value
            # Overwrite label file
            # im_label.setFileName('label_rpi_modif.nii.gz')
            im_label.save()

        # Bring template to subject space using landmark-based transformation
        sct.printv('\nEstimate transformation for step #0...', verbose)
        from msct_register_landmarks import register_landmarks
        warp_forward = ['template2subjectAffine.txt']
        warp_inverse = ['-template2subjectAffine.txt']
        try:
            register_landmarks(ftmp_template_label, ftmp_label, paramreg.steps['0'].dof, fname_affine=warp_forward[0], verbose=verbose, path_qc=param.path_qc)
        except Exception:
            sct.printv('ERROR: input labels do not seem to be at the right place. Please check the position of the labels. See documentation for more details: https://sourceforge.net/p/spinalcordtoolbox/wiki/create_labels/', verbose=verbose, type='error')

        # loop across registration steps
        for i_step in range(1, len(paramreg.steps)):
            sct.printv('\nEstimate transformation for step #'+str(i_step)+'...', verbose)
            # identify which is the src and dest
            if paramreg.steps[str(i_step)].type == 'im':
                src = ftmp_template
                dest = ftmp_data
                interp_step = 'linear'