Как преобразовать размер анизотропного вокселя (CT) в размер изотропного вокселя ⇐ Python

[Anonymous]

Я работаю над кодом для стандартизации размеров вокселей серии 3D Dicom на нескольких машинах (в рамках конвейера, но сейчас мы сосредоточимся на размерах вокселей). Он должен иметь возможность обрабатывать ввод как изотропных, так и анизотропных вокселей, а также сохранять исходный размер объектов, поскольку я буду выполнять измерения.
Код основан на Resample Scalar Volume от Slicer. используя в основном простой ITK, который хорошо работает для всех моих входов dicom. Исходный код написан на C++, но я модифицирую его на Python. Однако мой код привел к некоторым ошибкам. Я тестирую все свои результаты в Slicer, так как буду использовать это программное обеспечение для анализа. Пожалуйста, будьте осторожны, поскольку я новичок в программировании.
Я попробовал два разных подхода, оба из которых зашли в тупик. Буду рад получить совет по любому из них. В обоих подходах фактическое количество размеров вокселей и толщина среза выходных данных уже верны, но внешний вид неверен.
Подход 1: Вывод в виде DICOM. На выходе создаются только файлы .dcm, но по какой-то причине Slicer обнаружил преобразование, хотя файла преобразования нет. Фактический объем CT является правильным, но в результате преобразования объекты были измерены меньше, чем должны в плоскости xy. Я хочу избавиться от эффекта преобразования, так как боюсь, что он перейдет к следующей части конвейера. Это происходит как с изотропными, так и с анизотропными входами. Я уже пытался полностью избавиться от SetTransform из кода, но эффект трансформации все равно сохранялся. Попробовал отрегулировать направление и положение изображения, но не уверен, что это работает.
При активированном преобразовании (объекты меньшего размера в плоскости xy):
При активированном преобразовании (объекты меньшего размера в плоскости xy)< /p>
Преобразование деактивировано (правильный внешний вид):
Преобразование деактивировано (правильный внешний вид)
# 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 * (original_spacing / target_spacing) + 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
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)

# 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. Я подозреваю, что проблема 1 вызвана сохранением в процессе dicom, поэтому при втором подходе я попробовал сохранить в nifti. Это сработало хорошо, если мои входные данные изотропны, никаких странных преобразований, выглядит великолепно. Однако если мои входные данные анизотропны, изображение в плоскости xy выглядит нормально, но оно растянуто по оси Z, перевернуто, с странными чередующимися линиями . Попробовал отрегулировать направление и положение изображения, но не уверен, что это что-то даст, поскольку начало изображения равно 0 (в исходном изображении оно не было 0).
Я не уверен, что происходит с z. -ось, выглядит почти так, как будто она склеена с другим изображением/другой частью изображения, за исключением верхней части изображения (фактически нижней части пациента, поскольку изображение перевернуто). Абсолютный размер по оси Z и размер по высоте почти в 3 раза больше, чем должен быть.
Изображение оси Z:
Изображение оси Z
# 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 * (original_spacing / target_spacing) + 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'))



Подробнее здесь: https://stackoverflow.com/questions/791 ... voxel-size