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

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


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


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



    Содержание



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


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


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



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


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



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


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


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


    $$display$${μ_{X}-μ_{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) — это эрозия с последующей дилатацией. Применяется для удаления мелких шумовых объектов и для разделения объектов на несколько объектов.



    Алгоритм Ли и 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).
    • +20
    • 4,7k
    • 8
    Inobitec
    32,65
    Расширяем границы видимого и возможного для врачей
    Поделиться публикацией

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

    Комментарии 8

      +1
      Реально?
      В 2019 году, вместо того чтобы подготовить минимальную базу, за минимальные деньги, вы тратите время инженеров чтобы они вам нагенерили hand-craft признаков которые в 100% будут работать хуже чем нейронка?
        +2
        минимальная база и минимальные деньги это сколько на Ваш взгляд?
          0
          Конкретно по этой задаче, которую решает тут автор. Мне кажется, что для сопоставимых точностей хватит 50-100 КТшек. Скорее всего, даже 20-30.
          Разметку скорее всего проще сделать будет через набор срезов. Скорее всего на человека хватит 20-30 срезов. Остальное можно будет заинтерполировать. Для базы в 20 человек это 400 картинок на разметку. Разметить одну картинку ~30 секунд, скорее ощутимо меньше. Это ~3.5 часа на базу в 20 человек, на которой уже может начать работать.
          Разметчик по цене стоит в 4-5 раза дешевле программиста более-менее адекватного.

          Конечно, тут ещё должно быть время программиста выделено на написание софтины для разметки (на коленке и если данные 3д понятные — 1 день, добротно и если с данными в первый раз сталкиваться — неделя). Обучение тоже сколько-то съест. Но реально тут надо будет максимум U-net обучить для первой точности. Это пара дней для человека который разбирается в теме.
          Итого, идя по этому пути, можно решение сопоставимого качества получить за неделю где-то. Где 1-2 дня будет разметка.

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

          А вообще, конечно, существует и много других способов и подходов такое решать…
            +1
            Разметить одну картинку ~30 секунд, скорее ощутимо меньше. Это ~3.5 часа на базу в 20 человек, на которой уже может начать работать.


            Читаем первое предложение поста:

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


            Вы уверены что вы на сечении человека с первого раза отличите трахею от пищевода? Это все-таки не котиков от птичек отделять. В Mechanical Turk не отдашь. И может оказаться что разметчик стоит не так уж и дешево.

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

            Разметку скорее всего проще сделать будет через набор срезов. Скорее всего на человека хватит 20-30 срезов. Остальное можно будет заинтерполировать.

            Нейросеть потом сможет увидеть пространственную связь между срезами? Или для нее это будет 30 никак не связанных картинок? И в какой плоскости вы предлагаете нарезать слои?
              0
              К нейросетям много вопросов:
              1. На основе чего сделано предположение, что будет достаточно максимум 100 случаев? Это не говоря уже о том, что исследования в общем случае получить не так просто.
              2. Мы выполняем сегментацию 3-хмерных воксельных данных, а не отдельных срезов(есть принципиальная разница между сегментацией именно в 3-х измерениях и сегментацией по плоским срезам). Из этого следует во-первых то, что объем данных возрастает в несколько сотен или тысячу раз, а насколько я знаю, для нейросетей это проблема — даже обычным картинкам уменьшают разрешение, во-вторых просто уже архитектура сети скорее всего должна быть другой — не получится взять готовое решение. Вообще как то особо не получается найти примеры, чтобы нейросеть работала с воксельными данными.
              3. Я бы не сказал, что нейросети хорошо справляется именно с задачей сегментации. Та же U-Net, судя по результатам, приводимых в статьях, особой точностью не блещет — границы объектов получаются далеки от идеальных.
              4. Какой объем будет занимать результат обучения нейросети? Это должно работать в десктопной программе, которая не тащит с собой никаких особо больших баз.
            0
            можно использовать написанный алгоритм как некое начальное приближение. Прогнать датасет на этом алгоритме, получить грубо размеченные данные — и дальше вручную отобрать хорошие примеры / доразметить плохие. А потом конечно учить нейросети.
            +1
            Почему просто не использовать transfer function? Мы знаем плотность дыхательных органов. Соответсвенно можем на всем участке задать максимальную прозрачность, а на участке с плотностью дыхательных органов дать максимальное значение.
            Я так же писал подобные программы
              0
              Задание передаточной функции в конечном итоге соответствует просто выбору диапазона(ов) плотности. Как видно из текста статьи, задание диапазонов/порогов используется в качестве отдельных этапов алгоритма и целиком задачу не решает.

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

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