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

7. Задача выбора нужной серии из исходного набора данных

7.1. Постановка задачи

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

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

Основная проблема выделения серий состоит в том, что в тэгах DICOM-файлов далеко не всегда указаны параметры, по которым можно это сделать просто и быстро. Тэг «Series Description», который часто хранит название и параметры серии, является опциональным и может отсутствовать. Так же он может содержать неполную или некорректную информацию, т.к. допускается изменение его значения оператором.

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

Задачей автоматического выбора серий занимаются давно и применяют для этого разные подходы (например) – от анализа тэгов до обучения классификационных сетей. Однако, полного решения данной задачи со 100% точностью пока не найдено.

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

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

При решении каждой задачи необходимо уточнить (в том числе у заказчика) технические требования к МРТ-данным, на основе которых уже алгоритмом будет приниматься решение о корректности или дефектности того или иного исследования по наличию в нем нужных серий с нужными характеристиками.

Охватить все многообразие сложно, поэтому далее на примере 2D аксиальной T2W серии будет представлен общий подход к выбору серий.

В качестве метрик, оценивающих качество отбора, возьмем следующие:

Выбранные серии относятся к нужному типу – максимальная точность. Полнота не оценивается, т.к. количество серий искомого типа, но не подходящих для отбора, может быть достаточно велико.
Выбранные серии представлены во всех исследованиях. Отсутствие таких серий в данных может считаться нарушением протокола исследования. Это – некий аналог полноты, который здесь будет уравновешивать предыдущий критерий точности. Назовем его «полнота по исследованиям».

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

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

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

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

7.2. Выбор серий с изображениями. MRI IOD и Enhanced MRI IOD

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

Для этого вместо фильтрации по модальности (тэг (0008,0060) Modality) воспользуемся более надежным вариантом - DICOM-тэгом (0008,0016) SOPClassUID. Он имеет фиксированный набор значений для разных модальностей. Для каждого в стандарте DICOM описаны набор и структура используемых тэгов. Здесь мы ограничимся двумя значениями:

1.2.840.10008.5.1.4.1.1.4 (MR Image IOD) — старый и широко используемый формат;
1.2.840.10008.5.1.4.1.1.4.1 (Enhanced MR Image IOD) — более современный и удобный, но реже встречающийся формат.

В Enhanced MR Image IOD задача выделения некоторых типов серий решается тривиально с помощью одного единственного обязательного тэга (0008, 9209) Acquisition Contrast, который может принимать, в том числе, следующие значения: DIFFUSION, PROTON_DENSITY, STIR, T1, T2, UNKNOWN.

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

(0018, 9005) Pulse Sequence Name — коммерческое название ИП.
(0018, 9008) Echo Pulse Sequence — класс ИП. Содержит одно из значений SPIN, GRADIENT, BOTH.
(0018, 9018) Echo Planar Pulse Sequence — используется ли метод EPI. Содержит одно из значений YES, NO.
(0018, 9009) Inversion Recovery - используется ли один из методов IR, расположен на 3-м уровне вложенности внутри тэгов (5200, 9229) Shared Functional Groups Sequence и (0018, 9115) MR Modifier Sequence.
(0018, 9032) Geometry of k-Space Traversal — геометрия аксиального К-пространства. Содержит одно из значений RECTILINEAR, RADIAL, SPIRAL.
(0018, 0080) Repetition Time – TR, расположен на 3-м уровне вложенности внутри тэгов (5200, 9229) Shared Functional Groups Sequence и (0018, 9112) MR Timing and Related Parameters Sequence.
(0018, 0091) Echo Train Length - ETL, расположен на 3-м уровне вложенности внутри тэгов (5200, 9229) Shared Functional Groups Sequence и (0018, 9112) MR Timing and Related Parameters Sequence.
(0018, 9082) Effective Echo Time – TE, хранит значения для каждого слайса, расположен на 3-м уровне вложенности внутри тэгов (5200, 9230) Per-frame Functional Groups Sequence и (0018, 9114) MR Echo Sequence.

Перечень не полный и для тонкой настройки могут пригодиться и другие тэги — их можно посмотреть в описании DICOM-стандарта для этого формата.

Далее же сосредоточимся на решении поставленной задачи для данных в старом формате.

7.3. Определение параметров импульсных последовательностей в MR Image IOD

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

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

Ключевым атрибутом для определения класса ИП является обязательный для заполнения DICOM-тэг (0018,0020) ScanningSequence, который содержит набор из одного или нескольких значений:

SE – Spin Echo
GR – Gradient Recalled
IR – Inversion Recovery
EP – Echo Planar
RM – Research Mode

Значения SE и GR как раз и определяют класс ИП и в рассматриваемом тэге не могут встречаться вместе. И хотя существуют гибридные последовательности, смешивающие спиновое и градиентное эхо, в рамках рассматриваемой задачи и с учетом ограничений старого стандарта DICOM можно считать, что классы ИП спин-эха и градиентного эха не пересекаются. А если и да, то пересечение нас не интересует.

Значения IR и EP указывают на использование соответствующих методов. RM часто применяется для обозначения других методов, для которых отдельной аббревиатуры в стандарте DICOM не предусмотрено, например, для PROPELLER, причем, часто вместо SE и GR, что усложняет задачу определения класса ИП.

Поэтому для выделения, например, класса ИП спин-эхо, можно использовать два разных подхода – с максимальной точностью (только с наличием значения SE) или с максимальной полнотой (только без значения GR) и последующим уточнением. Во втором случае надо иметь в виду, что возможно отсутствие в данных всех трех значений – SE, GR и RM.

Два «соседних» тэга, из которых можно получить дополнительную информацию об ИП:

(0018,0021) SequenceVariant – обязательный, содержит одно или более значений.
(0018,0022) ScanOptions – обязательный, но может быть пустым или содержать другие значения, кроме перечисленных в стандарте DICOM.

Эти тэги далеко не всегда заполняются полностью (т.е. если значение отсутствует, не значит, что его нет), носят технический характер, и без глубокого понимания физического смысла МРТ и подходов к выбору настроек сканирования их сложно интерпретировать. Особенно это касается значений, добавляемых производителями томографов. Поэтому их стоит использовать с осторожностью, но они могут оказаться полезными, например, при поиске причин стабильных ошибок нейронных сетей.

Времена, описанные в разделе 3.1., представлены почти одноименными тэгами:

(0018,0080) RepetitionTime – TR, условно обязательный, хотя на практике заполняется всегда;
(0018,0081) EchoTime – TE, обязательный, но может быть пустым, хотя на практике заполняется всегда;
(0018,0082) InversionTime – TI, условно обязательный.

Следующий тэг, полезный для уточнения характера ИП – (0018,0091) EchoTrainLength, который упоминался ранее. Он обязательный, но может быть пустым. Если не пуст, то содержит одно целое число.

В сочетании с тэгами времен, данное значение помогает определять характер импульсной последовательности. Например, короткие TR и TE в сочетании с ETL равным единице соответствуют классической ИП T1 Spin Echo, которая часто применяется для получения T1W и DCE серий – она достаточно быстрая и обеспечивает хорошее качество. Для T2W серий нужны длинные TR и TE, поэтому для сокращения времени сканирования используются более быстрые последовательности, для которых ETL обычно гораздо больше единицы.

7.4. Определение геометрии данных

Начнем с тэга (0018,0023) MR AcquisitionType, который характеризует размерность К-пространства и содержит значение 2D или 3D (для Enhanced MR Image IOD возможно значение 1D, но это не наш случай). Хотя 3D-данные более удобны и информационно насыщенны для DS-задач, они не всегда включаются в протоколы и присутствуют далеко не во всех исследованиях. Поэтому есть смысл или выбирать только одно значение, или балансировать датасет по этим размерностям вместе с другими характеристиками данных.

Основные тэги с геометрическими характеристиками, часть которых уже упоминалась ранее:

(0028,0030) Pixel Spacing – вектор размеров пикселей в мм
(0018,0050) Slice Thickness – толщина слайса в мм
(0018,0088) Spacing Between Slices – расстояние между центрами слайсов в мм
(0028,0010) Rows и (0028,0011) Columns – размеры слайсов в пикселях

Тип проекции определяется по матрице направляющих косинусов в тэге (0020,0037) Image Orientation (Patient). Эта матрица хранится как одномерный вектор, первые три значения соответствуют косинусам первой строки слайса, вторые – первого столбца: [Xx, Xy, Xz, Yx, Yy, Yz]. Например, [1, 0, 0, 0, 1, 0] соответствует прямой аксиальной проекции.

Однако, там могут быть и отрицательные значения, если направление осей локальной и глобальной СК противоположные. А могут быть и меньше 1, если проекция «косая». Поэтому тип проекции лучше определять по модулю и с допуском, например | Xx | >= 0.8. Пороговое значение нужно определять по имеющимся данным. Для некоторых протоколов исследования «косые» проекции допускаются или даже рекомендуются.

7.5. Алгоритм выбора 2D аксиальной T2W серии

В качестве примера возьмем 2D аксиальную T2W серию. Для нее опишем все шаги отбора, начиная с выбора IOD.

Оставим только серии, для которых SOPClassUID равен 1.2.840.10008.5.1.4.1.1.4 – здесь рассмотрим только один IOD, старый, но можно прописать и вторую ветку выбора для Enhanced MR Image IOD.
Оставим только серии, для которых MR AcquisitionType равен «2D», поскольку в рассматриваемом случае такие данные встречаются в каждом исследовании.
Оставим только аксиальные проекции – нулевое и четвертое значения тэга Image Orientation (Patient) по модулю больше или равны 0.8.
Отбросим все серии, для которых тэг ScanningSequence содержит значения «GR» и «IR», т.к. подавляющее большинство T2W серий формируется и использованием ИП спин-эхо, а не градиентного эха. Если такие серии и встретятся – считаем их выбросами. И подавление сигналов от отдельных тканей так же не учитываем. Для дополнительной подстраховки на случай некорректных значений этого тэга можно добавить условие: InversionTime отсутствует или меньше или равен нулю (может быть установлен в -1).
Оставим только серии с длинным TR, т.е. больше 1500 – используем значение тэга RepetitionTime. Данный шаг требует контроля и, возможно, изменения порога, т.к. во многих исследованиях предпочитают выбирать быстрые ИП или особые методы для сокращения времени сканирования. На используемых для проверки алгоритма данных этот порог подошел.
Оставим только серии с длинным TE, т.е. больше 90 – используем значение тэга EchoTime. Данный шаг требует контроля и, возможно, изменения порога, т.к. во многих исследованиях предпочитают выбирать быстрые ИП или особые методы для сокращения времени сканирования. На используемых для проверки алгоритма данных этот порог был снижен до 65.
Оставим только серии с ETL больше 1 – тэг EchoTrainLength обязательный, но может быть пустым. Как уже говорилось ранее, для T2W серий с длинными TR и TE чаще всего используются быстрые ИП с достаточно высоким значением ETL.

На используемых для проверки алгоритма данных точность выделения T2W серий после этих шагов составила 98%. Причем, ошибочно были добавлены DWI серии, которые по рассмотренным параметрам логично могут считаться подмножеством T2. Полнота по исследованиям составила 99%, но ошибки представляли собой явное нарушение протокола исследования.

Стоит упомянуть, что ADC-серии часто копируют параметры DWI-серий, по которым рассчитаны, поэтому тоже могут попадать в формируемую выборку. При проверке этого не случилось, но такая ситуация вполне возможна.

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

Вернемся к точности и попробуем убрать DWI серии из выборки и подстраховаться насчет попадания в нее ADC.

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

Однако для данного типа серий характерен особый хаос в данных. В частности, стандарт DICOM для MR Image IOD вообще не предусматривает подобных тэгов, поэтому производители оборудования записывают эти данные куда попало. В более новом формате Enhanced MR Image IOD для описания диффузионных данных предусмотрен целый раздел:

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

Несколько примеров из практики по значениям тэга (0018, 9087) Diffusion b-value:

Тэг даже в старом формате MR Image IOD корректно заполнен для всех файлов DWI-серии и отсутствует в файлах других серий.
Тэг присутствует только для части файлов DWI-серии – для которых b-value больше нуля.
Тэг присутствует во всех файлах исследования, но для недиффузионных серий его значение равно нулю.
Тэг отсутствует во всех файлах всех серий, но информация о диффузии, включая b-фактор, сохранена в приватных тэгах, определенных производителем томографа.
Тэг отсутствует во всех файлах всех серий, информация о диффузии отсутствует в явном виде, но частично может быть восстановлена из некоторых необязательных описательных тэгов, например, (0008,103E) Series Description, (0020,4000) Image Comments, (0018,1030) Protocol Name, (0008,0008) Image Type и др.

Наиболее надежным для используемых данных оказался последний вариант – проверка по всем перечисленным тэгам: если тэг присутствует в данных, то его значение не должно включать следующие строки в любом написании: «DWI», «DIFFUSION», «FOCUS», «ADC», «Apparent» (последние два значения добавлены для исключения серий ADC).

После применения этого фильтра точность составила 100%.

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

Стоит также обозначить, что проверка выбора 2D аксиальной T2W серии осуществлялась для данных, подготовленных для решения единственной задачи по соответствующему протоколу исследования. Подобный алгоритм для других задач и данных, несомненно, будет нуждаться в адаптации как по пороговым значениям, так и, возможно, по используемым тэгам или другим характеристикам данных.

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