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


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


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



Содержание



Строение дыхательной системы


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


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



Шкала Хаунсфилда


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



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


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


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


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


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


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


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


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


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


Математическая морфология


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


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


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


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


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


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


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



Алгоритм Ли и RLE-компрессия


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


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



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


Пороговое преобразование базового объема


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



Отсечение внешнего воздуха


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



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



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



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



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



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



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



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




Выделение максимального по объему объекта


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



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



Закрытие сосудов внутри лёгких


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



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


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


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


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


Реализация алгоритма в среде MATLAB


Метод getRespiratoryOrgans


% Возвращает весь объем органов дыхания (объем легких и дыхательных путей)
% без разделения левого и правого легкого.

% V = базовый объем с данными по радиоплотности в единицах Хаунсфилда.
% cr = радиус морфологического закрытия сосудов.
% ci = количество итераций морфологического закрытия сосудов (например, 3 раза
% делают дилатацию и после этого 3 раза делают эрозию.

function RO = getRespiratoryOrgans(V,cr,ci)
% Пороговое преобразование базового объёма
% по порогу < -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


% Возвращает объем, связанный с краевыми поверхностями трехмерного объема
% (кроме верхней поверхности, поскольку легкие могут иметь
% соединение с верхней поверхностью).

% EAL = эродированный бинаризованный объем легких и воздуха.

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


% Возвращает самый большой объект в трехмерном объеме "V".

% O = воксели самого большого объекта.
% m = объем самого большого объекта.

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).