[recovery mode] Решение проблемы обнаружения центральной линии сосуда18.10.2018 16:04

Суть задачи

В процессе медицинской диагностики может возникнуть необходимость исследовать сосуды пациента. Такое исследование называется ангиографией. С появлением томографов в дополнение к классической ангиографии появились методы МРТ и КТ ангиографии, которые в отличие от традиционной ангиографии, дающей только плоскую картинку в одной проекции, позволяют получить полное трехмерное представление сосудов. Для проведения таких исследований пациенту в кровь вводится контраст — специальное вещество, делающее сосуды на снимках более яркими. В зависимости от предполагаемого диагноза, врач или оценивает общую картину, или пытается найти конкретные участки сосудов, в которых возникли проблемы. Если участок сосуда сужен и пропускает меньше крови, чем должен, то это место называется стенозом.

Одна из задач врача — найти стенозы и оценить, насколько они опасны. Задача же разработчика, как обычно, облегчить работу конечного пользователя. Для этого необходимо построить полную 3D модель стенок сосуда и провести их первичный анализ. Это является большой и интересной задачей, однако, в её основе лежит более простая и известная проблема — построение центральной линии сосуда.

Первая попытка

Перед продолжением чтения этой статьи желательно чуть-чуть ознакомиться с представлением данных, получаемых в результате работы томографов. Это можно сделать, прочитав нашу давно написанную статью про воксельный рендер DICOM Viewer изнутри. Если коротко: есть 3D-массив чисел, в каждом элементе которого хранится значение сигнала (интенсивность). Этот массив называется объемом. Сам элемент является вокселем, а его индексы в 3D-массиве будут являться 3D-координатами. Перед рендером каждого вокселя, его интенсивность обрабатывается функцией, и воксель приобретает определенные цвет, яркость и прозрачность.

Относительно задачи, первое с чем столкнулся именно я — это то, что наш рендер позволяет уже на визуальном уровне показать все сосуды. А именно, из непонятной «каши», как на рисунке слева, поигравшись с настройкам можно сделать вполне очевидные сосуды, как на рисунке справа:

Кажется, что задача решена: сосуды «вот они», на ум сразу приходит region growing алгоритм: мы знаем цвет и прозрачность нужных нам вокселей, значит можно итерационно перебирать их соседей, пока не будет проложен путь между указанными точками.

Сегментируем сосуд

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

Пример

Сосуд может загибаться или прилегать вплотную к самому себе:

Пример

Участки сосуда могут быть очень тонкими и даже прерываться:

Пример

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

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

Классический подход

Основной подход к проблеме обнаружения сосудов (vessel detection), заключается в вычислении собственных значений матрицы Гессе (Hessian matrix). Также в последнее время к этой задаче пробуют подключить нейронные сети. Однако, стандартное решение выглядит более надежным, потому что имеет упоминание в большем количестве литературы и используется не только для нахождения сосудов, но и для многих других задач (например, обнаружение объектов на астрономических снимках), поэтому от нейронных сетей пришлось отказаться.

В нашем случае матрица Гессе — это матрица, элементами которой являются частные производные второго порядка интенсивности в конкретном вокселе:

$H_{M}=\begin{pmatrix} \frac{\partial^2 f}{\partial^2 x} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\ \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial^2 y} & \frac{\partial^2 f}{\partial y \partial z} \\ \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial^2 z} \end{pmatrix}$

$\frac{\partial^2 f}{\partial^2 x},\frac{\partial^2 f}{\partial^2 y},\frac{\partial^2 f}{\partial^2 z},\frac{\partial^2 f}{\partial x \partial y},\frac{\partial^2 f}{\partial y \partial z},\frac{\partial^2 f}{\partial x \partial z}$ — соответственно, вторые частные производные.

После построения матрицы Гессе необходимо найти её собственные значения и собственные векторы решив следующее уравнение:

$\begin{vmatrix} \frac{\partial^2 f}{\partial^2 x}-\lambda_{1} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\ \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial^2 y}-\lambda_{2} & \frac{\partial^2 f}{\partial y \partial z} \\ \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial^2 z}-\lambda_{3} \end{vmatrix}=0$

$\lambda_{1},\lambda_{2},\lambda_{3}$ — собственные значения. Задача их нахождения сводится к решению кубического уравнения, но в нашем случае матрица симметрична, поэтому решение упрощается и может быть представлено небольшим псевдокодом:

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0) 
   eig1 = A(1,1)
   eig2 = A(2,2)
   eig3 = A(3,3)
else
   q = trace(A)/3
   p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
   p = sqrt(p2 / 6)
   B = (1 / p) * (A - q * I)
   r = det(B) / 2

   if (r <= -1) 
      phi = pi / 3
   else if (r >= 1)
      phi = 0
   else
      phi = acos(r) / 3
   end

   eig1 = q + 2 * p * cos(phi)
   eig3 = q + 2 * p * cos(phi + (2*pi/3))
   eig2 = 3 * q - eig1 - eig3
end

Где A — исходная матрица 3×3; I — единичная матрица 3×3; eig1, eig2, eig2 — собственные значения; trace () — возвращает след матрицы; det () — возвращает детерминант матрицы.

Теперь мы можем найти собственные векторы. Для этого в уравнение:

$\begin{vmatrix} \frac{\partial^2 f}{\partial^2 x}-\lambda & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\ \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial^2 y}-\lambda & \frac{\partial^2 f}{\partial y \partial z} \\ \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial^2 z}-\lambda \end{vmatrix}\cdot\begin{vmatrix} V_{x}\\V_{y}\\V_{z}\end{vmatrix}=0$

подставим одно из найденных собственных значений вместо $\lambda$ . Найденный вектор $\vec{V}=\{V_{x},V_{y},V_{z}\}$ и будет являться собственным вектором. Таких векторов получится три: $\vec{V_{1}},\vec{V_{2}},\vec{V_{3}}$ — по одному для каждого собственного значения. Собственные значения и векторы необходимо отсортировать в таком порядке ${\lvert}\lambda_{1}{\rvert}\leq{\lvert}\lambda_{2}{\rvert}\leq{\lvert}\lambda_{3}{\rvert}$ , при этом векторы тоже меняются местами, т.е. меняя местами значения $\lambda_{1}$ и $\lambda_{2}$ мы также меняем местами $\vec{V_{1}}$ и $\vec{V_{2}}$ . Если после сортировки выполняется условие $\lambda_{2}<0,\lambda_{3}<0$ , то считается, что воксель, в котором построили матрицу, принадлежит сосуду. При этом собственный вектор $\vec{V_{1}}$ будет указывать направление вдоль сосуда независимо от того как близко к стенке находится воксель:

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

$V_{\sigma}=\left(1-\frac{{\lvert}{\lvert}\lambda_{2}{\rvert}-{\lvert}\lambda_{3}{\rvert}{\rvert}}{{\lvert}\lambda_{2}{\rvert}+{\lvert}\lambda_{3}{\rvert}}\right)\cdot\left(\lambda_{1}\frac{2}{3}-\lambda_{2}-\lambda_{3}\right)\cdot\left(1-e^{-\frac{{\lambda_{1}}^2+{\lambda_{2}}^2+{\lambda_{3}}^2}{2}}\right)$

$V_{\sigma}$ — соответственно, значение сосудистости

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

Определяем направление сосуда в текущем вокселе,
Делаем перпендикулярный направлению срез сосуда и перемещаемся по градиенту в воксель среза с максимальной сосудистостью,
Определяем направление сосуда в вокселе с максимальной сосудистостью,
Делаем небольшой шаг в направлении сосуда и попадаем в новый воксель,
Возвращаемся в шаг 1.

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

Анимация алгоритма

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

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

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

$k(x)=\frac{1}{6}\begin{cases} (12-9B-6C)\lvert{x}\rvert^3+(-18+12B+6C)\lvert{x}\rvert^2+(6-2B),\quad if\>\lvert{x}\rvert<1\\(-B-6C)\lvert{x}\rvert^3+(6B+30C)\lvert{x}\rvert^2+(-12B-48C)\lvert{x}\rvert+(8B+24C),\\if\> 1\leq\lvert{x}\rvert<2 \\0,\quad otherwise\end{cases}$

где и имеют предопределенные значения. Если C=0 , то мы имеем дело с B-сплайном, если B=0 , то с кардинальным сплайном. В нашем случае B=0 , $C=\frac{1}{2}$ (Catmull-Rom filter). Тогда получаем:

$k(x)=\frac{1}{6}\begin{cases} 9\lvert{x}\rvert^3-15\lvert{x}\rvert^2+6,\quad if\>\lvert{x}\rvert<1\\-3\lvert{x}\rvert^3+15\lvert{x}\rvert^2-24\lvert{x}\rvert+12,\quad if\> 1\leq\lvert{x}\rvert<2 \\0,\quad otherwise\end{cases}$

Рассмотрим случай для 1D. Если даны значения в точках $p_{-1},p_{-0},p_{1},p_{2}$ c координатами -1, 0, 1, 2 , и мы интерполируем f(x) , где $x\in[0,1)$ , то мы можем вычислить вес для каждой вершины:

$\\w_{-1}(x)=\frac{1}{2}(x^2(2-x)-x)\\w_{0}(x)=\frac{1}{2}(x^2(3x-5)+2)\\w_{1}(x)=\frac{1}{2}(x^2(4-3x)+x)\\w_{2}(x)=\frac{1}{2}x^2(x-1)$

Итоговое проинтерполированное значение:

$f(x)=p_{-1}\cdot w_{-1}(x)+p_{0}\cdot w_{0}(x)+p_{1}\cdot w_{1}(x)+p_{2}\cdot w_{2}(x)$

Теперь рассмотрим наш полноценный случай для 3D пространства. Как нетрудно догадаться, мы будем интерполировать внутри куба размерностью 4×4х4 вокселя. За основу берем веса, рассмотренные для 1D случая:

$f(x_{0},y_{0},z_{0})=\sum_{i=-1}^{2}\sum_{j=-1}^{2}\sum_{k=-1}^{2}w_{i}(x_{0}-\lfloor x_{0}\rfloor)\cdot w_{j}(y_{0}-\lfloor y_{0}\rfloor)\cdot w_{k}(z_{0}-\lfloor z_{0}\rfloor)\cdot\\\cdot I(\lfloor x_{0}\rfloor+i,\lfloor y_{0}\rfloor+j,\lfloor z_{0}\rfloor+k)$

где $(x_{0},y_{0},z_{0})$ — дробные координаты точки, I(x,y,z) — интенсивность в вокселе с координатами (x,y,z) , $\lfloor\cdot\rfloor$ — округление в меньшую сторону.

Для кого-то может показаться интересным факт, что сумма весов всех вокселей равна 1:

$\sum_{i=-1}^{2}\sum_{j=-1}^{2}\sum_{k=-1}^{2}w_{i}(x-\lfloor x\rfloor)\cdot w_{j}(y-\lfloor y\rfloor)\cdot w_{k}(z-\lfloor z\rfloor)=1$

Детали

Вернемся к вычислению производных для матрицы Гессе. У нас есть координаты и интенсивность, как численно считается вторая частная производная все знают: основным является метод конечных разностей. Для точки (x,y,z) :

Формула

$\\\frac{\partial^2 f(x,y,z)}{\partial^2 x}={f(x+1,y,z)-2f(x,y,z)+f(x-1,y,z)}\\\frac{\partial^2 f(x,y,z)}{\partial^2 y}={f(x,y+1,z)-2f(x,y,z)+f(x,y-1,z)}\\\frac{\partial^2 f(x,y,z)}{\partial^2 z}={f(x,y,z+1)-2f(x,y,z)+f(x,y,z-1)}\\\frac{\partial^2 f(x,y,z)}{\partial x\partial y}=\frac{f(x+1,y+1,z)-f(x-1,y+1,z)-f(x+1,y-1,z)+f(x-1,y-1,z)}{4}\\\frac{\partial^2 f(x,y,z)}{\partial y\partial z}=\frac{f(x,y+1,z+1)-f(x,y-1,z+1)-f(x,y+1,z-1)+f(x,y-1,z-1)}{4}\\\frac{\partial^2 f(x,y,z)}{\partial x\partial z}=\frac{f(x+1,y,z+1)-f(x-1,y,z+1)-f(x+1,y,z-1)+f(x-1,y,z-1)}{4}$

Для устранения шума перед вычислением производных необходимо выполнить сглаживание Гаусса с некоторым значением $\sigma$ и только потом вычислять матрицу Гессе в конкретной точке уже по сглаженным значениям интенсивности. Сглаживание Гаусса в точке $(x_{0},y_{0},z_{0})$ :

$G(x_{0},y_{0},z_{0})=\sum_{i=-r}^{r}\sum_{j=-r}^{r}\sum_{k=-r}^{r}I(x_{0}+i,y_{0}+j,z_{0}+k)\cdot\frac{1}{(2\pi)^\frac{3}{2}\sigma^3}\cdot{e}^{-\frac{i^2+j^2+k^2}{2\sigma^2}}$

где I(x,y,z) возвращает значение интенсивности в точке (x,y,z) ; обычно принимает значение $\lceil{3\sigma}\rceil$ ( $\lceil{\cdot}\rceil$ — округление в большую сторону); $\sigma$ — коэффициент сглаживания.

Значения интенсивности в перпендикулярном срезе сосуда до сглаживания, после сглаживания, на третьем рисунке значения сосудистости:

При этом, чем большее значение $\sigma$ мы используем, тем более толстые сосуды мы пробуем найти. А поскольку нам нужны вообще все сосуды — и тонкие, и толстые, то нам нужно много значений $\sigma$ , и для каждого значения нам следует выполнить сглаживание. Мы остановились на диапазоне $\sigma$ , который включает в себя 10 значений.

Однако, с ростом $\sigma$ уменьшаются значениях вторых производных. Это, в свою очередь, приводит к тому, что сосудистость, рассчитанная с большой $\sigma$ , оказывается меньше сосудистости, рассчитанной в той же точке с маленькой $\sigma$ , даже если в действительности мы имеем дело с толстым сосудом. Получается, что независимо от реальной картины всегда преобладают структуры, напоминающие тонкие сосуды. Поэтому возникает вопрос: как соотнести друг с другом все полученные сосудистости для каждой из $\sigma$ ? Для этого необходимо выполнить нормализацию между результатами вычислений. В литературе обычно проводятся манипуляции с полученной сосудистостью, например:

$V=V_{\sigma}e^{k\sigma}$ , где — нормализованная сосудистость, $V_{\sigma}$ — сосудистость, полученная для сглаживания с $\sigma$ , подбирается экспериментально, хороший вариант 0.5.

Мы же воспользовались так называемой гамма-нормализацией: $H_{M}'=H_{M}\sigma^{2d\gamma}$ , где $H_{M}$ — матрица Гессе, — порядок производной, т.е. 2, $\gamma$ подбирается экспериментально, в нашем случае хорошо себя показал вариант 0.5. Тогда вся формула сводится к $H_{M}'=H_{M}\sigma^2$ .

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

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

- выбрать значение $\sigma$ и провести сглаживание в точке и в её соседях (они нужны для вычисления производных),
- посчитать вторые производные и построить матрицу Гессе,
- нормализовать матрицу Гессе, домножив ее на $\sigma^2$ , и найти собственные значения,
- вычислить сосудистость,
- если остались незадействованные $\sigma$ , то вернуться в начало, в противном случае переход в пункт 2.
Рассчитать собственные векторы матрицы Гессе, для которой сосудистость оказалась максимальной, и получить вектор-направление сосуда.

Совмещая этот алгоритм с алгоритмом построения центральной линии мы получим финальный результат:

Оптимизация

Описанный выше подход будет работать прекрасно, но есть несколько моментов. Во-первых, для повышения производительности сглаживание необходимо заранее выполнять сразу для всех вокселей и для всех $\sigma$ . Во-вторых, с учетом того, что исследования обычно имеют размер примерно 512×512х512 вокселей, размер памяти, требуемой для хранения результатов сглаживания, в среднем будет занимать около 5 Гб. Чтобы сократить количество расходуемой памяти мы воспользовались пирамидой (scale space pyramid).

Идея заключается в том, что раз уж после каждого сглаживания значения интенсивности в соседних вокселях размываются и становятся примерно равными, то нет смысла хранить их все. Т.е. чем больше $\sigma$ , тем меньше нам потом потребуется просчитанных вокселей, чтобы восстановить сглаживание по всему объему.

В общем случае работает это так. Нулевой слой объема является оригинальным. Вдоль каждой из осей объема несколько раз применяется сглаживающий фильтр (¼, ½, ¼), и после размерность объема уменьшается в два раза. Полученный объем будет принадлежать первому слою. Затем операция повторяется и получается второй слой. Потом третий и т.д. Каждому слою соответствует определенное значение $\sigma$ . Пример для 2D:

Нетрудно посчитать, что в 3D суммарное количество памяти будет равно $\frac{8}{7}N$ , где — количество памяти, требуемое для хранения оригинального объема. Однако, использование пирамиды таит в себе очень много трудностей и проблем, с которыми мы столкнулись, и которые тянут на отдельную статью, если про них рассказывать.

Итог

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

Для КТ исследований не удается обнаружить сосуды, которые проходят внутри костей (как, например, в позвоночнике)
Не всегда удается обнаружить сосуды, если они проходят вблизи вытянутых продолговатых костей, тканей или структур (в таких случаях кости, ткани и структуры сами могут ошибочно приниматься за сосуды)
Узким местом являются бифуркации (раздвоения сосудов)

Проблемы 1 и 2 решаются путем вычитания костей. Для этого необходимо два КТ исследования — с контрастом и без. Понятно, что единственным отличием таких исследований будут являться подсвеченные сосуды. Другими словами, если вычесть интенсивность каждого вокселя исследования без сосудов из интенсивности каждого вокселя исследования с сосудами, то в результате ненулевую интенсивность будут иметь только воксели, относящиеся к сосудам. Но поскольку два исследования нельзя сделать с абсолютно одинаковым положением тела пациента, то главной проблемой становится совмещение двух исследований в 3D-пространстве. Для этого используется трансформация поворотом и смещением. Ориентация в пространстве идет по костям, т.к. они представляют собой жесткие структуры. Для нахождения костей мы воспользовались очень интересным алгоритмом водопадов, основанных на графах (waterfall based on graphs).

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

На этом все, спасибо за внимание. На всякий случай ссылка на продукт.