各向异性CT体素转各向同性的Python代码问题排查

需求为标准化多设备3D DICOM序列的体素尺寸，需兼容各向同性/各向异性输入，且保留物体原始尺寸以用于测量。代码基于Slicer的Resample Scalar Volume模块，将C++代码转为Python版并使用SimpleITK实现，所有输出均在Slicer中测试。

尝试两种方案均遇到瓶颈：输出体素尺寸和层厚虽正确，但视觉显示异常：

• 方案1（输出DICOM序列）：无变换文件，但Slicer检测到变换，导致XY平面物体测量偏小，移除SetTransform后问题仍存在。
• 方案2（输出NIfTI）：各向同性输入时正常，但各向异性输入时Z轴拉伸、倒置且出现异常交替线条，Z轴尺寸偏大近3倍。

方案1：输出DICOM序列代码 #

# Function to resize voxel spacing using SimpleITK
def resize_voxel(itk_image, original_spacing, original_direction, target_spacing):
    original_size = itk_image.GetSize()
    target_size = [
        int(original_size[i] * (original_spacing[i] / target_spacing[i]) + 0.5)
        for i in range(3)
    ]

    resampler = sitk.ResampleImageFilter()
    resampler.SetOutputSpacing(target_spacing)
    resampler.SetSize(target_size)
    resampler.SetOutputDirection(itk_image.GetDirection())
    resampler.SetOutputOrigin(itk_image.GetOrigin())
    
    # Use an explicit identity transform to prevent unintended scaling or rotation
    transform = sitk.Euler3DTransform()
    transform.SetIdentity()
    resampler.SetTransform(transform)
    resampler.SetInterpolator(sitk.sitkBSpline)

    # Adjust to ensure absolute sizes are preserved with high precision
    resampler.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resampler.SetOutputPixelType(itk_image.GetPixelID())
    resampler.SetNumberOfThreads(4)

    resized_image = resampler.Execute(itk_image)
    return resized_image

# Function to save the 3D image as a DICOM series
def save_dicom_series(image_3d, dicom_names, output_folder, target_spacing, original_image_position_patient):
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    # Convert SimpleITK image back to NumPy array
    image_3d_array = sitk.GetArrayFromImage(image_3d)

    for i in range(len(image_3d_array)):
        slice_array = image_3d_array[i]
        output_dicom_path = os.path.join(output_folder, f"{i:03d}.dcm")

        # Load the original slice's DICOM for metadata and anonymize it
        original_dicom = anonymize_dicom(dicom_names[i])

        # Update the pixel data
        original_dicom.PixelData = slice_array.tobytes()
        original_dicom.Rows, original_dicom.Columns = slice_array.shape

        # Update the spacing information in the DICOM metadata
        original_dicom.PixelSpacing = [target_spacing[0], target_spacing[1]]
        original_dicom.SliceThickness = target_spacing[2]
        original_dicom.SpacingBetweenSlices = target_spacing[2]

        # Update Image Position Patient to reflect new slice location
        if original_image_position_patient is not None:
            new_position = original_image_position_patient.copy()
            new_position[2] = round(original_image_position_patient[2] + float(i) * target_spacing[2], 6)
            original_dicom.ImagePositionPatient = new_position

        # Update Slice Location and Instance Number
        original_dicom.InstanceNumber = i + 1
        original_dicom.SliceLocation = i * float(original_dicom.SliceThickness)

        # Debug output for updated metadata
        print(f"Saving slice {i + 1} with SliceThickness: {original_dicom.SliceThickness}, PixelSpacing: {original_dicom.PixelSpacing}, SpacingBetweenSlices: {original_dicom.SpacingBetweenSlices}, ImagePositionPatient: {original_dicom.ImagePositionPatient}")

        original_dicom.save_as(output_dicom_path)

# Main function to process the entire DICOM series
def process_dicom_series(patient_folder, output_folder):
    # Load the DICOM series
    image_3d, dicom_names, spacing, image_position_patient, image_direction = load_dicom_series(patient_folder)

    # Convert numpy array to SimpleITK image for processing
    itk_image = sitk.GetImageFromArray(image_3d)
    itk_image.SetSpacing(spacing)
    itk_image.SetDirection(image_direction)

    # Resize the voxel image to standardized voxel size of 180x180x180 micrometers
    target_spacing = (0.18, 0.18, 0.18)
    resized_voxel_image = resize_voxel(itk_image, spacing, image_direction, target_spacing)
    
    # Save the processed 3D image as a DICOM series
    save_dicom_series(resized_voxel_image, dicom_names, output_folder, target_spacing, image_position_patient)

方案2：输出NIfTI文件代码 #

# Function to resize voxel spacing using SimpleITK
def resize_voxel(itk_image, original_spacing, original_direction, target_spacing):
    original_size = itk_image.GetSize()
    target_size = [
        int(original_size[i] * (original_spacing[i] / target_spacing[i]) + 0.5)
        for i in range(3)
    ]

    resampler = sitk.ResampleImageFilter()
    resampler.SetOutputSpacing(target_spacing)
    resampler.SetSize(target_size)
    resampler.SetOutputDirection(original_direction)
    resampler.SetOutputOrigin(itk_image.GetOrigin())
    
    # Use an explicit identity transform to prevent unintended scaling or rotation
    transform = sitk.Euler3DTransform()
    transform.SetIdentity()
    resampler.SetTransform(transform)
    resampler.SetInterpolator(sitk.sitkBSpline)

    # Adjust to ensure absolute sizes are preserved with high precision
    resampler.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resampler.SetOutputPixelType(itk_image.GetPixelID())
    resampler.SetNumberOfThreads(4)

    resized_image = resampler.Execute(itk_image)
    return resized_image

# Function to save the 3D image as a NIfTI file using SimpleITK
def save_nifti(itk_image, output_path):
    sitk.WriteImage(itk_image, output_path)
    print(f"NIfTI volume successfully saved to {output_path}")

# Main function to process the entire DICOM series
def process_dicom_series(patient_folder, output_folder):
    # Load the DICOM series
    image_3d, dicom_names, spacing, image_position_patient, image_direction = load_dicom_series(patient_folder)

    # Convert numpy array to SimpleITK image for processing
    itk_image = sitk.GetImageFromArray(image_3d)
    itk_image.SetSpacing(spacing)
    itk_image.SetDirection(image_direction)

    # Resize the voxel image to standardized voxel size of 250x250x250 micrometers
    target_spacing = (0.25, 0.25, 0.25)
    resized_voxel_image = resize_voxel(itk_image, spacing, image_direction, target_spacing)
    
    # Save the processed 3D image as a DICOM series using SimpleITK
    save_nifti(resized_voxel_image, os.path.join(output_folder, 'output_volume.nii'))

内容的提问来源于stack exchange，提问作者Andrew S