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

в 4:06, , рубрики: clusterization, DICOM, DICOM Viewer, edge detection, inobitec, mri, Алгоритмы, ангиография, Блог компании Inobitec, инобитек, кластеризация, МРТ, обнаружение границ, обработка изображений, сосуды, Софт

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Разрез

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

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

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

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

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

1. Применяем сглаживание Гаусса с небольшим значением Продвинутый подход к обнаружению границ на примере стенок сосуда - 18 для подавления шума:
Для точки с координатами Продвинутый подход к обнаружению границ на примере стенок сосуда - 19:
Продвинутый подход к обнаружению границ на примере стенок сосуда - 20, где Продвинутый подход к обнаружению границ на примере стенок сосуда - 21 возвращает значение интенсивности в точке Продвинутый подход к обнаружению границ на примере стенок сосуда - 22; r обычно принимает значение Продвинутый подход к обнаружению границ на примере стенок сосуда - 23 (Продвинутый подход к обнаружению границ на примере стенок сосуда - 24 — округление в большую сторону); Продвинутый подход к обнаружению границ на примере стенок сосуда - 25 — коэффициент сглаживания.

2. В каждой точке находим градиент и его значение, вычисления проводятся со сглаженными значениями интенсивности:
Продвинутый подход к обнаружению границ на примере стенок сосуда - 26,
где Продвинутый подход к обнаружению границ на примере стенок сосуда - 27 — частные производные. Они находятся методом конечных разностей:

Продвинутый подход к обнаружению границ на примере стенок сосуда - 28,
где Продвинутый подход к обнаружению границ на примере стенок сосуда - 29 — интенсивность в точке Продвинутый подход к обнаружению границ на примере стенок сосуда - 30 после сглаживания.
Далее нам потребуются направление градиента Продвинутый подход к обнаружению границ на примере стенок сосуда - 31 (Продвинутый подход к обнаружению границ на примере стенок сосуда - 32 — нормализация вектора) и значение градиента Продвинутый подход к обнаружению границ на примере стенок сосуда - 33 (Продвинутый подход к обнаружению границ на примере стенок сосуда - 34 — длина вектора)

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

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

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

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

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

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

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

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

Результат

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Пример

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

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

Пример

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

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

Пример

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

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

Пример

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

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

Пример

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Автор: Антон Шаталов

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js