Продолжаем разбираться со особенностями МРТ-данных для обучения нейронных сетей. Содержание и первые части цикла статей здесь.

6. Некоторые библиотеки для работы с МРТ-изображениями

Основной акцент при рассмотрении всех аспектов работы с МРТ-изображениями делается на данных в формате DICOM, которые создаются непосредственно при проведении исследования на МР-томографе, а, значит, не содержат неизвестных искажений от применения различных преобразований.

Кроме этого формата при работе с МРТ-данными часто используются и другие – как промежуточные или более удобные для тех или иных задач или областей исследований – NRRD, NIFTI, MHA и др. Для каждого формата существуют свои библиотеки для Python, равно как и многие библиотеки поддерживают по несколько форматов. Здесь и сейчас не стоит цель рассмотреть все из них – остановимся только на четырех. И только на некоторых их особенностях – все остальное можно прочитать в документации.

6.1. SimpleITK

ITK (Insight ToolKit) – это кроссплатформенная библиотека с открытым исходным кодом, которая содержит обширный набор программных инструментов для анализа изображений. Предназначена для обработки, сегментации и регистрации научных изображений в двух, трех или более измерениях.

SimpleITK – это упрощенный программный интерфейс для базовых алгоритмов и структур данных ITK. Поддерживает несколько языков программирования, включая C++, Python, R, Java, C#, Lua, Ruby и TCL, а также более 15-ти форматов данных, включая DICOM и NRRD. Разработан сообществом ITK, в том числе, для биомедицинских наук.

Сайт ITK
Сайт SimpleITK
Документация: https://simpleitk.readthedocs.io/, https://simpleitk.org/doxygen/latest/html/
Исходники

Дальнейшие примеры кода будут использовать следующий вид импорта данной библиотеки:

import SimpleITK as sitk

В SimpleITK реализованы 2 API – объектный и процедурный. При этом часто функции просто дублируют вызовы объектного API соответствующих классов, обеспечивая их упрощенный вызов.

Например, считывание одиночного DICOM-файла может выполняться в процедурном стиле:

image = sitk.ReadImage(file_path)

который является аналогом следующего объектного кода:

reader = sitk.ImageFileReader()
reader.SetFileNames(file_name)
image = reader.Execute()

Многие другие задачи решаются аналогично, например, считывание серии из набора файлов:

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(list_of_file_names)
image_series = series_reader.Execute()

Но, например, загрузка тэгов DICOM-файла хотя и выполняется в объектном API, не возвращает данные, а сохраняет их внутри самого ридера. При этом далеко не все тэги становятся доступными:

reader = sitk.ImageFileReader()
reader.SetFileName(file_name)
reader.ReadImageInformation()

all_tags = reader.GetMetaDataKeys()
if reader.HasMetaDataKey(key):
   value = reader.GetMetaData(key)

SimpleITK базируется на двух фундаментальных понятиях – изображение (Image) и пространственное преобразование (Transform). Преобразования применяются к изображениям с помощью процедуры ресэмплинга (Resample).

Image

Изображение в SimpleITK – это набор точек на прямоугольной сетке, занимающий физическую область в пространстве глобальной системы координат – ровно так, как это определяется для вокселей. Это существенно отличается от многих других библиотек, которые рассматривают изображение как массив квадратных пикселей. Как следствие, в SimpleITK важную роль играет интервал между пикселями, который может быть анизотропен. При работе с изображениями SimpleITK следует помнить, что координаты каждого пикселя соответствуют его центру (точке сетки), а сам пиксель занимает пространство вокруг этой точки в соответствии с его реальными размерами.

Не менее значима и привязка изображения к физическому пространству, задаваемому глобальной системой координат, единицами измерения которой обычно подразумеваются (но нигде не определяются) миллиметры, а начало координат зафиксировано в пространстве относительно томографа или пациента. При этом, с каждым изображением связана своя система координат, определяемая векторами первой строки\первого столбца\и т.п. данных – она может быть не ортогональна.

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

Значением каждого пикселя для МРТ-данных является одно целое или вещественное число. Однако SimpleITK поддерживает и векторные значения, например, для RGB, RGBA и дру��их многоканальных изображений, а также комплексные значения. Последние могут вызывать проблемы и при работе с МРТ-изображениями, если тип загружаемых данных явно не указан.

Получение геометрических параметров для объекта класса Image выполняется методами этого класса:

image.GetOrigin(), image.GetSpacing(), image.GetDirection()
image.GetSize(), image.GetPixel((0, 0))
image.TransformIndexToPhysicalPoint((0, 0)) и др.

Изменение таких параметров лучше делать не вручную, а через применение пространственных преобразований.

Иногда требуется изменение только лишь значений пикселей с помощью алгоритма, который не реализован в библиотеке или неудобен\нестабилен. Также может возникнуть необходимость использовать только лишь значения как изображения в других библиотеках, например, при визуализации или подготовке батчей при обучении нейронных сетей. Для таких задач можно использовать конвертацию данных из SimpleITK в Numpy и обратно:

np_array = sitk.GetArrayFromImage(image)

new_image = sitk.GetImageFromArray(np_array)
new_image.SetOrigin(image.GetOrigin())
new_image.SetDirection(image.GetDirection())
new_image.SetSpacing(image.GetSpacing())

Стоит отметить, что порядок размерностей изображений в форматах SimpleITK и Numpy отличается: в SimpleITK – [x, y, z], в Numpy [z, y, x].

Resampling

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

new_image = sitk.Resample(
  image,                       # Изображение (объект класса Image) с системой координат О1
  referenceImage,              # Cетка (объект класса Image) с системой координат О2
  transform,                   # Пространственное преобразование, переводящее О1 в О2
  interpolator,                # Интерполятор, определяющий значение пикселей в системе О2
  defaultPixelValue = 0.0,      # Значение пикселей, которые не могут быть рассчитаны
  outputPixelType = sitkUnknown,# Тип данных нового изображения
)

Кроме нее можно пользоваться классом ResampleImageFilter объектного API.

SimpleITK включает три варианта для указания пространственной сетки:

Сетка, которая определена исходным изображением (т.е. фактически переноса нет, выполняется преобразование «на месте»);
Сетка, которая определена другим изображением;
Сетка, заданная через параметры: размер, начало координат, интервалы и матрица направляющих косинусов.

SimpleITK содержит множество методов интерполяции. Но наиболее часто используются два:

sitkLinear – используется для большинства задач и является компромиссом между точностью и вычислительной эффективностью;
sitkNearestNeighbor – используется для интерполяции сегментационных масок. Не вносит новые метки в результат, но может приводить к искажениям на границе масок. Вместо него можно использовать sitkLabelLinear или sitkLabelGaussian – оба более медленные (особенно второй), но более точные.

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

Transform

SimpleITK поддерживает два типа пространственных преобразований: с глобальной (неограниченной) пространственной областью и с ограниченной. В SimpleITK доступны различные глобальные 2D и 3D преобразования (сдвиг, вращение, растяжение и т.д.). Ограниченные преобразования определяются как тождественные за пределами своей области. К ним относятся преобразование B-сплайна, часто называемое деформацией свободной формы, и преобразование поля смещения. Также SimpleITK поддерживает составное преобразование либо с ограниченной, либо с глобальной областью. Оно представляет собой несколько преобразований, применяемых одно за другим в порядке стека (LIFO).

Для создания глобального преобразования достаточно выбрать нужный класс и указать размерность – 2 или 3. Преобразования имеют настраиваемые параметры, каждое свои, например:

affine = sitk.AffineTransform(2)
affine.SetTranslation((x_translation, y_translation))

Применение преобразования к изображению осуществляется через ресэмплинг, например:

new_image = sitk.Resample(image, referenceImage, affine, sitk.sitkLinear)

Преобразование без переноса можно выполнить с помощью альтернативной сигнатуры без указания параметра referenceImage:

new_image = sitk.Resample(image, affine, sitk.sitkLinear)

Если же нужен перенос без преобразования, то можно использовать базовый класс Transform как тождественное преобразование:

new_image = sitk.Resample(image, referenceImage, sitk.Transform(), sitk.sitkLinear)

Registration

Задача регистрации серий в SimpleITK решается с помощью нелинейной оптимизации, которая подбирает параметры преобразования с минимизацией метрики, определяющей схожесть исходного (moving) и целевого (fixed) изображений. Поддерживаются только изображения с типом данных sitkFloat32 или sitkFloat64.

Использование метрики сходства часто приводит к неудачной регистрации сильно отличающихся серий даже при близости пространственных характеристик, например, T2W и DWI с высоким b-фактором.

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

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

Регистрация изображений в SimpleITK выполняется функцией ImageRegistrationMethod, которая создает нужное преобразование, передаваемое в ресэмплинг. Подробнее о настройках преобразования регистрации можно посмотреть в документации.

Вообще говоря, повторные регистрации одних и тех же наборов данных дадут немного разные результаты. Это связано с реализацией вычисления метрики сходства, которая использует как рандомизацию, так и многопоточность. Можно установить нужные параметры для воспроизводимости результатов, но они сильно замедляют вычисления.

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

6.2. PyDicom

PyDicom – это библиотека для Python с открытым исходным кодом, предназначенная для чтения, изменения и записи файлов в формате DICOM.

Эта библиотека удобнее, чем SimpleITK в части доступа к тэгам – можно получить значения всех, которые есть в файле, хотя само чтение выполняется медленнее. Также можно сформировать и сохранить итоговый файл в разных вариантах формата DICOM – библиотека наряду с простыми и удобными функциями предоставляет возможность полного контроля над его структурой и содержанием. Это может оказаться полезным, например, при создании отчетов DICOM SR в структуре, требуемой заказчиком.

При необходимости можно загрузить только тэги без собственно данных изображения – это ускоряет первичную обработку исследований, например, при выборе файлов, принадлежащих только нужной серии.

tags = pydicom.dcmread(file_name, stop_before_pixels = True)

ds = pydicom.dcmread(file_name)
image = ds.pixel_array

При этом данные будут загружены в виде Numpy-массива размерностью [z, y, x] для 3D и [y, x] для 2D.

Можно просмотреть все тэги файла в удобном структурированном виде, просто напечатав загруженный файл:

Рисунок 44. Результат вывода на экран содержимого DICOM-файла [источник]

Доступ к каждому тэгу осуществляется либо по его коду, либо по названию (не для всех тэгов):

elem = ds[0x0008, 0x0016]
>>> elem
(0008, 0016) SOP Class UID                       UI: CT Image Storage

>>> elem.keyword
'SOPClassUID'

elem = ds['SOPClassUID']
>>> elem
(0008, 0016) SOP Class UID                       UI: CT Image Storage

>>> elem.value
'1.2.840.10008.5.1.4.1.1.2'

>>> ds.SOPClassUID
'1.2.840.10008.5.1.4.1.1.2'

Иногда бывает полезно преобразовать загруженный DICOM-файл в иерархический словарь для его обработки средствами самого Python:

ds.to_json_dict() – для получения основных тэгов
ds.file_meta.to_json_dict() – для получения метаинформации файла

Для создания нового DICOM-файла в полностью контролируемом режиме достаточно сформировать аналогичные словари тэгов в структуре, определенной стандартом и\или заказчиком, и записать их в нужные поля. Пример для бинарных данных с упаковкой их в побитную структуру:

dcm = pydicom.dataset.Dataset()
dcm.update(dcm.from_json(ds_json_dict))

meta = pydicom.Dataset()
meta.update(meta.from_json(file_meta_json_dict))
dcm.file_meta = meta

dcm.PixelData = pydicom.pixel_data_handlers.util.pack_bits(pixel_data)

pydicom.dcmwrite(file_name, dcm, write_like_original = False)

6.3. PyNrrd

Формат NRRD (Nearly Raw Raster Data) предназначен для научной и медицинской визуализации и обработки изображений с использованием N-мерных данных. Это простой и гибкий формат файла с заголовком, содержащим информацию о данных, которую можно прочитать, просто открыв файл NRRD в текстовом редакторе.

Этот формат поддерживает множество библиотек, включая SimpleITK, но иногда удобнее работать с ним с помощью нативной библиотеки PyNrrd:

Данный формат может оказаться полезным для быстрого и простого сохранения промежуточных данных. Например, при обучении нейронных сетей на 2D-изображениях (особенно, собранных из разных серий) формирование батчей «на лету» очень ресурсоемко, т.к. требует загрузки в память всех серий полностью для выбора из каждой единственного слайса. Альтернативой кэшированию может стать предварительное формирование датасета из собранных слайсов, сохраненных в формате NRRD со всей необходимой для обучения метаинформацией (например, размерами пикселей). Аналогично можно сохранять и сегментационные маски.

При этом необязательно соблюдать все требования формата – благо, поля заголовка являются необязательными. Достаточно записать только нужную информацию и сохранить вместе со слайсом. Можно даже и вовсе не использовать заголовок.

Хотя возможно чтение из NRRD-файла отдельно данных и отдельно заголовка, удобнее вызывать одну функцию, которая сама последовательно выполняет эти два шага:

image, header = nrrd.read(file_path)

При этом данные будут загружены в виде Numpy-массива размерностью [z, y, x] для 3D и [y, x] для 2D. Заголовок представляет собой OrderedDict с перечнем доступных полей.

Если требуется сохранение только данных, то можно записывать файл без заголовка. Минимум необходимых полей будет сгенерирован автоматически, и при чтении заголовок будет все равно возвращен.

nrrd.write(file_path, image)

Для сохранения в заголовке своих данных можно просто сформировать словарь и передать его как заголовок. Поля могут соответствовать требованию формата NRRD по структуре, а могут быть полностью кастомными. При сохранении заголовок будет дополнен таким же минимумом необходимых полей.

header = {'space directions': [[0.39, 0], [0, 0.39]], 'custom field': 123456}
nrrd.write(file_path, image, header)

При использовании NRRD-формата для создания промежуточного датасета рекомендуется устанавливать поле заголовка encoding в значение 'raw', т.к. по умолчанию используется gzip-сжатие, что сильно замедляет чтение таких файлов.

6.4. MONAI

MONAI – это набор активно развивающихся фреймворков PyTorch с открытым исходным кодом для глубокого обучения в медицинской визуализации, часть экосистемы PyTorch.

Сайт и документация: https://monai.io/
Исходники: https://github.com/Project-MONAI

MONAI создан NVIDIA и Королевским колледжем Лондона как инструмент для обмена передовым опытом в области ИИ в сфере медицинской визуализации. Фреймворки проекта охватывают весь жизненный цикл ИИ-приложений от аннотирования данных до развёртывания и оптимизации:

MONAI Core: специализированная платформа для обучения моделей ИИ в области медицинской визуализации.
MONAI Label: интеллектуальный инструмент для разметки изображений, позволяющий быстро аннотировать новые наборы данных.
MONAI Deploy App SDK: инструмент создания и внедрения ИИ-приложений.
MONAI Model Zoo: коллекция предобученных моделей медицинской визуализации в формате MONAI Bundle.

MONAI Architecture — Рисунок 45. Архитектура проекта MONAI [источник]

MONAI Core содержит большое количество различного инструментария – лосс-функции, метрики, оптимизаторы, аугментации, заготовки в виде отдельных слоев и блоков, а также уже готовые реализации разных архитектур.

Среди этих инструментов можно найти не только стандартные, но и специфические для медицинских данных. Основной акцент при этом делается на обучение на 3D-данных, что оправданно по сути, однако проблематично для небольших датасетов.

Подробное описание данного набора фреймворков здесь не предусмотрено, т.к. это был бы неоправданно большой объем – с необходимыми каждому темами можно познакомиться в документации. Но в дальнейшем эксперименты по возможности будут включать элементы MONAI.

МРТ для DataScience. Часть 6