Продвинутый подход к обнаружению границ на примере стенок сосуда

Интересная информация

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

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

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

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

Построение центральной линии

Самый трудный в плане реализации этап (по крайней мере, на него ушло больше всего времени). Метод основан на использование матрицы Гессе (Multi-Scale Vessel Segmentation Using Hessian Matrix Enhancement). Подробнее в уже упомянутой статье.

Формирование срезов

У нас есть только центральная линия, и нужны пространственно-зависимые интенсивности вокселей, с которыми можно удобно работать. Чтобы их получить собирается “стопка” из срезов. Для начала через фиксированные расстояния на центральной линии ставятся точки. Затем из каждой точки строится перпендикуляр . После находится второй перпендикуляр . Где — направление центральной линии в точке. Оба перпендикуляра нормализованы. В каждой точке центральной линии вектора образуют 2d систему координат. Таким образом формируются срезы:

Положение вокселя определим как , где , здесь — реальные координаты вокселя, k — номер среза. Обратная формула для реальных координат: . При переходе к новой системе координат пространство, формируемое срезами, упрощается:

То что нам нужно!

Построение внешней границы сосуда

Давайте взглянем на схему:

Нашу стопку срезов, полученную на предыдущем этапе, мы режем восемью плоскостями (по аналогии с разрезанием торта) и все вычисления проводим уже в пространстве плоскостей:

Разрез

Если отобразить нормализованные значения интенсивностей вокселей, попавших на режущую плоскость, то получим такую картину:

Для обнаружения границ сосуда используется классический подход (edge detection by gradient) совместно с поиском пути. Схема:

1. Применяем сглаживание Гаусса с небольшим значением для подавления шума:
Для точки с координатами :
, где возвращает значение интенсивности в точке ; r обычно принимает значение ( — округление в большую сторону); — коэффициент сглаживания.

2. В каждой точке находим градиент и его значение, вычисления проводятся со сглаженными значениями интенсивности:
,
где — частные производные. Они находятся методом конечных разностей:

,
где — интенсивность в точке после сглаживания.
Далее нам потребуются направление градиента ( — нормализация вектора) и значение градиента ( — длина вектора)

3. Направление градиента переводим в градусы или радианы:
(atan2() — функция арктангенса на C++, не путать с atan() ), затем угол округляем так, что он может иметь 4 значения с шагом 45 градусов, т.е. верх и низ считаются одним и тем же направлением:

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

5. Все оставшиеся воксели считаются границами. На основании значения градиента, порога кальция (который доступен не сразу) и близости к “вертикальному” центру, каждой точке присваивается определенная стоимость (чем ярче воксель, тем выше его приоритет при поиске пути):

В таком виде границы сосуда уже практически однозначно определены.

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

Результат

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

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

Хороший результат

Обнаружение порогов внутренней, внешней границ и порога кальция

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

Обнаружение порогов

Теперь рассмотрим алгоритм кластеризации expectation maximization (далее EM). Начнем с функции нормального распределения: она характеризуется математическим ожиданием и среднеквадратическим отклонением . Вот как они влияют на вид распределения:

Предположим, у нас есть данные (точки), которые пришли из “желтого” источника и из “синего” источника:

Тогда, используя стандартные формулы, мы легко находим математическое ожидание и среднеквадратическое отклонение для каждого источника. Формулы для источника “a”:

Но что делать, если мы знаем количество источников, но не знаем, какие точки какому источнику принадлежат? Что делать, если у нас такая картина?

Если бы кто-то пришел и сказал нам параметры распределений, то мы могли бы легко вычислить вероятность принадлежности каждой точки к каждому из источников. Вероятность того что точка принадлежит источнику “a”:

, где




А если нам надо вычислить параметры источника, зная вероятности:

Таким образом получается замкнутый круг: если бы мы знали параметры источников — мы бы вычислили, какая точка какому источнику принадлежит, но мы не знаем параметры. А если бы мы знали, какая точка какому источнику принадлежит — мы бы вычислили их параметры, но мы не знаем, какая точка какому источнику принадлежит. Балансирование между эти фактами именно то, чем занимается алгоритм EM.

На старте EM получает некоторые предопределенные параметры для источников, которые могут быть просто выбраны случайно. Очевидно, что если известны параметры, то мы можем вычислить вероятность принадлежности каждой точки каждому из источников. Теперь, когда вероятности известны, мы можем вычислить новые более точные параметры. Затем всё начинается сначала, но уже с новыми параметрами. После каждого цикла параметры источников становятся всё точнее.

Как мы можем использовать эти знания по отношению к сосудам? Давайте взглянем на структуру одного из них:

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

— жир
— стенка #1
— стенка #2
— контрастное вещество
— кальций

Распределение интенсивностей вокселей в каждом случае является нормальным. Т.е. у нас есть всё необходимое, чтобы, используя EM, найти параметры каждого источника.

Результаты получаются достаточно хорошие

Зеленая линия — гистограмма интенсивностей, красная — полученная математическая модель.

Теперь, когда мы знаем параметры каждого источника, мы можем вычислить пороги — значения интенсивностей, при пересечении которых, принадлежность вокселя меняется с одного источника на другой. Нас интересуют:

1. Порог внешней границы сосуда. Если интенсивность вокселя ниже этого значения, то считается, что он вообще не принадлежит сосуду;

2. Порог внутренней границы сосуда. Если интенсивность вокселя больше этого значения, то он
относится к просвету сосуда, т.е. к смеси крови и контрастного вещества;

3. Порог кальция. Если значение интенсивности вокселя больше этого значения, то он относится к кальцию.

Построение внутренней границы сосуда

Как всегда начнем со схемы, вычисления выполняются для каждого среза.

Если визуально отобразить данные согласно порогам, полученным на предыдущем этапе, мы получим такую картину:

Красный цвет — стенка сосуда. Зеленый цвет — просвет. Белый цвет — кальций.

Первое, что бросается в глаза — это “висячий” кальций, который не примыкает ни к одной из стенок. Это считается нормальным и возникает вследствии сглаживания, применяемого самим томографом.

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

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

1. Три квадрата одного знака и один противоположного, движение контура происходит по диагонали:

Пример

2. Два квадрата одного знака и два противоположного, причем квадраты с одинаковым знаком находятся по одну сторону, движение контура идет вертикально или горизонтально:

Пример

3. Два квадрата одного знака и два противоположного, квадраты с одинаковыми знаками размещены по разные стороны:

Пример

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

Пример

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

Пример

Конкретно первый и второй случаи:

Для каждого среза сосуда мы находим два основных контура — это контур внешней границы и контур внутренней границы. Внешний контур мы сразу “подрезаем” нашим другим внешним контуром, который мы нашли вначале статьи поиском путей. Это сделано, чтобы проигнорировать ответвления сосуда. Контуры кальция, которые находятся слишком далеко от внутренней стенки мы игнорируем, как-будто их вообще не существует, остальные находим и используем в дальнейшем. Если центр сосуда находится внутри кальция, то мы смещаем его в ближайшем к контуру кальция направлении, пока центр не окажется в просвете (в зеленой области). Такой обновленный центр я буду называть стартовой позицией.

Согласно схеме может быть два случая: простой и сложный.

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

Чтобы добиться нужного результата был разработан специальный алгоритм, в основу которого легли идеи, используемые в физических 2d движках, в частности polygons collisions solving и separate axis theorem.

Два базовых понятия, без которых не обойтись: для пересекающихся выпуклых многоугольников mtv-вектор (minimum translation vector) — это наименьшее смещение одного из многоугольников, после которого пересечение прекращается.

Ещё нам потребуются нормали многоугольника — в 2D они перпендикулярным граням и указывают “наружу”:

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

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

2. Если вершина стороннего контура попала внутрь стартового контура, то на ближайшую грань стартового контура применяется смещение пропорциональное mtv-вектору вершины;

3. Если вершина стартового контура попала внутрь стороннего контура, то на вершину стартового контура применяется смещение пропорциональное mtv-вектору вершины. Это вместе с предыдущим пунктом не позволяет контуру выйти за пределы границ других контуров;

4. Если для вершины не сработали случаи 2 и 3, то на нее применяется сила в усредненном направлении двух нормалей соседних граней. Это обеспечивает рост контура “вширь”.

Важно: нельзя смещать вершину или грань полностью на длину mtv-вектора — его умножают на некий коэффициент в пределах от 0.2 до 0.8. Коэффициенты для каждой силы в остальных случаях подбираются экспериментально.

Благодаря такому подходу мы находим просвет сосуда с учетом того, что кальций не примыкает вплотную к стенкам. Теперь просто объединяем результаты со всех срезов и получаем внутреннюю и внешнюю границы сосуда:

Кажущаяся неточность границы после стента вызвана аномальным раздвоением:

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