Алгоритм подавления полос на изображении как инструмент улучшения качества томографической реконструкции

    Вернемся к томографии, которой у нас в Smart Engines уделено большое внимание. Сегодня мы расскажем про алгоритм уменьшения выраженности полос на изображении. Полосы на томографической синограмме никому бы не мешали, ведь синограммы не предъявляются врачам или другим пользователям томографов, но эти полосы приводят к появлению концентрических окружностей на восстановленных изображениях (слева на рисунке). Основным инструментом борьбы с полосами в предлагаемом алгоритме является операция ведомой фильтрации (Guided Filtering). Мы расскажем как построить ведущее изображение для синограммы, рассчитать скорректированную синограмму и использовать ее в процедуре томографической реконструкции, чтобы получить восстановленное изображение без кольцевых артефактов (справа на рисунке).


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

    В томографических установках часто есть выделенная точка, относительно которой что-то вращается: либо объект, закрепленный в держателе на гониометре, вращается, а источник и детектор неподвижны; либо вокруг выделенной точки вращается система «источник-детектор». Это два принципиально разных подхода к организации процедуры сбора томографических проекций, проблемы есть в обоих случаях. Так откуда же на восстановленном изображении возникают артефакты типа колец и как уменьшить степень их выраженности? Результат реконструкции (горизонтальное сечение пористого объекта с артефактами в виде концентрических окружностей) приведен на рис. 1.


    Рис. 1 Результат реконструкции без применения процедуры подавления колец [1]

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


    Рис. 2 Томографическая проекция. Красным цветом выделено горизонтальное сечение, участвующее в построении синограммы

    Построим новое изображение — синограмму, собрав соответствующие строчки всех угловых проекций (рис. 3). i-я строка полученного изображения соответствует i-му проекционному углу. Т.е. каждый столбец содержит измерения одной и той же ячейки, соответствующие разным проекционным углам. Синограммой такое изображение называется не случайно. Легко заметить, что в центральной области оно состоит из синусоид.


    Рис. 3 Участок синограммы, собранной из строчек томографических проекций

    На изображении, особенно на светлых левом и правом краях, где отсутствует тень объекта, просматриваются вертикальные полосы. Наличие вертикальных полос на синограмме является причиной возникновения концентрических окружностей на восстановленном изображении. Причин появления вертикальных полос может быть несколько. Разный отклик ячеек детектора на одинаковый поступающий сигнал — одна из них. Производители детекторов стараются провести компенсацию этого эффекта в момент выпуска детектора на рынок. Компенсировать деградацию, происходящую во время жизненного цикла прибора, может периодически обновляемая, так называемая, pixel map. Ее создание — процедура затратная, поскольку требует наличия откалиброванного источника. Т.е. либо пользователь должен иметь собственный такой источник, либо вынужден обращаться в компании, предоставляющие такого рода услуги. Альтернатива — использовать алгоритмы подавления вертикальных полос. Второй возможный источник появления полос на синограмме — это сшивка участков изображения. Дело в том, что томографируемый объект не всегда умещается целиком в поле вида детектора. Человечество неумолимо движется в сторону увеличения пространственного разрешения метода томографии. Хочется томографировать большие объекты, например, человеческую голову (размер по вертикали несколько десятков сантиметров), с НАНОМЕТРОВЫМ разрешением. Легко вычислить сколько пикселей должна иметь матрица, чтобы зарегистрировать желаемую проекцию. Сейчас задачу пытаются решать, сшивая зарегистрированные участки частей объекта, снятые с перекрытием. При сшивке возникают схожие артефакты. Еще один источник полос — нестабильность самого пучка, т.е. изменение интенсивности пучка от проекции к проекции. Какова бы ни была причина появления вертикальных полос, при реконструкции они порождают кольцевые артефакты, которые обычно убирают постпроцессингом реконструированных изображений. Мы же будем бороться с кольцами фильтрацией вертикальных полосок.

    Описание алгоритма и экспериентов


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


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

    В методе подавления вертикальных полос алгоритм Guided Filtering (рис. 5 [2]) используется как основной инструмент.



    Рис. 5 Принципиальная схема фильтрации [2]

    Основой Guided Filtering является наличие ведущего и ведомого изображений. Мы хотим построить ведущее изображение, на котором синусоиды будут проявлены хорошо, а выраженность вертикальных полос ослабнет. Первым шагом рассчитаем производную по горизонтальному направлению (рис. 6), т.е. по направлению, перпендикулярному направлению полос.


    Рис. 6 Производная по горизонтальному направлению от прологарифмированной синограммы

    На увеличенном участке изображения (рис. 7) шум, обусловленный нестабильностью пучка во время измерения, проявляется в виде разрывов вертикальных полос.


    Рис. 7 Увеличенный участок изображения рис. 6

    Для каждого столбца выполним операцию одномерной свертки, чтобы уменьшить вклад высокочастотной шумовой компоненты (рис. 8).


    Рис. 8 Результат свертки

    Увеличенный участок приведен на рис. 9.


    Рис. 9 Увеличенный участок изображения

    Однако, мы все еще на пути построения ведущего изображения. Применим построчно операцию кумулятивного суммирования к изображению, представленному на рис. 8. Результат увел нас из пространства производных, сохранив низкочастотные зависимости (рис. 10).


    Рис. 10 Результат построчного применения операции кумулятивного суммирования

    Вычтем полученный результат из логарифмированной синограммы, завершив процесс построения ведущего изображения (рис. 11). Нам осталось выполнить операцию фильтрации (рис. 5).


    Рис. 11 Ведущее изображение

    Результат выполнения операции с окном (9,1) и E=0.00001 приведен на рис. 12.


    Рис. 12 Результат выполненной операции фильтрации


    Рис. 13 Разница между входным изображением и результатом фильтрации

    На рис. 14 приведены результаты томографической реконструкции с использованием нефильтрованных (слева) и фильтрованных (справа) проекций.


    Рис. 13 Результаты томографической реконструкции

    Заключение


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

    Список используемых источников
    1. Е.Е. Берловская, А.В. Бузмаков, А.С. Ингачева и др. Алгоритм подавления ортотропных артефактов регистрации изображений в рентгеновском и терагерцовом диапазонах. // Информационные процессы. 2019, т.19, ном. 2, стр. 1-9.
    2. kaiminghe.com/eccv10/eccv10ppt.pdf
    Smart Engines
    Обработка изображений, распознавание в видеопотоке

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

      0
      Спасибо за статью! Я тоже как раз страдаю от полос на рентгеновском изображении ))

      Сейчас более-менее успешно использую один доморощенный алгоритм, хотя изобретать ничего не следовало, так как разных алгоритмов известно много — но все они требуют отнюдь не моего уровня понимания математики и работают явно небыстро. В вашей статье кажется предлагается что-то сравнительно производительное и кажется этого алгоритма я еще не видел. Тем не менее я не все понял, расскажите подробнее (для дураков):
      1) По какому основанию логарифмируете исходное изображение?
      2) Как именно вычисляется производная поперек полос?
      3) Как именно производится свертка вдоль полос?
      4) «Применим построчно операцию кумулятивного суммирования к изображению» — как именно? Кумулятивная сумма по строке — это же по логике просто сумма значений в строке, какой в ней смысл?
        0
        Просим прощения, промазали. Ответ — чуть ниже.
        0
        Добрый день.
        1) Без разницы, но чтобы легче было обсуждать константы – натуральный.
        2) По двухточечному паттерну, не центрированному.
        3) Без разницы. Наивной дискретной сверткой, например. Если вопрос на самом деле был о ядре, то оно – гауссово.
        4) Кумулятивной суммой в нашей школе называется дискретный аналог неопределенного интеграла. Это не полная сумма, а последовательность частичных сумм.
          0
          Спасибо за ответ!
          2) Ну, я догадался, что используется конечная разность, вопрос в том, куда записывается (Yi-Yi-1) — в i-й столбец, или в i+1-й? Ответ наверное очевиден (и может быть вообще нет разницы), но я решил переспросить.
          3) Ну то есть просто берете гауссово ядро и каждый пиксель им обрабатываете? Радиус подбирается опытным путем очевидно?
          4) Ну, да, кумулятивная сумма так и определяется. Я просто пытаюсь понять, как вы используете эти суммы. То есть в i-й пиксель вы записываете сумму все производных в той же строке до позиции i?
            0
            Добрый день.
            2) Важно, чтобы операции кумулятивной суммации и разностного дифференцирования были согласованы.
            3) Да.
            4) Да.
              0
              Спасибо! Попробую попробовать на своих рентгенограммах, посмотрю на результат.
          0

          Не упрек, а уточнение. За статью спасибо. Этот метод еще бы и на gpu перевести… Да в реальном времени.

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

            Стоит все же грамотно использовать термины:


            Для каждого столбца выполним операцию одномерной свертки, чтобы уменьшить вклад высокочастотной шумовой компоненты (рис. 8).

            Видимо, речь все же о фильтре высоких частот. Свертка может быть с разными ядрами и вовсе не все из них будут фильтром высоких частот! Выше в комментариях говорится про гауссово ядро, но и для гауссова ядра далеко не для всех параметров получится фильтр высокой частоты (в пределе можно и свертку с константой получить).


            Применим построчно операцию кумулятивного суммирования к изображению, представленному на рис. 8.

            На картинке слева направо значения нарастают и потом убывают (слева черная область, в центре белая и справа опять черная). Ну, конечно, это может быть кумулятивная сумма по убывающей радиальной координате… если судить по исходной картинке — на кумулятивную сумму по возрастанию горизонтальной координаты не похоже. Что подразумевается автором — остается догадываться.

              0

              А может, речь и вовсе про фильтр низких частот, чтобы максимально исключить высокочастотную компоненту? Сначала подумалось про полосовой фильтр в области высоких частот (раз написано «уменьшить вклад»). В общем, это не воспроизводимо.

                0
                Заметки на Хабр нами обычно выкладываются по материалам научных публикаций. При этом текст сокращается и упрощается. Цель публикаций на Хабре — познакомить широкий круг людей с идеями и подходами, поэтому многие детали имплементации опускаются. К сожалению, при этом в тексте иногда появляются явные огрехи. Спасибо, что указываете на них, это позволит нам стать лучше. Последний же Ваш комментарий гораздо важнее. Для тех, кто хочет воспроизвести наш алгоритм (мы даже не надеялись на такую удачу) мы рекомендуем обратиться к нашей исходной научной статье: www.jip.ru/2019/200-207-2019.pdf. Там, мы надеемся, мы написали все необходимые детали. Если же мы (и рецензенты) что-то упустили, обращайтесь, мы ничего не утаим.
                  0

                  О, отлично! Хорошо бы эту ссылку в статью добавить, а то там в списке материалов совсем не то...

                    0

                    Посмотрел статью — все сводится к вызову неведомой функции ImGuidedFilter, которая в OpenCV отсутствует и что и как она делает, в статье не сказано. Сказано только, что результат этой функции просто супер. Мда.

                      0
                        0

                        Я тоже на эту функцию посмотрел, но название отличается ImGuidedFilter != guidedFilter и функция не из ядра OpenCV, а из дополнительных модулей — обо всем этом принято явно говорить в статьях, если уж зачем-то нужно переименовать библиотечную функцию и использовать дополнительные модули. На просторах интернета для OpenCV есть еще несколько отличающихся реализаций guidedFilter, кстати. В матлабе есть одноименная функция ImGuidedFilter, может, код у авторов вообще на матлабе написан, а для статьи несколько строчек переписали на питоне...

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

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