Решение проблемы обнаружения центральной линии сосуда

Суть задачи


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


Одна из задач врача — найти стенозы и оценить, насколько они опасны. Задача же разработчика, как обычно, облегчить работу конечного пользователя. Для этого необходимо построить полную 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}&lt;0,\lambda_{3}&lt;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. Определяем направление сосуда в текущем вокселе,
  2. Делаем перпендикулярный направлению срез сосуда и перемещаемся по градиенту в воксель среза с максимальной сосудистостью,
  3. Определяем направление сосуда в вокселе с максимальной сосудистостью,
  4. Делаем небольшой шаг в направлении сосуда и попадаем в новый воксель,
  5. Возвращаемся в шаг 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}$

где B и C имеют предопределенные значения. Если 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); r обычно принимает значение \lceil{3\sigma}\rceil (\lceil{\cdot}\rceil — округление в большую сторону); \sigma — коэффициент сглаживания.

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


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

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

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

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

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

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

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

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



Оптимизация


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

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

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

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

Итог


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

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

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

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

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

Inobitec

101,00

Расширяем границы видимого и возможного для врачей

Поделиться публикацией

Похожие публикации

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

    Не вполне понимаю, каким образом сократили расход памяти? Неужели древообразная ссылочная структура с листовыми узлами на разных уровнях? Мне это не кажется хорошим решением — ссылки сами по себе занимают объём, дико фрагментируют память и т.д.
    Мы разрабатывали т.н. "блочные модели" для геологии (по сути, тот же 3D массив). Расход памяти там сокращается за счёт своей системы поблочного свопинга и сжатия данных быстрыми алгоритмами. Пирамида только ускоряет обработку данных, но не особо сокращает расход памяти.


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

    Есть такое) Не так страшна реализация пирамиды, как реализация алгоритмов на ней.

      +1
      Не вполне понимаю, каким образом сократили расход памяти?
      Пирамида нужна только для хранения результатов сглаживания. Грубо говоря, если сделать в лоб, то нужно памяти: sizeof(float)*512x512x512 (для сглаживания с sigma1) + sizeof(float)*512x512x512 (для сглаживания с sigma2) +… + sizeof(float)*512x512x512 (для сглаживания с sigma10). Ответ там же:
      Идея заключается в том, что раз уж после каждого сглаживания значения интенсивности в соседних вокселях размываются и становятся примерно равными, то нет смысла хранить их все. Т.е. чем больше sigma, тем меньше нам потом потребуется просчитанных вокселей, чтобы восстановить сглаживание по всему объему.
      В самой пирамиде все хранится по аналогии с рисунком без всяких тонкостей: для sigma1 512х512х512, для sigma2 256х256х256, для sigma3 128х128х128 и т.д. Т.е. сокращение памяти идет за счет того, что само сглаживание имеет некоторые особенности. Scale space pyramid — довольно стандартное решение, можно посмотреть в интернете/википедии. Сжатие самих данных — это уже следующий шаг, но он не потребовался.
        0

        В компьютерной графике Scale Space Pyramid традиционно называется MIP maps.

          0
          Мипмапы используются для избегания шума, который возникает, если отрисовать текстуру в уменьшенном виде «пиксель-в-пиксель» на экране. В нашем случае пирамида нужна, чтобы потом получаемые значения производных после сглаживания максимально были приближены к этим же значениям, но получаемым без использования пирамид. Тут никак не участвуют расстояние до камеры или разрешение экрана, но зато нужна точность при вычислении производных, что привносит свои тонкости в построение пирамиды. Хотите, проведите аналогию, но, как мне кажется, она будет довольно грубой.
          0

          А, ясно. Вы сокращаете память с n*sizeof(float)*512^3 (n>1, n in Z) до (8/7)*sizeof(float)*512^3. Я подумал, что с sizeof(float)*512^3 до величины, меньшей этого значения.

        +3
        Уж не знаю, как вы так нагрешили раньше, но конкретно эта статья — огонь! Обязательно пишите еще.
          0
          А теперь все то же самое на бис для ультразвука :) CT/MRT конечно проще…
            –1
            Делаем небольшой шаг в направлении сосуда и попадаем в новый воксель

            Если пользователь отмечает две точки сосуда, то вы можете методом поиска кратчайшего пути найти весь сосуд между ними. Это ещё решит проблему с "разрывами" сосудов. Мы подобное делали на сходных объёмах данных, правда в 2D изображениях.
            Может быть будет возможно отказаться от сглаживания гауссианами, я так понимаю, что сосудистость резко возрастает к центру, а алгоритму поиска пути только и нужно, чтобы сосуд резко выделялся на фоне. Тогда линия пойдёт ровно по центру.

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


              Индекс похожести на дорогу (чем чернее, тем лучше)


              Увеличенный фрагмент дороги


                0
                Ух ты ж елки-палки, вот где специалисты-то водятся! Я вообще-то сервис-инженер, в частности по исследовательским рентгеновским томографам для мелких лабораторных животных. От человеческих томографов наши отличаются разрешением в единицы микрон на воксель и небольшим сканируемым объемом (труба 80мм диаметром). Иногда пользователи вопрос задают типа «а в какой программе можно мышкой щелкнуть и весь скелет автоматически выделить», да и самому иногда результаты обрабатывать приходится. А программа-то вот она, и не буржуйская :) Сколько оно у вас примерно стоит? А то исследователи обычно небогатые, по сравнению с медиками.
                  0
                  У нас несколько версий. Одна из них вообще бесплатная, но функционал также ограничен. Т.к. я не в курсе всех тонкостей, то лучше всего написать на почту (Наш сайт->Компания->Контакты). Ее почти круглосуточно мониторит специальный человек, и он с радостью отвечает на любые вопросы. В идеале надо подробно написать какие именно функционал и возможности вам нужны. Вариант два: просто скачать версию с сайта (она будет полноценной в течение месяца), а когда кончится триальный период, точно также спросить на почте «сколько будет стоить версия с этим и вот тем?».
                  0
                  А нельзя адаптировать под ваши нужды вот эту методу?
                  houdinigubbins.wordpress.com/2017/07/22/fiedler-vector-and-the-medial-axis
                    0
                    Нет, там на входе сетка и центральная линия ищется за счет контуров. С нашими входными данными делается ровно наоборот: только после построения центральной линии можно найти контуры (стенки) сосуда. Также задача найти контур в объеме нетривиальна и нормально не решается через классический edge detection.
                    0
                    Замечательная статья. Метод знаком, но с технической стороны, без названия. Не могли бы Вы дать источник метода?
                      0
                      Спасибо! У метода нет общепринятого названия, как вариант: Multi-Scale Vessel Segmentation Using Hessian Matrix Enhancement, ключевые слова Vessel Hessian Matrix. Полноценного источника также не существует: везде дается теория, типа «сглаживайте и стройте матрицу Гессе», но о конкретных реализациях или умалчивают, или они различаются. В итоге нужно перелопатить 100500 сайтов и pdf-фок, чтобы найти рабочее решение для каждого этапа вычислений.

                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                    Самое читаемое