Как я могу выполнить регистрацию 3D-изображения с помощью SimpleElastix из Python?

Question

Я пытаюсь зарегистрировать два 3d тома. Попытка этого может быть найдена здесь. Сначала код генерирует два разных тома, каждый из которых содержит ровно одну сферу радиуса 4. Затем я пытаюсь выровнять их, используя карту параметров перевода по умолчанию. Однако, как видно из последней строки (и из графиков при локальном запуске), результирующий том вообще не выровнен с фиксированным объемом. При попытке выполнить ту же процедуру в 2d на этот раз результирующее изображение выглядит правильно выровненным с фиксированным изображением, как можно видеть здесь. Я использую API SimpleElastix неправильно? Я просмотрел репозиторий SimpleElastix на Github, но не смог найти примеров регистрации 3D-изображений (по крайней мере, не используя тома, сгенерированные в Python, а затем преобразовав их в изображения ITK).

Код из 3d примера:

vol1 = np.zeros((50, 50, 50))
for x in range(vol.shape[0]):
    for y in range(vol.shape[1]):
        for z in range(vol.shape[2]):
            vol1[x, y, z] = np.linalg.norm(np.subtract([x, y, z], [5, 3, 2])) < 4

vol2 = np.zeros((50, 50, 50))
for x in range(vol.shape[0]):
    for y in range(vol.shape[1]):
        for z in range(vol.shape[2]):
            vol1[x, y, z] = np.linalg.norm(np.subtract([x, y, z], [20, 30, 10])) < 4

img_a = sitk.GetImageFromArray(vol1)
img_b = sitk.GetImageFromArray(vol2)

parameterMap = sitk.GetDefaultParameterMap('translation')

itk_filter = sitk.ElastixImageFilter()
itk_filter.LogToConsoleOn()
itk_filter.SetFixedImage(img_a)
itk_filter.SetMovingImage(img_b)
itk_filter.SetParameterMap(parameterMap)
itk_filter.Execute()

result_vol = sitk.GetArrayFromImage(itk_filter.GetResultImage())

np.max(np.abs(vol1 - result_vol))

Код из 2-го примера:

vol1 = np.zeros((50, 50))
for x in range(vol1.shape[0]):
    for y in range(vol1.shape[1]):
        vol1[x, y] = np.linalg.norm(np.subtract([x, y], [20, 20])) < 4

vol2 = np.zeros((50, 50))
for x in range(vol2.shape[0]):
    for y in range(vol2.shape[1]):
        vol2[x, y] = np.linalg.norm(np.subtract([x, y], [4, 5])) < 4

img_a = sitk.GetImageFromArray(vol1)
img_b = sitk.GetImageFromArray(vol2)

parameterMap = sitk.GetDefaultParameterMap('translation')

itk_filter = sitk.ElastixImageFilter()
itk_filter.LogToConsoleOn()
itk_filter.SetFixedImage(img_a)
itk_filter.SetMovingImage(img_b)
itk_filter.SetParameterMap(parameterMap)
itk_filter.Execute()

result_vol = sitk.GetArrayFromImage(itk_filter.GetResultImage())

np.max(np.abs(vol1 - result_vol))

Answer 1

Возможно, регистрация не удалась, потому что сферы не перекрываются. Вы можете попытаться сгладить изображения немного больше, используя, например, sitk.DiscreteGaussian(), Однако в этом искусственном примере весьма вероятно, что простое перемещение сфер ближе друг к другу приведет к хорошей регистрации.

Обратите внимание, что сложно регистрировать двоичные изображения, поскольку в изображениях такого типа очень мало информации о градиенте: градиенты отличны от нуля только на границе и равны нулю везде, потому что области имеют плоскую интенсивность.