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

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

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

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

_bau76jpkcf03bu-lwhmziqw-4c.jpeg

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

yxjuspnur3gwm3ghxhqv50upj9i.png

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

jvblho_x-0woitx02gbfjuxxhju.jpeg

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

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

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

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

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

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

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


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

n5o95wjj7avpy456i33a4swu_y4.png

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

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

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

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

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

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

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

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

bvkdfvdijxabpqjpbgrlw5pwdz4.gif7wfth0jdj6img3nvemmktkysboq.gif

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

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

ifkrh6br6sgi65g8sowexfxxsis.gif9frrncdin9axwhdwgad8vndnlqy.gif

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

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

dilvssyv2do2rvjlbb3xov798rw.gif

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

c_esdje-omegb-i31quh2tslaoi.png

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

mum1bidj1at0nvceitxh5jtzqs4.png

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

fuzb3bzxnvxvofjz-mxoelly5w4.gif

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

ufjambvoavwiks-09hms7wc02lg.png

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

furq26oo1ozjtnxrzu1ahknvjk4.gif

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

yji7nhiveqpxbky3_qgj9b6n0zq.png

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

vntunsho-og3pm1at7dqi3_ltna.gif

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

frhxdfz7ax9fz-ps7x9zjrlgiqo.png

64hd2ivbfu5g-d2yjjqkx5kv8pe.png

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

7lgtbdfvf0st73wqzlacvax1r-0.png5dce8urp2qvopyat905ex9y_qqs.png

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

jrmim2de5siht6zo2nhoz8l1e5e.png

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

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

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


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

0zd50js0yyqxagy_sgsp2wjyu8o.png


Метод 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).

© Habrahabr.ru