Хабрахабр

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

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

Также я предполагал, что морфологические операции – это только подготовка изображения к более сложным алгоритмам. Я предполагал, что без нейронной сети удастся получить точность не выше 70%. Но в результате обработки тех, хоть и немногочисленных 40 образцов томографических данных, что есть на руках, алгоритм выделил легкие без ошибок, причём после теста на первых пяти случаях алгоритм уже не претерпевал значительных изменений и с первого применения правильно отработал на остальных 35 исследованиях без изменения настроек.

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

Выделяют верхний и нижний отделы дыхательных путей. Дыхательная система включает в себя дыхательные пути и лёгкие. Всё, что выше гортани – верхний отдел, а остальное – нижний.
Перечислим дыхательные органы:
Носовая полость: — нос, гайморовы пазухи и т.д.
Глотка — канал, по которому перемещается пища и воздух.
Гортань – отвечает за образование голоса. Точкой разделения между нижними и верхними дыхательными путями является точка пересечения пищевых и дыхательных каналов. Находится на уровне шейных позвонков C4-C6.
Трахея – трубка, соединяющая гортань и бронхи.
Бронхи – дыхательные каналы, основная часть которых находится внутри лёгких.
Лёгкие – основной дыхательный орган.

Годфри Хаунсфилд — британский инженер-электрик, который вместе с американским теоретиком Алланом Кормаком разработал компьютерную томографию, за что получил Нобелевскую премию в 1979 году.

Шкала Хаунсфилда — количественная шкала рентгеновской плотности, которая измеряется в единицах Хаунсфилда, обозначаемых HU.

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

Рентгеновская плотность вычисляется по формуле:

$$display$$-μ_{water} \over μ_{water}-μ_{air}} \times 1000$$display$$

где $μ_X, μ_{water}, μ_{air}$ — линейные коэффициенты ослабления для измеряемого вещества, воды и воздуха.

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

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

  • Воздух: -1000 HU.
  • Дыхательные органы: от -950 до -300 HU.
  • Кровь (без контрастирования сосудов): от 0 до 100 HU.
  • Кости: от 100 до 1000 HU.

Ссылки на Википедию: шкала Хаунсфилда, Годфри Хаунсфилд, коэффициент ослабления.

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

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

К основным морфологическим операциям относят:

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

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

Применяется для закрытия отверстий внутри объектов и для объединения рядом стоящих объектов. Морфологическое закрытие (closing) — это дилатация с последующей эрозией.

Применяется для удаления мелких шумовых объектов и для разделения объектов на несколько объектов. Морфологическое открытие (opening) — это эрозия с последующей дилатацией.

Данный алгоритм изначально был придуман для поиска кратчайшего пути. Для выделения объектов в бинаризованном воксельном объёме используется алгоритм Ли. Его суть состоит в параллельном перемещении во всех возможных направлениях из начальной точки. Но мы используем его для выделения и перемещения объектов из одного трехмерного массива вокселей в другой. Для трехмерного случая возможны 26 либо 6 направлений движения из заданного вокселя (если воксель не находится с краю изображения).

Его суть заключается в том, что последовательности единиц и нулей заменяются цифрой, равной количеству элементов в последовательности. Для оптимизации по быстродействию был применен алгоритм кодирования длин серий (run-length encoding). Это позволяет уменьшить количество обращений к памяти. Например, строка “00110111” может быть заменена как: “2;2;1;3”.

Ссылки на Википедию: алгоритм Ли, алгоритм RLE.

Воксели воздуха имеют рентгеновскую плотность в промежутке от -1100 до -900 HU, а воксели дыхательных органов от -900 до -300 HU. С помощью томографа получены данные о рентгеновской плотности в каждой точке пространства. В итоге получим бинаризованный воксельный объём, содержащий только дыхательные органы и воздух. Поэтому можем убрать все лишние воксели, имеющие рентгеновскую плотность больше -300 HU.

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

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

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

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

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

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

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

Далее вычитаем внешний воздух из всего воздуха и дыхательных органов и получаем внутренний воздух и дыхательные органы.

Дыхательные органы — это отдельный объект. Далее выделим дыхательные органы как максимальный по объему объект. Связи между легкими и воздухом внутри желудочно-кишечного тракта нет.

Иначе в некоторых случаях может не оказаться связи между двумя лёгкими в результате низкого разрешения. Стоит заметить, что важен правильный выбор порога рентгеновской плотности на начальном шаге порогового преобразования. Поэтому следует повысить порог до -300 HU. Например, если считать, что воксели дыхательных органов имеют рентгеновскую плотность от -500 HU и менее, то в случае, приведённом ниже, выделение дыхательных органов как крупнейшего по объёму объекта приведёт к ошибке, так как отсутствует связь между двумя лёгкими.

Рентгеновская плотность сосудов составляет около -100.. Для захвата сосудов внутри лёгких применим морфологическое закрытие, то есть дилатацию с последующей эрозией с тем же радиусом. 100 HU.

Но в этом и нет необходимости.
Крупные кровеносные пути не закрылись. Цель данной операции была в уничтожении множества мелких отверстий внутри легких для упрощения дальнейшей сегментации легких.

В итоге получаем следующий алгоритм сегментации дыхательных органов:

  1. Пороговое преобразование базового объёма по порогу < -300 HU.
  2. Морфологическая эрозия радиусом 3 мм для разделения внешнего и внутреннего воздуха.
  3. Выделение внешнего воздуха по признаку соседства с граничными боковыми плоскостями воксельной сцены.
  4. Морфологическая дилатация внешнего воздуха для восстановления потерянной в результате эрозии информации.
  5. Вычитание внешнего воздуха из всего воздуха и дыхательных органов для получения внутреннего воздуха и дыхательных органов.
  6. Выделение максимального по объёму объекта.
  7. Морфологическое закрытие сосудов внутри лёгких.

Метод getRespiratoryOrgans

function RO = getRespiratoryOrgans(V,cr,ci)
% Загрузка из файла 3D матрицы со значениями
% рентгеновской плотности.
matFilePath=pwd+"/test_volumes/volume.mat";
load(matFilePath,'V');
% Пороговое преобразование базового объёма
% по порогу < -300 HU.
AL=~imbinarize(V,-300);
% Морфологическая эрозия радиусом 3 мм для
% разделения внешнего и внутреннего воздуха.
SE=strel('sphere',3);
EAL=imerode(AL,SE);
% Выделение внешнего воздуха по признаку соседства
% с граничными боковыми плоскостями воксельной сцены.
EA=getExternalAir(EAL);
% Морфологическая дилатация внешнего воздуха для
% восстановления потерянной в результате эрозии информации.
DEA=EA;
for i=1:4 DEA=imdilate(DEA,SE); DEA=DEA&AL;
end
% Вычитание внешнего воздуха из всего воздуха и дыхательных
% органов для получения внутреннего воздуха и дыхательных органов.
IAL=AL-DEA;
% Выделение максимального по объёму объекта.
RO=getMaxObject(IAL);
% Морфологическое закрытие сосудов внутри лёгких.
RO=closeVoxelVolume(RO,3,2);

Метод getExternalAir

function EA = getExternalAir(EAL)
% Функция bwlabeln сегментирует объекты: воксели одного
% объекта приравнивает к единице, другого – к двойке и т.д.
V=bwlabeln(EAL);
% Запрашиваем такие характеристики объектов, как ограничительный бокс
% и список всех вокселей объекта.
R=regionprops3(V,'BoundingBox','VoxelList');
n=height(R);
% Создаём 3-D матрицу для хранения вокселей внешнего воздуха.
s=size(EAL);
EA=zeros(s,'logical');
% Произведём перебор всех найденных объектов в цикле
% с целью нахождения объектов, принадлежащих внешнему воздуху.
for i=1:n % Определим координаты x и y, принадлежащие самым крайним % вокселям объекта. x0=R(i,1).BoundingBox(1); y0=R(i,1).BoundingBox(2); x1=x0+R(i,1).BoundingBox(4); y1=y0+R(i,1).BoundingBox(5); % Если крайние воксели объекта соприкасаются с боковыми % плоскостями сцены, то копируем все воксели данного объекта % в матрицу EA. if (x0 < 1 || x1 > s(1)-1 || y0 < 1 || y1 > s(2)-1) % Преобразуем данные о координатах вокселей объекта к % матричному типу: [[x1 y1 z1][x2 y2 z3] … [xn yn zn]]. mat=cell2mat(R(i,2).VoxelList); ms=size(mat); % Заполняем матрицу, содержащую воксели внешнего воздуха. for j=1:ms(1) x=mat(j,2); y=mat(j,1); z=mat(j,3); EA(x,y,z)=1; end end
end

Метод getMaxObject

### Метод getMaxObject
function [O,m] = getMaxObject(V)
% Сегментируем объекты. V=bwlabeln(V);
% Запрашиваем информацию об объёме объектов и координатах
% вокселей объектов.
R=regionprops3(V,'Volume','VoxelList');
% Определяем индекс максимального по объёму объекта.
v=R(:,1).Volume;
[m,i]=max(v);
% Создаём 3-D матрицу для хранения вокселей крупнейшего
% объекта.
s=size(V);
O=zeros(s,'logical');
% Переносим воксели крупнейшего объекта в новую матрицу. mat=cell2mat(R(i,2).VoxelList);
ms=size(mat);
for j=1:ms(1) x=mat(j,2); y=mat(j,1); z=mat(j,3); O(x,y,z)=1;
end

Исходный код можно скачать по ссылке.

Следующими статьями планируются:

  1. сегментация трахеи и бронхов;
  2. сегментация легких;
  3. сегментация долей легких.

Будут рассматриваться такие алгоритмы, как:

  1. дистанционное преобразование (distance transform);
  2. преобразование ближайших соседей (nearest neighbor transform, также известный как feature transform);
  3. вычисление собственных значений матрицы Гессе для сегментации плоских 3D объектов;
  4. сегментация методом водораздела (watershed segmentation).
Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть