EA017302B1 - Способ подавления шума серий цифровых рентгенограмм - Google Patents

Способ подавления шума серий цифровых рентгенограмм Download PDF

Info

Publication number
EA017302B1
EA017302B1 EA201101502A EA201101502A EA017302B1 EA 017302 B1 EA017302 B1 EA 017302B1 EA 201101502 A EA201101502 A EA 201101502A EA 201101502 A EA201101502 A EA 201101502A EA 017302 B1 EA017302 B1 EA 017302B1
Authority
EA
Eurasian Patent Office
Prior art keywords
frame
noise
motion
estimation
image
Prior art date
Application number
EA201101502A
Other languages
English (en)
Other versions
EA201101502A1 (ru
Inventor
Сергей Васильевич МЕРКУРЬЕВ
Original Assignee
Закрытое Акционерное Общество "Импульс"
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Закрытое Акционерное Общество "Импульс" filed Critical Закрытое Акционерное Общество "Импульс"
Priority to EA201101502A priority Critical patent/EA017302B1/ru
Priority to EP12186390.6A priority patent/EP2579207A3/en
Priority to CN2012103719525A priority patent/CN103177423A/zh
Priority to JP2012216792A priority patent/JP2013127773A/ja
Priority to US13/645,518 priority patent/US20130089247A1/en
Priority to KR1020120110554A priority patent/KR20130038794A/ko
Publication of EA201101502A1 publication Critical patent/EA201101502A1/ru
Publication of EA017302B1 publication Critical patent/EA017302B1/ru

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10016—Video; Image sequence
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10116—X-ray image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20021—Dividing image into blocks, subimages or windows
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20036—Morphological image processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20201—Motion blur correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20216—Image averaging
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Picture Signal Circuits (AREA)

Abstract

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

Description

Область техники
Изобретение относится к области обработки цифровых изображений и может быть использовано для решения задач обработки цифровых изображений, полученных с помощью высокоэнергетического излучения, в том числе рентгеновского. Более конкретно, данное изобретение предназначено для подавления шума серий цифровых рентгенограмм.
Предшествующий уровень техники
В настоящее время в клинической практике применяются различные методы обработки серий цифровых рентгенограмм (медицинских изображений): повышение резкости, выделение анатомических структур, субтракция и др. К важнейшим методам улучшения качества серий медицинских изображений относятся методы подавления шумов. Так, в ангиографии в реальном времени осуществляют управление медицинским инструментарием и исследуют состояние сосудистой системы с помощью введения через катетер контрастного вещества. С целью уменьшения лучевой нагрузки на пациента и медперсонал применяют малые дозы облучения, что приводит к получению медицинских изображений с низким значением соотношения сигнал/шум [Скал с1. а1., 1таде 8сс.|испес ЕШсгшд ίη ОнапШт-ЫтНсй Νοίδο \νί11ι ЛррБсаίίοηδ ίο Бо\\-Оо5с Ииогоксору, ΙΕΕΕ Тгапк. Мсй1еа1 1тадшд, νοί. 12, Цкис 3, 610-621, 1993]. Поэтому в цифровой радиологии представляется актуальной задача разработки методов фильтрации шумов, способных сохранить диагностически значимую информацию, работать в реальном времени при высоких значениях разрешения, частоты смены кадров и существенном уровне шума, зависимого от интенсивности сигнала.
Можно выделить несколько проблем, с которыми сталкивается исследователь при разработке алгоритма фильтрации серий цифровых медицинских изображений. Первая проблема связана с оценкой шума изображений, уровень которого существенно зависит от интенсивности сигнала. Вторая проблема обусловлена необходимостью компенсации возможного резкого изменения интенсивности изображений, возникающего при использовании автоматической системы подстройки интенсивности или вследствие нестабильности рентгеновской трубки. Третья проблема связана с учётом изменений интенсивности последовательности кадров, вызванных движением: перемещением стола или системы трубка/детектор, смещением пациента или физиологическими процессами в его организме, током контрастного вещества в сосудах. Наконец, четвёртая проблема состоит в необходимости соблюдения баланса между качеством работы фильтра и требованием выполнения фильтрации в реальном времени. Данная проблема включает также выбор метода фильтрации: пространственно-временное усреднение в объёме стека изображений, временная фильтрация в пространстве преобразования (пирамидальное преобразование, вейвлетпреобразование и т.п.), комбинированные методы усреднения и т.п. Опишем каждую из этих проблем более подробно.
Рассмотрим проблему оценки шума, зависящего от интенсивности сигнала. Практически каждый качественный метод фильтрации шума, как одиночных цифровых изображений, так и их серий, использует информацию о спектральных свойствах шума. Весьма широкий класс методов шумоподавления использует в качестве параметра величину дисперсии шума. В публикациях разных авторов продемонстрировано, что чем точнее известно значение дисперсии шума, тем качественнее результат работы алгоритма шумоподавления [Воксо с1. а1, ИоЦс Ксйисйоп Гог СГа 1тадс Бспкога Εχρίοίΐίη^ Ηνδ ΒΟκινίοΓ. Бспкога 9(3), 1692-1713, 2009; Εοί, Ргас11са1 ОспоЦшд оГ СБррсй ог Оусгсхроксй Νοίδν 1тадс8, Ргос. 16111 Бигорсап 81§па1 Ргосскк. СопГ., Ευδ^Ο, 2008].
Рассмотрим природу шума изображений в цифровой рентгенографии. Детектор измеряет интенсивность ослабленного (прошедшего через объект исследования) рентгеновского излучения. Каждая ячейка детектора накапливает за время экспонирования в среднем Ν электронов путём поглощения фотонов. Количество Ν накопленных электронов можно моделировать распределённой по закону Пуассона случайной величиной
Р^ = я) = ехр(-1У)—, н>0.
п!
Случайные флуктуации числа поглощённых фотонов называют шумом фотонов.
В современных детекторах основным источником шума является шум фотонов. К дополнительным источникам шума относят шумы детекторной системы: шум чтения, тепловой шум, шумы усилителей, шумы квантования и другие [Ошо М., ИоЦс, Νοίδη ИоЦс]. Общий эффект данных источников шума моделируют распределённой по Гауссу случайной величиной [Воксо с1. а1., ИоЦс Ксйисйоп Гог СГа 1тадс Бспкогк Εχρίοίΐίπβ Ηνδ Всйауюг, Бспкога 9(3), 1692-1713, 2009; Εοί с1. а1., Ргасйса1 ро188ошап-даи881ап пощс тойсНпд апй Г1111пд Гог кшДс-ипадс гачг-йаЩ 1тадс Ргоссккшд, ΙΕΕΕ Тгапкасйопк оп 17, 1737-1754, 2008]. Согласно широко используемой модели при линейных электронных схемах дисперсия шума фотонов и дополнительных источников шумов линейно зависит от величины полезного сигнала [Яне, Цифровая обработка изображений, М.: Техносфера, 2007] σ2(/(ρ)) = «/(μ) + ί>, (1) где Ι(ρ) - уровень интенсивности сигнала в пикселе ρ = (χ, у) изображения.
Таким образом, шум цифровых изображений линейно зависит от интенсивности сигнала. Проблеме
- 1 017302 оценки шума, зависящего от сигнала, посвящено достаточно много публикаций [Воксо с1. а1., Νοίδο Вейисйои Гог С Га 1таде 8епкогк Εχρίοίΐίημ Ηνκ ВсЕауюг. 8епкогк 9(3), 1692-1713, 2009; Ροί еГ. а1., Ргасйса1 роккошап-даикыап поке тойеОпд апй ГООпд Гог ктдкйтаде гате-йаГа, 1таде Ргосеккшд, ΙΕΕΕ ТгапкасГ1опк оп 17, 1737-1754, 2008; РогкГпег, 1таде Ргергосеккшд Гог РеаГиге ЕхГгасйоп ίη Б|д0а1 1п1еп511у, Со1ог апй Вапде 1тадек, 1п 8рппдег ЬесГиге №1ек оп Εηγ11ι 8аепсек, 1998; Непке1 еГ. а1., ВоЬикГ апй РакГ Ебйтайоп оГ 81дпа1-БерепйепГ №ке ш Мей1са1 Х-Вау 1таде 8едиепсек, 8ргшдег, 2006; Ь1и еГ. а1., ЛиГотайс екйтайоп апй гето\га1 оГ поке Ггот а ктд1е 1таде, ΙΕΕΕ ТгапкасОопк оп Райет Апа1ук1к апй МасЫпе ШГеШдепсе, 30(2), 299-314, 2008; 8а1тей еГ. а1., 8щпа1-йерепйеп1 №ке С’11агас1епхайоп Гог МаттодгарЫс 1таде§ Бепокшд. ΙΜΕΚΟ ТС4 Зутрокшт (IМΕΚΟТС4 '08), Риен/е, Йа1у, 2008; \Уаед1Е ОггекйдаОопк шГо Г1е №ке СйагасГейкйск оГ Б|дЦ|/ей Лейа1 1таде8, 1п: 1пГ. Агсй. Рог РйоГодг. Лпй ВетоГе 8епкшд, Уо1. 32-2, 1998]. Так, в работе [Непке1 еГ. а1., ВоЬикГ апй РакГ ΕкГ^шаГ^οп оГ 8щпа1-Берепйеп1 №ке ш Мейюа1 Х-Вау 1таде 8едиепсек, 8ргшдег, 2006] описан непараметрический метод покадровой оценки шума серий цифровых медицинских изображений, при этом особый акцент делают на разработке алгоритма, способного работать в реальном времени. В публикации [Ро1 еГ. а1., Ргас11са1 роккошап-даиккап поке тойеОпд апй ГиОпд Гог ктдкйтаде гате-йаГа, 1таде Ргосеккшд, ΙΕΕΕ ТгапкасОопк оп 17, 1737-1754, 2008] рассматривают двухпараметрический подход к оценке шума цифровых снимков. Шум исходного цифрового изображения (полученного непосредственно с детектора и не прошедшего через нелинейные преобразования, такие как гамма-коррекция и т.п.) моделируют аддитивной по отношению к сигналу случайной величиной:
Л(р) = /(р) + сг2(7(р))«(р), (2) где 1п(р) - уровень интенсивности сигнала в пикселе р наблюдаемого зашумлённого изображения, а2(1(р)) - зависимость дисперсии шума от интенсивности сигнала вида (1), п(р)е^0,1) - стандартная нормальная случайная величина. Предлагают способ построения модельных кривых дисперсии шума, учитывающий нелинейности в работе сенсора, приводящие к недоэкспонированию и переэкспонированию, т.е. к нарушению линейного закона (1) на краях динамического диапазона.
Вторая проблема в разработке динамического фильтра заключается в необходимости учёта и компенсации изменения интенсивности фильтруемой серии изображений, т.н. мерцаний. Данная проблема подобна задаче компенсации фликкер-эффекта, возникающей как необходимый первый этап в цепочке различных методов обработки кадров видеопоследовательности [Бе1оп, Мо\зе апй У1йео 8са1е-Т1ше Ε^иа1^ζаГ^οη. ΙΕΕΕ ТгапкасОопк оп 1таде Ргосеккшд, 15(1):241-248, 2006; РШе еГ. а1., Α №\ν ВоЬикГ Тес1пк|ие Гог 81аЬШ/1пд ВйдЫпекк Р1исГиаГюпк ш Опаде 8едиепсе, Ш 2пй ^огккйор оп 81а0кОса1 МеГйойк ш У1йео Ргосеккшд Ш сопщпсОоп \νί11ι Ε^ν 2004. 8ргшдег, 2004; ^опд, Iшр^ονеά Рйскег Вешονа1 Тйгоидй МоОоп УесГогк Сотрепкайоп, Ш ТЫгй [п(егпа(1опа1 СопРегепсе оп Опаде апй ОгарЫск, рр. 552-555, 2004]. Для серии цифровых рентгенограмм мерцания могут быть вызваны использованием автоматической системы подстройки интенсивности, либо появится вследствие нестабильной работы рентгеновской трубки. Отсутствие учёта и компенсации мерцаний при временном усреднении серии изображений отрицательно сказывается на качестве работы фильтра. Так, временное усреднение стационарных по движению кадров без компенсации фликкер-эффекта приводит к существенному искажению оцениваемых значений текущего кадра. При наличии движения и его компенсации на основе методик типа поиска соответствия блоков пикселей (далее блоков), значительное изменение интенсивности изображений в серии приводит к грубым ошибкам в оценке движения и, как следствие, к плохой его компенсации.
Можно выделить две главные стороны в оценке и компенсации мерцаний. С одной стороны, в общем случае фликкер-эффект имеет нелокальную природу по пространству отдельно взятого кадра серии изображений. С другой стороны, он проявляется только при рассмотрении динамики процесса, т. е. нескольких подряд идущих кадров. Эти особенности составляют главную трудность в компенсации фликкер-эффекта, состоящую в определении того, что локальные изменения в контрасте/интенсивности кадра являются следствием движения или обусловлены прыгнувшей интенсивностью излучения рентгеновской трубки. Естественный способ моделирования фликкер-эффекта, широко описанный в литературе [Бе1оп, Мо\ае апй У1йео 8са1е-Т1те ΙΕΕΕ ТгапкасОопк оп Опаде Ргосеккшд, 15(1):241-248,
2006; Рй1е еГ. а1, Α №\ν ВоЬикГ Тес1шк|ие Гог 81аЫН/1пд ВпдЫпекк Р1исГиаГюпк ш Опаде 8едиепсе, Ш 2пй ХУогккОор оп 81а0кОса1 МеГйойк ш У1йео Ргосеккшд Ш сопщпсОоп νίίΠ Εί’ί’ν 2004. 8ргшдег, 2004; \Уопд еГ.а1, Iшр^ονеά Рйскег Вешονа1 Тйгоидй МоОоп УесГогк Сотрепкайоп, Ш ТЫгй !п1егпа11опа1 СопРегепсе оп Опаде апй ОгарЫск, рр. 552-555, 2004], основан на следующей формуле:
1(р)=А(р)10(р)+В(р), (3) где Е(р) - интенсивность исходного изображения в пикселе р; А(р), В(р) - гладкие, медленно изменяющиеся функции, определённые в области задания изображения; Ιφ) - интенсивность наблюдаемого изображения. Данная модель предполагает линейную зависимость интенсивности исходного изображения от интенсивности наблюдаемого изображения, а также учитывает отмеченную выше пространственную неоднородность фликкер-эффекта. Значительное количество подходов к оценке и компенсации мерцаний основано на разбиении изображения на сетку перекрывающихся блоков и поблочную оценку параметров функций А(р), В(р) с возможной отбраковкой выбросов, при этом для более аккуратной оценки параметры модели находят совместно с оценкой движения.
- 2 017302
Третья проблема, возникающая при разработке качественного метода фильтрации шумов серий изображений, заключается в необходимости учёта изменений в последовательности кадров, обусловленных движением. Движение является обязательным атрибутом любой диагностически значимой серии цифровых медицинских изображений. Усреднение нестационарных по движению кадров приводит к появлению артефактов размытия на результирующем изображении. На практике широко применяются два класса подходов к решению проблемы учёта движения при пространственно-временном усреднении кадров видеопоследовательности. Подходы первого типа основаны на детектировании регионов изображения, соответствующих движению и уменьшению силы временного подавления шума в этих регионах |Аи.Гг1сН(1д с1. а1., Х-Кау Ниогоксору 8райо-Тсшрога1 РШсппд \νί11ι 0Ь]сс1 Эс1се1юп. ΙΕΕΕ Тгаик. Мсй1са1 1шадшд 14 (1995) 733-746; Нси8с1 с!. а1., Мойои апй №15с ЭйссИои Гог 8райо-Тсшрога1 РШсппд оГ МсЙ1са1 Х-Кау 1шадс 8сс.|испсс5. Вюшсй, Тссй./Вюшсй. Еид. 50-1 (2005) 1106-1107; Коигай, Мойои Эйссйои аий ЕШшайои. 1и Воу1к, А.С., сй.: НаийЬоок оГ 1шадс аий У1йсо Ргосскщид. 2ий сй. ЕЕсуюг Асайсшк Ргскк (2005) 253-274]. Подходы второго типа оценивают траектории перемещения пикселей между последовательными изображениями, для чего осуществляют учёт и компенсацию движения [ВгаПсаи с! а1., Иикс Ксйисйои Е111сг8 Гог Эуиашю 1шадс 8сс.|испсс5: А гс\зс\\·. Ргос. ΙΕΕΕ, уо1. 83, рр. 1272-1292, 1995]. Стационарные по движению изображения можно явно сформировать после оценки движения в результате применения процедуры компенсация движения.
Простейшие из алгоритмов первого типа достаточно легко реализовать для выполнения в реальном времени. Для таких методов детектирование движения сводится, например, к пороговой обработке разности между текущим усредняемым кадром и ссылочным кадром. Если для некоторого пикселя текущего кадра данная разность выше определённого порогового значения, то временное усреднение заменяют, например, усреднением в пространстве текущего кадра. Величину порога естественно привязать к уровню шума, полученному на этапе оценки шума. Следует отметить, что подобные методы учёта движения страдают определёнными недостатками: появляется хорошо заметный импульсный шум, возникают негладкости на границах движущихся объектов. В свою очередь, алгоритмы второго типа, в которых производится оценка с последующей компенсацией движения (например, на основе поиска соответствия блоков в изображениях), считаются более качественными [Вгайсаи с! а1., ЫоЦс Ксйисйои П11сг5 Гог Эуиашю 1шадс 8сс.|испсс5: А гс\зс\\·. Ргос. ΙΕΕΕ, уо1. 83, рр. 1272-1292, 1995]. Однако есть определённые трудности и при их использовании: не все виды движения можно эффективно скомпенсировать [Яне, Цифровая обработка изображений, М.: Техносфера, 2007], кроме того, в силу относительно высокой вычислительной требовательности данных методов существуют определённые сложности при их реализации в реальном времени.
Четвёртая проблема состоит в необходимости обеспечения баланса между качеством работы фильтра и требованием выполнения в реальном времени. Необходимость выполнения обработки в реальном времени серии медицинских изображений при высокой частоте и большом размере кадров (например, 30 кадров в секунду при разрешении 512x512 пикселей, 15 кадров в секунду при разрешении 1024x1024 пикселей в кадре) на современных рентгенодиагностических комплексах вступает в определённое противоречие с желанием обеспечить высокое качество результата фильтрации. Такие операции, как оценка шума или оценка движения, сами по себе являются достаточно требовательными к вычислительным ресурсам даже с учётом возможностей современных компьютеров.
Раскрытие изобретения
Целью заявляемого изобретения является улучшение качества серий цифровых медицинских изображений (цифровых рентгенограмм).
Технический результат заявляемого способа заключается в повышении качества получаемого изображения за счет подавления шумов.
Технический результат в заявляемом способе подавления шума серий цифровых рентгенограмм, включающем получение серии цифровых рентгенограмм, покадровую оценку дисперсии шума, зависящего от сигнала, оценку и компенсацию движения, оценку и компенсацию фликкер-эффекта, рекурсивное усреднение кадров достигается тем, что при оценке дисперсии шума осуществляют морфологическое удаление значений пикселей снимка шума, соответствующих резким изменениям в исходном снимке; получают табличную зависимость шума от интенсивности сигнала с помощью робастной кусочнолинейной аппроксимации интервальных оценок дисперсии шума; вычисляют карту шума в виде попиксельной оценки дисперсии шума исходного цифрового изображения; при оценке и компенсации движения используют пирамидальную схему поиска соответствия блоков, блоки оценки движения рассматривают с перекрытием, при компенсации движения перекрытия блоков пикселей усредняют; учёт фликкерэффекта осуществляют делением ссылочного и текущего кадров на соответствующие поля средних значений, которые получают линейной НЧ-фильтрацией данных кадров; при компенсации фликкер-эффекта приводят локальную среднюю яркость ссылочного кадра к локальной средней яркости текущего, для чего делят ссылочный кадр на его скомпенсированное по движению поле средних значений и умножат на поле средних значений текущего кадра; при рекурсивном усреднении осуществляют коррекцию артефактов компенсации движения с помощью смешивания текущего отфильтрованного кадра с его исход
- 3 017302 ной версией на основе вычисления их пиксельного и структурного подобия.
Возможен вариант исполнения заявляемого способа, в котором при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка (отбрасывание границ), используют разность между текущим и скомпенсированным по движению ссылочным кадром.
Наиболее близким к используемому в настоящем изобретении методу оценки шума цифровых медицинских изображений (рентгенограмм) является способ, изложенный в публикации |Непке1 с1. а1., ВоЬик! апй Рак! Екйтайоп οί 81дпа1-Оерепйеп1 Νοίδο ίη Мейюа1 Х-Вау 1таде 8сс.|испес5. 8ртшдет, 2006], в соответствии с которым для построения по исходному изображению оценки зависящего от сигнала шума достаточно осуществить следующие этапы:
оценить полезный сигнал с помощью низкочастотной фильтрации исходного изображения и вычислить разность между исходным изображением и его оценкой, получив тем самым изображение шума;
отбросить тем или иным способом значения пикселей изображения шума, соответствующие резким изменениям (границы, одиночные горячие пиксели) в исходном изображении;
разбить диапазон интенсивностей оценочного изображения на интервалы и для каждого такого интервала накопить значения пикселей изображения шума, соответствующие пикселям оценочного изображения;
вычислить дисперсию шума в каждом интервале интенсивности по накопленным в данном интервале значениям пикселей изображения шума.
Настоящее изобретение использует изложенные в упомянутой выше публикации принципы оценки шума. Наиболее существенные отличия состоят в том, что для каждого кадра серии изображений выполняют следующие операции:
осуществляют удаление значений пикселей изображения шума, соответствующих резким изменениям в исходном изображении, с помощью морфологического выделения значений пикселей изображения шума, соответствующих границам на исходном изображении;
выполняют робастную локальную линейную аппроксимацию интервальных оценок дисперсии шума, в результате которой получают табличную функцию, описывающую зависимость шума от интенсивности сигнала;
вычисляют на основе оценочного изображения и найденной табличной функции карту шума в виде изображения, являющегося попиксельной оценкой дисперсии шума исходного цифрового изображения.
Наиболее близким к используемому в настоящем изобретении методу учёта и компенсации фликкер-эффекта является метод, описанный в публикации |\Уопд е1 а1., 1тртоуей Ейскет Вешоуа1 Тйгоидй Μοίίοη Уес1огк Сотрепкайоп. Ιη ТЫтй 1п1етпайопа1 СопГегепсе оп 1таде апй СтарЫск, рр. 552-555, 2004]. В этом методе мерцания моделируют гладкой функцией, мультипликативно искажающей значения интенсивностей кадров видеопоследовательности. Оценку данной функции предлагают осуществлять совместно с оценкой движения на основе поиска соответствия блоков, предлагают использовать при этом кусочно-постоянную искажающую функцию. Величину фликкер-эффекта в каждом блоке оценки движения текущего кадра полагают равной отношению средних значений этого блока и соответствующего блока ссылочного кадра, найденного в результате компенсации движения. Полученное кусочнопостоянное поле оценок фликкер-эффекта подвергают пороговой обработке и сглаживанию для удовлетворения условиям гладкости.
В настоящем изобретении фликкер-эффект также моделируют в виде мультипликативной гладкой функции. Однако, компенсацию фликкер-эффекта осуществляют на основе приведения локальных средних предыдущего кадра к локальным средним текущего кадра посредством деления скомпенсированного по движению предыдущего отфильтрованного кадра (ссылочного кадра) на скомпенсированное по движению его поле средних значений и умножением на поле средних значений текущего кадра. Поле средних значений кадра получают с помощью его усреднения скользящим средним с использованием окна, апертура которого равна размеру блока оценки движения. Соответствующим образом модифицируется и карта шума ссылочного изображения.
Движение в заявляемом изобретении оценивают и компенсируют на основе поиска соответствия блоков в изображениях серии |\Уапд е1. а1., Рак! Л1дотййтк Гог 111е Екйтайоп оГ Мойоп Уес1огк. 1ЕЕЕ Тгапк. 1таде Ртосеккшд, уо1. 8, по. 3, рр. 435-438, Матсй 1999]. Учёт фликкер-эффекта, оказывающего существенное влияние на качество оценки векторов движения, осуществляют на основе простейшего подхода, при котором на этапе оценки движения текущее и ссылочное изображения серии делят на их поля средних значений, что позволяет осуществить своего рода нормировку по яркости. Помимо учёта эффекта мерцаний интенсивности кадров оценка движения осуществляется с использованием пирамидальной схемы. Причина этого состоит в том, что при фильтрации в реальном времени серии медицинских изображений особую трудность представляет оценка и компенсация быстрых движений, например, в исследованиях сердца (коронарографии), когда используется сильное приближение к объекту исследования. В результате, для успешной компенсации движения, в особенности на низкой частоте смены кадров, область поиска подходящих блоков становится слишком велика. Поэтому для решения проблемы возрастающей вычислительной сложности при компенсации быстрых движений в изобретении используют
- 4 017302 пирамидальный алгоритм оценки движения, при котором текущий кадр и ссылочный кадр, который должен быть скомпенсирован по движению, раскладывают на несколько уровней в пирамиду изображений (не более трёх уровней, зависит от типа исследования);
оценку векторов перемещения блоков осуществляют на низкочастотных составляющих пирамиды поиском по небольшой области, при выборе радиуса которой учитывают число уровней в пирамиде;
при переходе с нижнего на верхний уровень пирамиды размер блока компенсации движения масштабируют и уточняют вектор движения данного блока поиском в подобласти малого радиуса.
Применение пирамидальной схемы позволяет существенно расширить область поиска подходящих блоков на ссылочном кадре. Для уменьшения эффекта блочности, типичного при компенсации движения на основе использования прямоугольных блоков, в изобретении блоки оценки движения рассматривают с перекрытиями [Отсйатб е1. а1., Отег1арреб В1оск Сотрепка!юп: Ап Екбтабоп-Тйеогебс Арртоасй. ΙΕΕΕ Ттапкасбопк Оп 1таде Ртосекктд, Уо1. 3, Ыо. 5, 8ер!етЬет 1994, рр. 693-699]. Типичный размер перекрытия составляет четверть радиуса блока оценки движения. При построении скомпенсированного по движению кадра пиксели, соответствующие областям перекрытия блоков, усредняют с использованием гладкой весовой функции.
Отличительной особенностью заявляемого изобретения является использование оригинальной методики компенсация артефактов алгоритма учета движения. Грубые ошибки компенсации движения, которые могут возникнуть из-за слишком малого радиуса области поиска, перекрытия объектов или резкой смены содержимого кадров, уменьшают с помощью смешивания отфильтрованного кадра с усреднённой в пространстве изображения его версией. Веса смешивания вычисляют на основе пиксельного и структурного подобий текущего изображения скомпенсированному по движению изображению. В качестве меры пиксельного подобия используют билатеральное расстояние [Тотаы е1. а1., Вйа1ета1 ЕШетшд Гог Сгау апб Со1ог 1тадек, 1п Ргос. 6111 ΙπΙ. СопГ. Сотри!ег Уыоп. №\ν Бе1Ы, 1пб1а, 1998, рр. 839-846]. Для вычисления структурного подобия применяют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров [Виабек е1. а1., Ыоп1оса1 1таде апб тоу1е бепоыпд. 1п!. 1. сотри!. У18юп, 76(2):123-139, 2008]. При вычислении финальной меры подобия кадров структурное и пиксельное подобия смешивают по линейному закону.
В настоящем изобретении с целью обеспечения высокого уровня подавления шума и удовлетворения требованию выполнения в реальном времени используют рекурсивный фильтр Калмана с автоматической оценкой зависящего от сигнала шума, учётом возможного резкого изменения интенсивности последовательных кадров, оценкой и компенсацией движения, уменьшением артефактов компенсации движения и временного пересглаживания. При разделении данных задач между многоядерным процессором обычного на сегодняшний день настольного компьютера и графической картой с технологией Ы\зб1а СИБА удалось получить результат, удовлетворяющий требованию фильтрации с высоким качеством в реальном времени.
Подробное описание изобретения
Заявляемое техническое решение, возможность его технической реализации и достижение указанного технического результата поясняется фиг. 1-4. На фиг. 1 показано устройство для получения рентгенограмм. На фиг. 2 показан график табличной функции, описывающей зависимость стандартного отклонения шума от интенсивности сигнала для кадра серии изображений сосудистого исследования, изображённого на фиг. 4. На фиг. 3 приведена карта стандартного отклонения шума изображения фиг. 4. На фиг. 4 показан исходный кадр серии изображений. На фиг. 5 приведён результат фильтрации данной серии изображений с 80-процентным уровнем подавления шума (параметр г в формуле (12) равен 0.8). Для лучшего восприятия фиг. 4 и 5 был применён фильтр повышения резкости. На фиг. 6 показано разностное изображение кадров фиг. 4 и 5.
Получение рентгенограмм осуществляют, например, с помощью устройства, показанного на фиг. 1. Оно содержит рентгеновскую трубку 1, которая испускает пучок рентгеновского излучения 2. Пучок рентгеновского излучения 2 поступает на детектор 3. Детектор 3 содержит сцинтилляционный экран (на чертеже не показан) и матрицу фоточувствительных элементов (на чертеже не показаны). Сцинтилляционный экран оптически связан с поверхностью активной области матрицы фоточувствительных элементов. Пучок рентгеновского излучения 2 попадает на детектор 3, сцинтилляционный экран преобразует его в видимый свет, который преобразуется сенсорами детектора в цифровую форму.
Предлагаемый в настоящем изобретении фильтр осуществляет подавление шума в несколько этапов. Далее рассмотрим каждый из этих этапов более подробно.
Этап 1. Покадровая оценка дисперсии шума, зависящего от интенсивности сигнала. На данном этапе вначале осуществляют оценку изображения с помощью линейной низкочастотной фильтрации текущего кадра [Непке1 е!. а1, РоЬик! апб Еак! Екбтабоп оГ 81дпа1-Берепбеп! Ыо18е ш Меб1са1 Х-Рау 1таде 8едиепсек, 8ргтдег, 2006]. Для обеспечения жёстких требований к скорости выполнения в приложениях реального времени целесообразно использовать простейшую линейную фильтрацию (например, биномиальным фильтром). На основе полученного сглаженного изображения строят изображение шума - разницу между исходным и фильтрованным изображениями.
Оценка изображения простейшими фильтрами неидеальна - границы оказываются пересглаженны
- 5 017302 ми. В результате при взятии разности между исходным изображением и сглаженным изображением получают изображение шума, содержащее помимо шумовых пикселей в гладких областях определённое количество пикселей, соответствующих резким изменениям (например, границы анатомических структур - далее нешумовые пиксели). Данные пиксели могут существенно исказить оцениваемое значение дисперсии и должны быть исключены из расчёта. Для этого применяются различные методики [Εοί с1 а1., Ртас!юа1 рощкошап-даикйап поще тобейпд апб ГбОпд Гог кшд1е-1таде га\\-ба1а. 1таде Ргосеккшд, ΙΕΕΕ Тгапкасбопк оп 17, 1737-1754 (Ос1оЬет 2008), 8а1теп е! а1., 81дпа1-берепбеп1 Ыо1ке С11агас1епха1юп Гог МаттодтарЫс 1тадек Оепоыпд. ΙΜΕΚΟ ТС4 8утрокшт (1МЕКОТС4 '08), Пгепхе. 1!а1у, 8ер1етЬег 2008], как правило, сводящиеся к пороговой обработке сглаженных производных исходного изображения, при этом величина порога определяется локальной оценкой отношения сигнал/шум. В областях изображения, содержащих большое количество деталей, такая оценка оказывается обычно неудовлетворительной. Поэтому в настоящем изобретении при отбрасывании границ предлагают более простой подход к выделению границ, не требующий вычисления производных и локальных оценок стандартного отклонения. Суть данного морфологического подхода к отбрасыванию нешумовых пикселей состоит в следующем:
1) изображение шума разбивают на две составляющие - бинарные изображения положительных и отрицательных изменений;
2) для выделения областей, соответствующих границам, полученные изображения подвергают морфологическим операциям эрозии с последующей дилатацией;
3) обработанные изображения объединяют для получения одного бинарного изображения - карты границ исходного изображения.
С целью более полного сохранения тонких структур морфологические операции эрозии и дилатации выполняют масками малого размера (окно 2x2).
При вычислении интервальных оценок дисперсии шума находят минимальную и максимальную интенсивности оценочного изображения (края диапазона интенсивностей), выбирают шаг разбиения, равный, например, 32 градациям серого цвета. Далее для каждого пикселя оценочного изображения находят интервал, которому принадлежит значение данного пикселя, и соответствующее значение пикселя изображения шума используют для вычисления оценки дисперсии шума в данном интервале (при этом исключают пиксели, соответствующие границам). При вычислении интервальной оценки дисперсии шума можно применить различные формулы, например обычную несмещённую оценку, или робастную оценку по формуле медианы абсолютных отклонений |Непке1 е! а1., КоЬик! апб Еай Εκί^таί^οи оГ 81дпа1Эерепбеп1 Ыо1ке т Мебка1 Х-Кау 1таде Зециепсек, 8ргтдег, 2006; Ео1 е! а1., Ргасбса1 рощкошап-даикйап поще тобейпд апб Гйбпд Гог йпд1е-1таде га^-ба!а, 1таде Ртосекйпд, ΙΕΕΕ Ттапкасбопк оп 17, 1737-1754, 2008]. На выходе данной процедуры получают табличную функцию, описывающую зависимость дисперсии шума от интенсивности сигнала.
Неточности, возникающие при построении карты границ, могут привести к грубым ошибкам при вычислении интервальных оценок дисперсии. Поэтому данные оценки уточняют, используя технику итеративного отбрасывания выбросов [Непке1 е!. а1., КоЬик! апб Еак! Εй^та!^οи оГ 8|дпа1-0ерепбеп1 ШЬе т Меб1са1 Х-Кау 1таде Зециепсек, 8ргшдет, 2006]: для каждого интервала итеративно исключают значения шумовых пикселей, превышающие по абсолютной величине порог, равный трём стандартным отклонениям шума, с последующим пересчётом оценки дисперсии шума в данном интервале.
После того как вычислены интервальные оценки дисперсии шума, осуществляют оценку зависимости дисперсии шума от интенсивности. При параметрической оценке в принципе можно построить тем или иным способом (например, методом наименьших квадратов, минимизации функции правдоподобия, направленной оптимизации) оценку параметров модели шума (1). Возможен также учёт нелинейностей сенсора [Ео1 е! а1., Ртасбса1 рохккошап-даикйап поще тобейпд апб Гйбпд Гог ктд1е-1таде та^-ба!а, 1таде Ртосеккшд, ΙΕΕΕ Ттапкасбопк оп 17, 1737-1754, 2008]. Однако, как отмечено в работе [Непке1 е! а1., КоЬик! апб Еак! Εκ!^таί^οи оГ 8|дпа1-Эерепбеп1 ШЬе ш Мебюа1 Х-Кау 1таде Зециепсек, 8ргтдег, 2006], построение параметрической модели, адекватно описывающей зависимость дисперсии шума от интенсивности сигнала, в силу ряда факторов может вызвать серьезные затруднения. К данным факторам можно отнести нелинейности в работе сенсора, нелинейные предобработки исходных снимков (например, логарифмирование). Поэтому более выгодным с точки зрения удобства реализации и универсальности применения представляется подход, при котором на основе интервальных оценок дисперсии строится непараметрическая оценка зависимости шума от интенсивности.
В настоящем изобретении используют непараметрический подход к построению искомой зависимости, при котором по полученным интервальным оценкам дисперсии шума создают интерполирующую табличную функцию.
Данную табличную функцию формируют на основе робастной локальной линейной аппроксимации интервальных оценок дисперсии. Использование робастных методов позволяет дополнительно снизить влияние выбросов (грубых ошибок в интервальных оценках дисперсии), в то время как локальность аппроксимации обеспечивает повторение сложного хода кривой, описывающей зависимость шума от интенсивности. Таким образом, полученная табличная функция каждой интенсивности исходного изобра
- 6 017302 жения ставит в соответствие оценку дисперсии шума. Точками входа в таблицу могут служить, например, интенсивности оценочного изображения.
На практике, в случае параметрической оценки шума, может быть использован подход, при котором применяется преобразование, стабилизирующее дисперсию шума исходного изображения [81атск с1. а1, 1таде Ртосеккшд апб Эа1а Лпа1у81к: ТНе МиШксак Лрртоасй. СатЬпФде Ишуегайу Ргекк, 1998; Яне, Цифровая обработка изображений, М.: Техносфера, 2007]. При этом проблема фильтрации шума, зависящего от сигнала, сводится к задаче подавления аддитивного не зависимого от сигнала шума заданной дисперсии. В настоящем изобретении осуществляют непараметрическую оценку шума, поэтому при построении карты шума предлагают использовать следующий подход: на основе оценочного изображения и интерполирующей таблицы построить карту шума - изображение, каждый пиксель которого оценивает среднеквадратическое отклонение шума в соответствующем пикселе исходного изображения. Карта шума даёт попиксельную оценку шума с достаточной для практического применения точностью. Возможен вариант исполнения заявляемого способа, в котором при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка (отбрасывание границ), используют разность между текущим и скомпенсированным по движению ссылочным кадром.
Этап 2. Оценка и компенсация движения и фликкер-эффекта. В настоящем изобретении оценку движения осуществляют на основе поиска соответствия блоков, с учётом возможной вариативности интенсивности последовательных кадров. Для этого до применения процедуры поиска векторов движения текущий и ссылочный кадры делят на их поля средних значений. Поля средних получают фильтрацией соответствующих изображений усредняющим по окрестности фильтром, апертура которого равна размеру используемого блока оценки движения, например 16x16 пикселей. При оценке и компенсации движения блоки оценки движения рассматривают с перекрытием, равным, например, четверти размера блока. При поиске соответствия блоков используют пирамидальный алгоритм оценки движения, при котором текущий кадр и ссылочный кадр, который должен быть скомпенсирован по движению по отношению к текущему кадру, раскладывают на несколько уровней в пирамиду изображений;
оценку векторов перемещения блоков осуществляют на низкочастотных пирамидах изображений;
при переходе с нижнего на верхний уровень пирамиды размеры блоков и их перекрытий масштабируют, а сам вектор движения данного блока уточняют поиском в области малого радиуса;
при построении скомпенсированного по движению кадра области перекрытия блоков усредняют.
Оценку движения от ссылочного кадра к текущему осуществляют на низкочастотных изображениях пирамиды. Используют не более трёх уровней разложения, как правило два, поскольку при более глубоком разложении при оценке движения мелких объектов (тонкие сосуды) можно допустить грубые ошибки в оценке вектора движения. При переходе с нижнего на верхний уровень пирамиды размеры блока и их перекрытия удваивают и осуществляют уточняющий поиск вектора движения по области небольшого размера, например 3x3 пикселя. При расчёте квадрата евклидовой метрики как меры подобия блоков используют технику раннего прекращения вычислений [^апд е1 а1., Еак! Л1дотййт8 йот 1Пе Екйтайои ой МоДои УесЮгк. 1ЕЕЕ Тгапк. 1таде Ртосеккшд, уо1. 8, по. 3, рр. 435-438, Матсй 1999]. Кроме того, применяют технику отсеивания неподвижных блоков из расчёта оценки вектора движения: если выборочное стандартное отклонение блока оценки движения не превосходит величины к <Лр), где к - значение порога (здесь для данного параметра используют значение от 2 до 3), - величина стандартного отклонения шума кадра в центральном пикселе блока оценки движения.
После оценки движения формируют скомпенсированный ссылочный кадр, при этом для уменьшения эффекта блочности перекрытия блоков усредняют. Усреднение перекрытий блоков осуществляют с использованием гладкого весового окна. Далее с целью компенсации фликкер-эффекта приводят локальные средние ссылочного кадра к локальным средним текущего кадра. Для этого делят скомпенсированный ссылочный кадр на его скомпенсированное по движению поле средних и умножают на поле средних значений текущего кадра согласно следующей формуле:
в которой ДчОО - скомпенсированный по движению ссылочный кадр;
Цр) . значение поля средних текущего кадра серии изображений в пикселе р;
Л-ι (/?) . скомпенсированное по движению поле средних значений ссылочного изображения;
Ζ-ι(ρ). результирующее, скомпенсированное по движению ссылочное изображение, локальная средняя интенсивность которого приведена к локальной средней интенсивности текущего изображения. Соответствующему преобразованию подвергают и карту дисперсии шума ссылочного кадра.
Этап 3. Рекурсивное усреднение с коррекцией артефактов компенсации движения
Грубые ошибки в компенсации движения, которые могут возникнуть из-за слишком малого радиуса области поиска, перекрытия объектов или резкой смены содержимого кадров, в настоящем изобретении уменьшают с помощью смешивания отфильтрованного кадра с усреднённой в пространстве изображения
- 7 017302 его версией. Веса смешивания вычисляют на основе пиксельного и структурного подобий текущего изображения скомпенсированному по движению изображению. В качестве меры пиксельного подобия используют билатеральное расстояние [Тотав1 с1. а1, Вйа!ега1 ННегшд Гог Сгау аиб Со1ог 1шадев, Ιη Ргос. 6111 Ιηΐ. СоиГ. Сошри1ег УЫоп. Ыете Бе1Ы, 1пб1а, 1998, рр. 839-846]. Для вычисления структурного подобия применяют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров [Виабев е! а1., Ыои1оса1 ипаде апб тоу1е бепоыпд. Ιηΐ. I. сотри!. У15юп, 76(2):123-139, 2008].
Пусть 1г(р), 1с(р) обозначают значения интенсивностей ссылочного и текущего кадров в пикселе р соответственно. Тогда пиксельное подобие кадров вычисляют согласно следующей пороговой функции:
где кр, Тр - параметры, определяющие характерную форму функции; а2(1с) - величина дисперсии шума текущего кадра в пикселе р. В свою очередь, структурное подобие между этими кадрами вычисляют по формуле
1+е(*,1ИЛ)-/Ч/Л2Р(/Д-/ ) ’ где кв, Тв - параметры, определяющие форму пороговой функции; Р(1с), Р(1в) - квадратные блоки пикселей (патчи) соответствующих кадров с центром в пикселе с номером р; параметры ||Р(1с)-Р(1г)||2 квадрат евклидова расстояния между блоками пикселей.
Для вычисления результирующей меры подобия кадров, объединяющей пиксельное и структурное сходство кадров, используют соотношение (5)
Параметр λ определяет преобладание структурного веса над весом, вычисленным по интенсивности соответствующих пикселей. Данная форма вычисления смешанного веса, при значениях параметра λ, близких к 1 (например, 0,9), позволяет при смешивании кадров уменьшить влияние одиночных импульсных значений интенсивности пикселей, одновременно восстанавливая правильное структурное содержание кадра (подобное содержанию исходного текущего неотфильтрованного кадра).
Временное усреднение осуществляют рекурсивно с помощью одномерной Калмановской фильтрации для скомпенсированного по движению ссылочного кадра (кадр-предиктор) и текущего кадра (наблюдаемый кадр) в соответствии со следующими формулами 1.
Оценка интенсивности по Калману:
(6) £(/>)= Ъ(с)+ *,(М(р)-/дМ (7)
Δ,(ρ)=Δμ(ρ)-^,(ρ)·Α%(ρ), (8)
2). Компенсация артефактов:
(7) (8)
3). Регулируемая сила подавления шума:
ра), скомпенсированная по движению и нормированная по интенсивности, данное значение берут из карты шума; А,(р) - величина дисперсии шума наблюдаемого кадра (текущего кадра);
А-ι (р) . интенсивность ссылочного кадра, скомпенсированного по движению и интенсивности (см. формулу (4)); 1|(р) - интенсивность наблюдаемого кадра;
Др) - рекурсивная оценка интенсивности наблюдаемого кадра;
- рекурсивная оценка дисперсии наблюдаемого кадра;
\У(р) - мера структурного и пиксельного подобия рекурсивно-усреднённого кадра наблюдаемому кадру (см. формулу (5));
1е |(р) - оценка интенсивности с помощью усреднения в пространстве наблюдаемого кадра биномиальным фильтром размера 3x3 пикселя; г - регулируемая сила подавления шума. В качестве меры ошиб ки наблюдаемого кадра используют рассчитанную карту дисперсии шума Д1(р), при этом карту дисперсии шумакадра-предиктора компенсируют по движению и интенсивности на основе найденных
на этапе оценки движения векторов и рассчитанных полей средних значений.
Лучший вариант осуществления изобретения
Для оценки дисперсии шума текущего кадра вначале получают оценочное изображение и изображение шума, для чего исходное изображение 1(х, у) фильтруют следующим низкочастотным линейным биномиальным фильтром размера 3x3:
Я1=[1 2 1]/4, Н~НТ-Н'.
При этом получают сглаженное изображение 1е(х, у) = 1*Н. Далее вычисляют изображение шума Щх, у) = 1(х, у) - 1е(х, у).
При отбрасывания границ из расчётов дисперсии на основе изображение шума Ые(х,у) формируют изображения положительных и отрицательных изменений
Для выделения крупных областей, соответствующих границам, данные бинарные изображения последовательно подвергают морфологическим операциям эрозии и дилатации, используя маски размера 2x2:
Затем эти изображения объединяют для получения карты границ исходного изображения Е(х, у) = В-е(х, у)пВ+е(х, у).
Для вычисления интервальных оценок дисперсии шума находят минимальную 1т1п и максимальную 1тах интенсивности изображения 1е(х, у), выбирают шаг йМ и разбивают диапазон интенсивностей на интервалы М1 с шагом ЬМ (ЬМ = 32). Для каждого пикселя изображения 1е(х, у) находят интервал, которому принадлежит значение этого пикселя, а соответствующее значение пикселя изображения Ые(х, у) используют для вычисления оценки дисперсии шума σ2(ί) в данном интервале М1, при этом исключают значения пикселей, для которых значение карты границ Е(х, у) = 1. При вычислении интервальной оценки дисперсии шума применяют формулу несмещённой оценки дисперсии
где №е(]) - значение пикселя изображения шума Ые(х, у) из интервала М1, п1 - общее количество накопленных значений пикселей изображения шума в интервале - среднее значений шумовых пикселей в интервале М1.
Полученные интервальные оценки дисперсии σ2(ί) уточняют таким образом, что для каждого интервала М1 итеративно исключают значения шумовых пикселей, превышающие по абсолютной величине порог, равный трём стандартным отклонениям шума, с последующим пересчётом оценки дисперсии шума в данном интервале:
где ЫДк+Ц - набор значений шумовых пикселей в интервале М1 на итерации к. В дальнейших расчётах используются лишь те интервалы, в которых накоплено достаточное количество значений шумовых пикселей (например, не меньше 500). Кроме того, из рассмотрения исключаются те интервалы,
ДР среднее значение е в которых существенно отличается от нуля (отклоняется от нуля более чем на половину шага сетки интервалов йМ), поскольку в таких интервалах с высокой вероятностью доминируют остаточные значения нешумовых пикселей.
Для оценки зависимости дисперсии шума от интенсивности по полученным интервальным оценкам дисперсии шума (отмечены кружками на фиг. 2) создают интерполирующую табличную функцию. Данную табличную функцию формируют на основе робастной локальной линейной аппроксимации интервальных оценок дисперсии. Для этого на сетке интервалов М1 выбирают шаг 11| и радиус г аппроксимации (которые можно сделать зависимыми от количества точек п1 в интервалах М1). Значения полученной на предыдущем этапе табличной функции дисперсии шума от интенсивности сигнала аппроксимируют по следующему правилу:
- 9 017302
где к - номер узла аппроксимации; тк - центр интервала М1(к-й1) ;
параметры а, Ь вычисляют на основе минимума суммы абсолютных отклонений |Рге55 с1 а1.. Νυтепса1 Рсс1рс5 ίη С: ТНе Ай οί ΞοίοηΙίΠο С'отриПпд. 8есопб Εάίΐίοη. Νο\ν Уогк: СатЬпбде итуегзйу Ргс55. 1999] к/г/+г! σ\])-α·ιη -ό|->ηιίη ауЬ '
Полученную в результате такого процесса таблицу значений σ2(*Λ)} интерполируют на всём диапазоне интенсивностей [1тт.1тах]. т.е. интерполирующую таблицу искомой зависимости дисперсии шума от интенсивности сигнала (сплошная линия фиг. 2).
При построении карты шума на основе оценочного изображения и интерполирующей таблицы строят карту шума - изображение. каждый пиксель которого оценивает дисперсию шума в соответствующем пикселе исходного изображения (фиг. 3).
Оценку и компенсацию движения осуществляют на основе поиска соответствия блоков. с учётом возможной вариативности интенсивности последовательных кадров. Для этого перед оценкой движения текущий и ссылочный кадры делят на их поля средних значений. Данные поля средних получают фильтрацией соответствующих изображений усредняющим по окрестности фильтром. размер которого равен размеру используемого блока оценки движения. например 16x16 пикселей. При оценке движения блоки рассматривают с перекрытием. равным. например. 4 пикселям. При поиске соответствия блоков используют пирамидальную схему:
текущий кадр и ссылочный кадр. который должен быть скомпенсирован по движению по отношению к текущему кадру. раскладывают на несколько уровней в пирамиду изображений;
оценку векторов перемещения блоков осуществляют на низкочастотных пирамидах изображений;
при переходе с нижнего на верхний уровень пирамиды размеры блоков и их перекрытий масштабируют. а сам вектор движения данного блока уточняют поиском в области малого радиуса;
при построении скомпенсированного по движению кадра области перекрытия блоков усредняют.
Оценку движения от ссылочного кадра к текущему осуществляют на низкочастотных изображениях пирамиды. Используют не более трёх уровней разложения. как правило. два уровня. поскольку при более глубоком разложении. в особенности при оценке движения мелких объектов (сосуды). на грубых уровнях разрешения тонкие структуры пропадут. что приведёт к грубым ошибкам в оценке векторов движения. При переходе с нижнего на верхний уровень пирамиды размеры блока и их перекрытия удваивают и осуществляют уточняющий поиск вектора движения по области небольшого размера. например 3x3 пикселя. При расчёте квадрата евклидовой метрики как меры подобия блоков используют технику раннего прекращения вычислений |\Уа1щ е1. а1.. Еай А1догййт8 1ог 1йе ΕδΙίιηηΙίοη οί Μοίίοη Уес1ог5. ΙΕΕΕ Тгащ. 1таде Ргосеккшд. νοί. 8. ηο. 3. рр. 435-438. Магсй 1999]. Кроме того. применяют технику отсеивания неподвижных блоков из расчёта оценки вектора движения: если выборочное стандартное отклонение блока оценки движения не превосходит величины ' σ(ρ), где к - значение порога (здесь для данного параметра используют значение от 2 до 3), - величина стандартного отклонения шума кадра в центральном пикселе блока оценки движения. После оценки движения формируют скомпенсированный ссылочный кадр. при этом для уменьшения эффекта блочности перекрытия блоков усредняют. Аналогичным образом компенсируют карту шума ссылочного кадра. Усреднение перекрытий блоков осуществляют с использованием гладкого весового окна. например косинусного.
Для компенсации фликкер-эффекта приводят локальные средние ссылочного кадра к локальным средним текущего кадра. Для этого делят скомпенсированный ссылочный кадр на его скомпенсированное по движению поле средних и умножают на поле средних значений текущего кадра согласно формуле (4). Карту шума ссылочного кадра подвергают аналогичной нормировке.
Рекурсивная фильтрация текущего кадра
Временное усреднение скомпенсированных по движению и интенсивности кадров осуществляют с помощью рекурсивной фильтрации. используя одномерный (попиксельный) фильтр Калмана. Расчёт ведут по формулам (6)-(12). Для подавления артефактов временного усреднения. обусловленных ошибками в компенсации движения. используют смешивание фильтрованного нормализованного кадра с его исходной версией (формулы (5). (9). (10)). Веса смешивания (5) находят на основе вычисления яркостного и структурного подобия текущего кадра скомпенсированному по движению и интенсивности ссылочному кадру. при этом в качестве меры структурного подобия используют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров. Расчёт структурного подобия вычисляют скользящим блоком. размер которого выбирают равным. например. 11x11 пикселей. Расчёт пиксельного подобия вычисляют как расстояние в пространстве интенсивностей между текущим кадром и ссылочным. скомпенсированным по движению кадром. Силу подавления шума г в формулах (11). (12) на реальных
- 10 017302 сериях берут равной, например, г = 0,8, что соответствует пятикратному подавлению уровня шума. Фиг. 4-6 демонстрируют результат применения фильтра к реальной серии (ангиографическое исследование): фиг. 4 показывает исходный кадр серии, фиг. 5 - фильтрованный кадр, фиг. 6 - разность между кадрами фиг. 4-5; для улучшения восприятия исходный и отфильтрованный кадры были подвергнуты преобразованию повышения резкости.

Claims (2)

  1. ФОРМУЛА ИЗОБРЕТЕНИЯ
    1. Способ подавления шума серий цифровых рентгенограмм, включающий получение серии цифровых рентгенограмм, покадровую оценку дисперсии шума, зависящего от сигнала, оценку и компенсацию движения, оценку и компенсацию фликкер-эффекта, рекурсивное усреднение кадров, отличающийся тем, что при покадровой оценке дисперсии шума, зависящего от интенсивности сигнала, осуществляют морфологическое удаление значений пикселей снимка шума, соответствующих границам в исходном снимке; получают табличную зависимость шума от интенсивности сигнала с помощью робастной кусочно-линейной аппроксимации интервальных оценок дисперсии шума; вычисляют карту шума в виде попиксельной оценки дисперсии шума исходного цифрового изображения; при оценке и компенсации движения и фликкер-эффекта используют пирамидальную схему поиска соответствия блоков; блоки оценки движения рассматривают с перекрытием, при компенсации движения перекрытия блоков усредняют с использованием гладкого весового окна; учёт фликкер-эффекта осуществляют делением ссылочного и текущего кадров на соответствующие поля средних значений, которые получают линейной НЧфильтрацией данных кадров; для компенсации фликкер-эффекта приводят локальную среднюю яркость ссылочного кадра к локальной средней яркости текущего, для чего делят ссылочный кадр на его скомпенсированное по движению поле средних значений и умножают на поле средних значений текущего кадра; при рекурсивном усреднении с коррекцией артефактов компенсации движения осуществляют смешивание текущего отфильтрованного кадра с его исходной версией на основе вычисления пиксельного и структурного подобия данных кадров.
  2. 2. Способ по п.1, отличающийся тем, что при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка, используют разность между текущим и скомпенсированным по движению ссылочным кадром.
EA201101502A 2011-10-07 2011-10-07 Способ подавления шума серий цифровых рентгенограмм EA017302B1 (ru)

Priority Applications (6)

Application Number Priority Date Filing Date Title
EA201101502A EA017302B1 (ru) 2011-10-07 2011-10-07 Способ подавления шума серий цифровых рентгенограмм
EP12186390.6A EP2579207A3 (en) 2011-10-07 2012-09-27 Method of noise reduction in digital x-ray frames series
CN2012103719525A CN103177423A (zh) 2011-10-07 2012-09-28 对数字x射线帧系列降噪的方法
JP2012216792A JP2013127773A (ja) 2011-10-07 2012-09-28 デジタルｘ線フレームシリーズにおけるノイズ低減方法
US13/645,518 US20130089247A1 (en) 2011-10-07 2012-10-05 Method of Noise Reduction in Digital X-Ray Frames Series
KR1020120110554A KR20130038794A (ko) 2011-10-07 2012-10-05 디지털 엑스레이 프레임들 시리즈에서의 잡음 감소의 방법

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
EA201101502A EA017302B1 (ru) 2011-10-07 2011-10-07 Способ подавления шума серий цифровых рентгенограмм

Publications (2)

Publication Number Publication Date
EA201101502A1 EA201101502A1 (ru) 2012-10-30
EA017302B1 true EA017302B1 (ru) 2012-11-30

Family

ID=47136908

Family Applications (1)

Application Number Title Priority Date Filing Date
EA201101502A EA017302B1 (ru) 2011-10-07 2011-10-07 Способ подавления шума серий цифровых рентгенограмм

Country Status (6)

Country Link
US (1) US20130089247A1 (ru)
EP (1) EP2579207A3 (ru)
JP (1) JP2013127773A (ru)
KR (1) KR20130038794A (ru)
CN (1) CN103177423A (ru)
EA (1) EA017302B1 (ru)

Families Citing this family (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6071444B2 (ja) * 2012-11-07 2017-02-01 キヤノン株式会社 画像処理装置及びその作動方法、プログラム
CA2901847A1 (en) 2013-03-15 2014-09-18 Yale University Techniques for processing imaging data having sensor-dependent noise
KR102104403B1 (ko) * 2013-05-28 2020-04-28 한화테크윈 주식회사 단일영상 내의 안개 제거 방법 및 장치
CN104240184B (zh) * 2013-06-08 2017-09-26 通用电气公司 噪声标准差的估算方法和系统
BR112016007072A2 (pt) * 2013-10-01 2017-08-01 Agfa Healthcare método para a redução de ruido em uma sequência de imagens
JP6254813B2 (ja) * 2013-10-10 2017-12-27 キヤノン株式会社 Ｘ線画像処理装置および方法、並びにｘ線撮像装置、プログラム
EP3055836B1 (en) * 2013-10-11 2019-03-20 Mauna Kea Technologies Method for characterizing images acquired through a video medical device
CN104182948B (zh) * 2013-12-23 2015-07-22 上海联影医疗科技有限公司 一种相关性噪声的估计方法
US10043243B2 (en) * 2016-01-22 2018-08-07 Siemens Healthcare Gmbh Deep unfolding algorithm for efficient image denoising under varying noise conditions
WO2017132600A1 (en) * 2016-01-29 2017-08-03 Intuitive Surgical Operations, Inc. Light level adaptive filter and method
CN106355559B (zh) * 2016-08-29 2019-05-03 厦门美图之家科技有限公司 一种图像序列的去噪方法及装置
EP3438922B1 (en) * 2017-07-27 2019-11-20 Siemens Healthcare GmbH Method for processing at least one x-ray image
GB2567668B (en) * 2017-10-20 2022-03-02 Vivo Mobile Communication Co Ltd Deflickering of a series of images
TWI680675B (zh) * 2017-12-29 2019-12-21 聯發科技股份有限公司 影像處理裝置與相關的影像處理方法
US10891720B2 (en) 2018-04-04 2021-01-12 AlgoMedica, Inc. Cross directional bilateral filter for CT radiation dose reduction
US11080898B2 (en) 2018-04-06 2021-08-03 AlgoMedica, Inc. Adaptive processing of medical images to reduce noise magnitude
FR3083415B1 (fr) * 2018-06-29 2020-09-04 Electricite De France Traitement d'un bruit impulsionnel dans une sequence video
CN110730280B (zh) * 2018-07-17 2021-08-31 瑞昱半导体股份有限公司 噪声等化方法与噪声去除方法
CN110008958B (zh) * 2018-08-01 2021-04-16 金华市灵龙电器有限公司 直流无刷电机控制机构
IT201900001225A1 (it) * 2019-01-28 2020-07-28 General Medical Merate S P A Metodo predittivo per controllare un'apparecchiatura radiologica e apparecchiatura radiologica che lo implementa
KR102304124B1 (ko) * 2019-02-21 2021-09-24 한국전자통신연구원 학습기반 3d 모델 생성 장치 및 방법
JP7403994B2 (ja) * 2019-08-21 2023-12-25 富士フイルムヘルスケア株式会社 医用画像処理装置および医用画像処理方法
DE102019122667A1 (de) * 2019-08-22 2021-02-25 Schölly Fiberoptic GmbH Verfahren zur Unterdrückung von Bildrauschen in einem Videobildstrom, sowie zugehöriges medizinisches Bildaufnahmesystem und Computerprogrammprodukt
CN111242831B (zh) * 2020-01-20 2022-11-08 暨南大学 一种基于Zernike矩抗几何攻击的可逆鲁棒水印方法
CN113709324A (zh) * 2020-05-21 2021-11-26 武汉Tcl集团工业研究院有限公司 一种视频降噪方法、视频降噪装置及视频降噪终端
CN113160072B (zh) * 2021-03-19 2023-04-07 聚融医疗科技(杭州)有限公司 一种基于图像金字塔的鲁棒自适应帧相关方法及系统
WO2022265875A1 (en) * 2021-06-18 2022-12-22 Subtle Medical, Inc. Systems and methods for real-time video denoising
CN114114265A (zh) * 2021-12-13 2022-03-01 中国地质大学（北京） 星载全极化sar电离层闪烁误差补偿中的平滑窗设计方法
CN114638757B (zh) * 2022-03-17 2024-08-02 苏州工业园区智在天下科技有限公司 降噪方法、装置、终端及计算机可读存储介质
CN115661006B (zh) * 2022-12-29 2023-03-10 湖南国天电子科技有限公司 一种海底地貌图像去噪方法
CN116249018B (zh) * 2023-05-11 2023-09-08 深圳比特微电子科技有限公司 图像的动态范围压缩方法、装置、电子设备和存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1392047A2 (en) * 2002-08-22 2004-02-25 Samsung Electronics Co., Ltd. Digital document processing for image enhancement
RU2322694C2 (ru) * 2006-03-09 2008-04-20 Общество с ограниченной ответственностью "Комэксп" Способ обработки изображений
KR20080044737A (ko) * 2006-11-16 2008-05-21 주식회사 메디슨 초음파 영상 처리 방법
RU2342701C1 (ru) * 2007-08-15 2008-12-27 Российская Федерация, от имени которой выступает Министерство обороны Российской Федерации Способ комплексирования цифровых многоспектральных полутоновых изображений
US20090195535A1 (en) * 2008-02-05 2009-08-06 Docomo Communications Laboratories Usa, Inc. Methods for fast and memory efficient implementation of transforms

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5923775A (en) * 1996-04-04 1999-07-13 Eastman Kodak Company Apparatus and method for signal dependent noise estimation and reduction in digital images
DE60141961D1 (de) * 2001-09-10 2010-06-10 Texas Instruments Inc Verfahren und Vorrichtung zur Bewegungsvektorenabschätzung
US7187794B2 (en) * 2001-10-18 2007-03-06 Research Foundation Of State University Of New York Noise treatment of low-dose computed tomography projections and images
US8107535B2 (en) * 2003-06-10 2012-01-31 Rensselaer Polytechnic Institute (Rpi) Method and apparatus for scalable motion vector coding
US20070092007A1 (en) * 2005-10-24 2007-04-26 Mediatek Inc. Methods and systems for video data processing employing frame/field region predictions in motion estimation
US8965078B2 (en) * 2009-02-20 2015-02-24 Mayo Foundation For Medical Education And Research Projection-space denoising with bilateral filtering in computed tomography
ATE547775T1 (de) * 2009-08-21 2012-03-15 Ericsson Telefon Ab L M Verfahren und vorrichtung zur schätzung von interframe-bewegungsfeldern
US20110134315A1 (en) * 2009-12-08 2011-06-09 Avi Levy Bi-Directional, Local and Global Motion Estimation Based Frame Rate Conversion
CN102014240B (zh) * 2010-12-01 2013-07-31 深圳市蓝韵实业有限公司 一种实时医学视频图像去噪方法
US8538114B2 (en) * 2011-06-06 2013-09-17 Kabushiki Kaisha Toshiba Method and system utilizing parameter-less filter for substantially reducing streak and or noise in computer tomography (CT) images

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1392047A2 (en) * 2002-08-22 2004-02-25 Samsung Electronics Co., Ltd. Digital document processing for image enhancement
RU2322694C2 (ru) * 2006-03-09 2008-04-20 Общество с ограниченной ответственностью "Комэксп" Способ обработки изображений
KR20080044737A (ko) * 2006-11-16 2008-05-21 주식회사 메디슨 초음파 영상 처리 방법
RU2342701C1 (ru) * 2007-08-15 2008-12-27 Российская Федерация, от имени которой выступает Министерство обороны Российской Федерации Способ комплексирования цифровых многоспектральных полутоновых изображений
US20090195535A1 (en) * 2008-02-05 2009-08-06 Docomo Communications Laboratories Usa, Inc. Methods for fast and memory efficient implementation of transforms

Also Published As

Publication number Publication date
EP2579207A2 (en) 2013-04-10
EA201101502A1 (ru) 2012-10-30
US20130089247A1 (en) 2013-04-11
KR20130038794A (ko) 2013-04-18
JP2013127773A (ja) 2013-06-27
EP2579207A3 (en) 2013-07-10
CN103177423A (zh) 2013-06-26

Similar Documents

Publication Publication Date Title
EA017302B1 (ru) Способ подавления шума серий цифровых рентгенограмм
US10292672B2 (en) Radiographic image processing device, method, and recording medium
US10282823B2 (en) Simulating dose increase by noise model based multi scale noise reduction
WO2013125984A2 (ru) Способ подавления шума цифровых рентгенограмм
US7689055B2 (en) Method and apparatus for enhancing image acquired by radiographic system
US7831097B2 (en) System and method for image reconstruction
US10194881B2 (en) Radiographic image processing device, method, and recording medium
RU2441281C1 (ru) Способ оценки шума цифровых рентгенограмм
US9996910B2 (en) Radiographic image processing device, method, and recording medium
US20150081262A1 (en) Method and system for statistical modeling of data using a quadratic likelihood functional
CN111899188B (zh) 一种神经网络学习的锥束ct噪声估计与抑制方法
Choi et al. Speckle noise reduction in ultrasound images using SRAD and guided filter
JP2015024155A (ja) 画像処理装置
US20080285881A1 (en) Adaptive Image De-Noising by Pixels Relation Maximization
KR20130083205A (ko) 양전자 방출 단층 촬영에서의 영상 보정 방법 및 장치
US20070140582A1 (en) Systems and Methods For Reducing Noise In Image Sequences
CN102890817A (zh) 医学图像的处理方法及系统
Ullherr et al. Using measurements of the spatial SNR to optimize phase contrast X-ray imaging
KR20100121273A (ko) 엑스선 영상 처리 장치 및 방법
CN106875342B (zh) 一种计算机断层图像处理方法和装置
TW201417769A (zh) 一種影像品質的改善處理方法及其造影系統
Auvray et al. Multiresolution parametric estimation of transparent motions and denoising of fluoroscopic images
Vashistha et al. ParaPET: non-invasive deep learning method for direct parametric PET reconstruction using histoimages
Malpica et al. Endocardial tracking in contrast echocardiography using optical flow
Kim Restoration of Chest X-ray by Kalman Filter

Legal Events

Date Code Title Description
MM4A Lapse of a eurasian patent due to non-payment of renewal fees within the time limit in the following designated state(s)

Designated state(s): MD

MM4A Lapse of a eurasian patent due to non-payment of renewal fees within the time limit in the following designated state(s)

Designated state(s): AM AZ BY KZ KG TJ TM