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

СПОСОБ ПОДАВЛЕНИЯ ШУМА СЕРИЙ ЦИФРОВЫХ РЕНТГЕНОГРАММ Особенностью предлагаемого способа подавления шума серий цифровых рентгенограмм является то, что при покадровой оценке дисперсии шума, зависящего от интенсивности сигнала,осуществляют морфологическое удаление значений пикселей снимка шума, соответствующих границам в исходном снимке. Получают табличную зависимость шума от интенсивности сигнала с помощью робастной кусочно-линейной аппроксимации интервальных оценок дисперсии шума. Вычисляют карту шума в виде попиксельной оценки дисперсии шума исходного цифрового изображения. При оценке и компенсации движения и фликкер-эффекта используют пирамидальную схему поиска соответствия блоков. Блоки оценки движения рассматривают с перекрытием, при компенсации движения перекрытия блоков усредняют с использованием гладкого весового окна. Учт фликкер-эффекта осуществляют делением ссылочного и текущего кадров на соответствующие поля средних значений, которые получают линейной НЧ-фильтрацией данных кадров. Для компенсации фликкер-эффекта приводят локальную среднюю яркость ссылочного кадра к локальной средней яркости текущего, для чего делят ссылочный кадр на его скомпенсированное по движению поле средних значений и умножают на поле средних значений текущего кадра. При рекурсивном усреднении с коррекцией артефактов компенсации движения осуществляют смешивание текущего отфильтрованного кадра с его исходной версией на основе вычисления пиксельного и структурного подобия данных кадров. Технический результат заявляемого способа заключается в повышении качества получаемого изображения за счет подавления шумов. Подавление шумов осуществляется в 3 последовательных этапа: покадровая оценка дисперсии шума, зависящего от сигнала; оценка и компенсация движения и фликкерэффекта; рекурсивное усреднение кадров с коррекцией артефактов компенсации движения. Меркурьев Сергей Васильевич (RU) Пронин В.О. (RU)(71)(73) Заявитель и патентовладелец: ЗАКРЫТОЕ АКЦИОНЕРНОЕ ОБЩЕСТВО "ИМПУЛЬС" (RU) 017302 Область техники Изобретение относится к области обработки цифровых изображений и может быть использовано для решения задач обработки цифровых изображений, полученных с помощью высокоэнергетического излучения, в том числе рентгеновского. Более конкретно, данное изобретение предназначено для подавления шума серий цифровых рентгенограмм. Предшествующий уровень техники В настоящее время в клинической практике применяются различные методы обработки серий цифровых рентгенограмм (медицинских изображений): повышение резкости, выделение анатомических структур, субтракция и др. К важнейшим методам улучшения качества серий медицинских изображений относятся методы подавления шумов. Так, в ангиографии в реальном времени осуществляют управление медицинским инструментарием и исследуют состояние сосудистой системы с помощью введения через катетер контрастного вещества. С целью уменьшения лучевой нагрузки на пациента и медперсонал применяют малые дозы облучения, что приводит к получению медицинских изображений с низким значением соотношения сигнал/шум [Chan et. al., Image Sequence Filtering in Quantum-Limited Noise with Applications to Low-Dose Fluoroscopy, IEEE Trans. Medical Imaging, vol. 12, issue 3, 610-621, 1993]. Поэтому в цифровой радиологии представляется актуальной задача разработки методов фильтрации шумов, способных сохранить диагностически значимую информацию, работать в реальном времени при высоких значениях разрешения, частоты смены кадров и существенном уровне шума, зависимого от интенсивности сигнала. Можно выделить несколько проблем, с которыми сталкивается исследователь при разработке алгоритма фильтрации серий цифровых медицинских изображений. Первая проблема связана с оценкой шума изображений, уровень которого существенно зависит от интенсивности сигнала. Вторая проблема обусловлена необходимостью компенсации возможного резкого изменения интенсивности изображений,возникающего при использовании автоматической системы подстройки интенсивности или вследствие нестабильности рентгеновской трубки. Третья проблема связана с учтом изменений интенсивности последовательности кадров, вызванных движением: перемещением стола или системы трубка/детектор,смещением пациента или физиологическими процессами в его организме, током контрастного вещества в сосудах. Наконец, четвртая проблема состоит в необходимости соблюдения баланса между качеством работы фильтра и требованием выполнения фильтрации в реальном времени. Данная проблема включает также выбор метода фильтрации: пространственно-временное усреднение в объме стека изображений,временная фильтрация в пространстве преобразования (пирамидальное преобразование, вейвлетпреобразование и т.п.), комбинированные методы усреднения и т.п. Опишем каждую из этих проблем более подробно. Рассмотрим проблему оценки шума, зависящего от интенсивности сигнала. Практически каждый качественный метод фильтрации шума, как одиночных цифровых изображений, так и их серий, использует информацию о спектральных свойствах шума. Весьма широкий класс методов шумоподавления использует в качестве параметра величину дисперсии шума. В публикациях разных авторов продемонстрировано, что чем точнее известно значение дисперсии шума, тем качественнее результат работы алгоритма шумоподавления [Bosco et. al, Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi, Practical Denoising of Clipped or Overexposed Noisy Images, Proc. 16th EuropeanSignal Process. Conf., EUSIPCO, 2008]. Рассмотрим природу шума изображений в цифровой рентгенографии. Детектор измеряет интенсивность ослабленного (прошедшего через объект исследования) рентгеновского излучения. Каждая ячейка детектора накапливает за время экспонирования в среднем N электронов путм поглощения фотонов. Количество N накопленных электронов можно моделировать распределнной по закону Пуассона случайной величиной Случайные флуктуации числа поглощнных фотонов называют шумом фотонов. В современных детекторах основным источником шума является шум фотонов. К дополнительным источникам шума относят шумы детекторной системы: шум чтения, тепловой шум, шумы усилителей,шумы квантования и другие [Gino M., Noise, Noise, Noise]. Общий эффект данных источников шума моделируют распределнной по Гауссу случайной величиной [Bosco et. al., Noise Reduction for Cfa Imagemodeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]. Согласно широко используемой модели при линейных электронных схемах дисперсия шума фотонов и дополнительных источников шумов линейно зависит от величины полезного сигнала [Яне, Цифровая обработка изображений, М.: Техносфера, 2007] где I(p) - уровень интенсивности сигнала в пикселе p = (x, y) изображения. Таким образом, шум цифровых изображений линейно зависит от интенсивности сигнала. Проблеме-1 017302 оценки шума, зависящего от сигнала, посвящено достаточно много публикаций [Bosco et. al., Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi et. al., PracticalSequences, Springer, 2006] описан непараметрический метод покадровой оценки шума серий цифровых медицинских изображений, при этом особый акцент делают на разработке алгоритма, способного работать в реальном времени. В публикации [Foi et. al., Practical poissonian-gaussian noise modeling and fittingfor single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008] рассматривают двухпараметрический подход к оценке шума цифровых снимков. Шум исходного цифрового изображения (полученного непосредственно с детектора и не прошедшего через нелинейные преобразования, такие как гамма-коррекция и т.п.) моделируют аддитивной по отношению к сигналу случайной величиной: где In(p) - уровень интенсивности сигнала в пикселе p наблюдаемого зашумлнного изображения,2(I(p - зависимость дисперсии шума от интенсивности сигнала вида (1), n(p)N(0,1) - стандартная нормальная случайная величина. Предлагают способ построения модельных кривых дисперсии шума,учитывающий нелинейности в работе сенсора, приводящие к недоэкспонированию и переэкспонированию, т.е. к нарушению линейного закона (1) на краях динамического диапазона. Вторая проблема в разработке динамического фильтра заключается в необходимости учта и компенсации изменения интенсивности фильтруемой серии изображений, т.н. мерцаний. Данная проблема подобна задаче компенсации фликкер-эффекта, возникающей как необходимый первый этап в цепочке различных методов обработки кадров видеопоследовательности [Delon, Movie and Video Scale-TimeMotion Vectors Compensation, In Third International Conference on Image and Graphics, pp. 552-555, 2004]. Для серии цифровых рентгенограмм мерцания могут быть вызваны использованием автоматической системы подстройки интенсивности, либо появится вследствие нестабильной работы рентгеновской трубки. Отсутствие учта и компенсации мерцаний при временном усреднении серии изображений отрицательно сказывается на качестве работы фильтра. Так, временное усреднение стационарных по движению кадров без компенсации фликкер-эффекта приводит к существенному искажению оцениваемых значений текущего кадра. При наличии движения и его компенсации на основе методик типа поиска соответствия блоков пикселей (далее блоков), значительное изменение интенсивности изображений в серии приводит к грубым ошибкам в оценке движения и, как следствие, к плохой его компенсации. Можно выделить две главные стороны в оценке и компенсации мерцаний. С одной стороны, в общем случае фликкер-эффект имеет нелокальную природу по пространству отдельно взятого кадра серии изображений. С другой стороны, он проявляется только при рассмотрении динамики процесса, т.е. нескольких подряд идущих кадров. Эти особенности составляют главную трудность в компенсации фликкер-эффекта, состоящую в определении того, что локальные изменения в контрасте/интенсивности кадра являются следствием движения или обусловлены "прыгнувшей" интенсивностью излучения рентгеновской трубки. Естественный способ моделирования фликкер-эффекта, широко описанный в литературеImage and Graphics, pp. 552-555, 2004], основан на следующей формуле: где I0(p) - интенсивность исходного изображения в пикселе p; A(p), B(p) - гладкие, медленно изменяющиеся функции, определнные в области задания изображения; I(p) - интенсивность наблюдаемого изображения. Данная модель предполагает линейную зависимость интенсивности исходного изображения от интенсивности наблюдаемого изображения, а также учитывает отмеченную выше пространственную неоднородность фликкер-эффекта. Значительное количество подходов к оценке и компенсации мерцаний основано на разбиении изображения на сетку перекрывающихся блоков и поблочную оценку параметров функций A(p), B(p) с возможной отбраковкой выбросов, при этом для более аккуратной оценки параметры модели находят совместно с оценкой движения.-2 017302 Третья проблема, возникающая при разработке качественного метода фильтрации шумов серий изображений, заключается в необходимости учта изменений в последовательности кадров, обусловленных движением. Движение является обязательным атрибутом любой диагностически значимой серии цифровых медицинских изображений. Усреднение нестационарных по движению кадров приводит к появлению артефактов размытия на результирующем изображении. На практике широко применяются два класса подходов к решению проблемы учта движения при пространственно-временном усреднении кадров видеопоследовательности. Подходы первого типа основаны на детектировании регионов изображения, соответствующих движению и уменьшению силы временного подавления шума в этих регионах(2005) 253-274]. Подходы второго типа оценивают траектории перемещения пикселей между последовательными изображениями, для чего осуществляют учт и компенсацию движения [Brailean et al., NoiseReduction Filters for Dynamic Image Sequences: A review. Proc. IEEE, vol. 83, pp. 1272-1292, 1995]. Стационарные по движению изображения можно явно сформировать после оценки движения в результате применения процедуры компенсация движения. Простейшие из алгоритмов первого типа достаточно легко реализовать для выполнения в реальном времени. Для таких методов детектирование движения сводится, например, к пороговой обработке разности между текущим усредняемым кадром и ссылочным кадром. Если для некоторого пикселя текущего кадра данная разность выше определнного порогового значения, то временное усреднение заменяют,например, усреднением в пространстве текущего кадра. Величину порога естественно привязать к уровню шума, полученному на этапе оценки шума. Следует отметить, что подобные методы учта движения страдают определнными недостатками: появляется хорошо заметный импульсный шум, возникают негладкости на границах движущихся объектов. В свою очередь, алгоритмы второго типа, в которых производится оценка с последующей компенсацией движения (например, на основе поиска соответствия блоков в изображениях), считаются более качественными [Brailean et al., Noise Reduction Filters for Dynamic Image Sequences: A review. Proc. IEEE, vol. 83, pp. 1272-1292, 1995]. Однако есть определнные трудности и при их использовании: не все виды движения можно эффективно скомпенсировать [Яне,Цифровая обработка изображений, М.: Техносфера, 2007], кроме того, в силу относительно высокой вычислительной требовательности данных методов существуют определнные сложности при их реализации в реальном времени. Четвртая проблема состоит в необходимости обеспечения баланса между качеством работы фильтра и требованием выполнения в реальном времени. Необходимость выполнения обработки в реальном времени серии медицинских изображений при высокой частоте и большом размере кадров (например, 30 кадров в секунду при разрешении 512512 пикселей, 15 кадров в секунду при разрешении 10241024 пикселей в кадре) на современных рентгенодиагностических комплексах вступает в определнное противоречие с желанием обеспечить высокое качество результата фильтрации. Такие операции,как оценка шума или оценка движения, сами по себе являются достаточно требовательными к вычислительным ресурсам даже с учтом возможностей современных компьютеров. Раскрытие изобретения Целью заявляемого изобретения является улучшение качества серий цифровых медицинских изображений (цифровых рентгенограмм). Технический результат заявляемого способа заключается в повышении качества получаемого изображения за счет подавления шумов. Технический результат в заявляемом способе подавления шума серий цифровых рентгенограмм,включающем получение серии цифровых рентгенограмм, покадровую оценку дисперсии шума, зависящего от сигнала, оценку и компенсацию движения, оценку и компенсацию фликкер-эффекта, рекурсивное усреднение кадров достигается тем, что при оценке дисперсии шума осуществляют морфологическое удаление значений пикселей снимка шума, соответствующих резким изменениям в исходном снимке; получают табличную зависимость шума от интенсивности сигнала с помощью робастной кусочнолинейной аппроксимации интервальных оценок дисперсии шума; вычисляют карту шума в виде попиксельной оценки дисперсии шума исходного цифрового изображения; при оценке и компенсации движения используют пирамидальную схему поиска соответствия блоков, блоки оценки движения рассматривают с перекрытием, при компенсации движения перекрытия блоков пикселей усредняют; учт фликкерэффекта осуществляют делением ссылочного и текущего кадров на соответствующие поля средних значений, которые получают линейной НЧ-фильтрацией данных кадров; при компенсации фликкер-эффекта приводят локальную среднюю яркость ссылочного кадра к локальной средней яркости текущего, для чего делят ссылочный кадр на его скомпенсированное по движению поле средних значений и умножат на поле средних значений текущего кадра; при рекурсивном усреднении осуществляют коррекцию артефактов компенсации движения с помощью смешивания текущего отфильтрованного кадра с его исход-3 017302 ной версией на основе вычисления их пиксельного и структурного подобия. Возможен вариант исполнения заявляемого способа, в котором при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка (отбрасывание границ), используют разность между текущим и скомпенсированным по движению ссылочным кадром. Наиболее близким к используемому в настоящем изобретении методу оценки шума цифровых медицинских изображений (рентгенограмм) является способ, изложенный в публикации [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006], в соответствии с которым для построения по исходному изображению оценки зависящего от сигнала шума достаточно осуществить следующие этапы: оценить полезный сигнал с помощью низкочастотной фильтрации исходного изображения и вычислить разность между исходным изображением и его оценкой, получив тем самым изображение шума; отбросить тем или иным способом значения пикселей изображения шума, соответствующие резким изменениям (границы, одиночные "горячие" пиксели) в исходном изображении; разбить диапазон интенсивностей оценочного изображения на интервалы и для каждого такого интервала накопить значения пикселей изображения шума, соответствующие пикселям оценочного изображения; вычислить дисперсию шума в каждом интервале интенсивности по накопленным в данном интервале значениям пикселей изображения шума. Настоящее изобретение использует изложенные в упомянутой выше публикации принципы оценки шума. Наиболее существенные отличия состоят в том, что для каждого кадра серии изображений выполняют следующие операции: осуществляют удаление значений пикселей изображения шума, соответствующих резким изменениям в исходном изображении, с помощью морфологического выделения значений пикселей изображения шума, соответствующих границам на исходном изображении; выполняют робастную локальную линейную аппроксимацию интервальных оценок дисперсии шума, в результате которой получают табличную функцию, описывающую зависимость шума от интенсивности сигнала; вычисляют на основе оценочного изображения и найденной табличной функции карту шума в виде изображения, являющегося попиксельной оценкой дисперсии шума исходного цифрового изображения. Наиболее близким к используемому в настоящем изобретении методу учта и компенсации фликкер-эффекта является метод, описанный в публикации [Wong et al., Improved Flicker Removal ThroughMotion Vectors Compensation. In Third International Conference on Image and Graphics, pp. 552-555, 2004]. В этом методе мерцания моделируют гладкой функцией, мультипликативно искажающей значения интенсивностей кадров видеопоследовательности. Оценку данной функции предлагают осуществлять совместно с оценкой движения на основе поиска соответствия блоков, предлагают использовать при этом кусочно-постоянную искажающую функцию. Величину фликкер-эффекта в каждом блоке оценки движения текущего кадра полагают равной отношению средних значений этого блока и соответствующего блока ссылочного кадра, найденного в результате компенсации движения. Полученное кусочнопостоянное поле оценок фликкер-эффекта подвергают пороговой обработке и сглаживанию для удовлетворения условиям гладкости. В настоящем изобретении фликкер-эффект также моделируют в виде мультипликативной гладкой функции. Однако, компенсацию фликкер-эффекта осуществляют на основе приведения локальных средних предыдущего кадра к локальным средним текущего кадра посредством деления скомпенсированного по движению предыдущего отфильтрованного кадра (ссылочного кадра) на скомпенсированное по движению его поле средних значений и умножением на поле средних значений текущего кадра. Поле средних значений кадра получают с помощью его усреднения скользящим средним с использованием окна,апертура которого равна размеру блока оценки движения. Соответствующим образом модифицируется и карта шума ссылочного изображения. Движение в заявляемом изобретении оценивают и компенсируют на основе поиска соответствия блоков в изображениях серии [Wang et. al., Fast Algorithms for the Estimation of Motion Vectors. IEEETrans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]. Учт фликкер-эффекта, оказывающего существенное влияние на качество оценки векторов движения, осуществляют на основе простейшего подхода, при котором на этапе оценки движения текущее и ссылочное изображения серии делят на их поля средних значений, что позволяет осуществить своего рода нормировку по яркости. Помимо учта эффекта мерцаний интенсивности кадров оценка движения осуществляется с использованием пирамидальной схемы. Причина этого состоит в том, что при фильтрации в реальном времени серии медицинских изображений особую трудность представляет оценка и компенсация быстрых движений, например, в исследованиях сердца (коронарографии), когда используется сильное приближение к объекту исследования. В результате, для успешной компенсации движения, в особенности на низкой частоте смены кадров, область поиска подходящих блоков становится слишком велика. Поэтому для решения проблемы возрастающей вычислительной сложности при компенсации быстрых движений в изобретении используют-4 017302 пирамидальный алгоритм оценки движения, при котором текущий кадр и ссылочный кадр, который должен быть скомпенсирован по движению, раскладывают на несколько уровней в пирамиду изображений(не более трх уровней, зависит от типа исследования); оценку векторов перемещения блоков осуществляют на низкочастотных составляющих пирамиды поиском по небольшой области, при выборе радиуса которой учитывают число уровней в пирамиде; при переходе с нижнего на верхний уровень пирамиды размер блока компенсации движения масштабируют и уточняют вектор движения данного блока поиском в подобласти малого радиуса. Применение пирамидальной схемы позволяет существенно расширить область поиска подходящих блоков на ссылочном кадре. Для уменьшения эффекта блочности, типичного при компенсации движения на основе использования прямоугольных блоков, в изобретении блоки оценки движения рассматривают с перекрытиями [Orchard et. al., Overlapped Block Compensation: An Estimation-Theoretic Approach. IEEETransactions On Image Processing, Vol. 3, No. 5, September 1994, pp. 693-699]. Типичный размер перекрытия составляет четверть радиуса блока оценки движения. При построении скомпенсированного по движению кадра пиксели, соответствующие областям перекрытия блоков, усредняют с использованием гладкой весовой функции. Отличительной особенностью заявляемого изобретения является использование оригинальной методики компенсация артефактов алгоритма учета движения. Грубые ошибки компенсации движения,которые могут возникнуть из-за слишком малого радиуса области поиска, перекрытия объектов или резкой смены содержимого кадров, уменьшают с помощью смешивания отфильтрованного кадра с усредннной в пространстве изображения его версией. Веса смешивания вычисляют на основе пиксельного и структурного подобий текущего изображения скомпенсированному по движению изображению. В качестве меры пиксельного подобия используют билатеральное расстояние [Tomasi et. al., Bilateral Filteringfor Gray and Color Images, In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]. Для вычисления структурного подобия применяют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров [Buades et. al., Nonlocal image and movie denoising. Int. J. comput. Vision,76(2):123-139, 2008]. При вычислении финальной меры подобия кадров структурное и пиксельное подобия смешивают по линейному закону. В настоящем изобретении с целью обеспечения высокого уровня подавления шума и удовлетворения требованию выполнения в реальном времени используют рекурсивный фильтр Калмана с автоматической оценкой зависящего от сигнала шума, учтом возможного резкого изменения интенсивности последовательных кадров, оценкой и компенсацией движения, уменьшением артефактов компенсации движения и временного пересглаживания. При разделении данных задач между многоядерным процессором обычного на сегодняшний день настольного компьютера и графической картой с технологиейNvidia CUDA удалось получить результат, удовлетворяющий требованию фильтрации с высоким качеством в реальном времени. Подробное описание изобретения Заявляемое техническое решение, возможность его технической реализации и достижение указанного технического результата поясняется фиг. 1-4. На фиг. 1 показано устройство для получения рентгенограмм. На фиг. 2 показан график табличной функции, описывающей зависимость стандартного отклонения шума от интенсивности сигнала для кадра серии изображений сосудистого исследования, изображнного на фиг. 4. На фиг. 3 приведена карта стандартного отклонения шума изображения фиг. 4. На фиг. 4 показан исходный кадр серии изображений. На фиг. 5 приведн результат фильтрации данной серии изображений с 80-процентным уровнем подавления шума (параметр r в формуле (12) равен 0.8). Для лучшего восприятия фиг. 4 и 5 был применн фильтр повышения резкости. На фиг. 6 показано разностное изображение кадров фиг. 4 и 5. Получение рентгенограмм осуществляют, например, с помощью устройства, показанного на фиг. 1. Оно содержит рентгеновскую трубку 1, которая испускает пучок рентгеновского излучения 2. Пучок рентгеновского излучения 2 поступает на детектор 3. Детектор 3 содержит сцинтилляционный экран (на чертеже не показан) и матрицу фоточувствительных элементов (на чертеже не показаны). Сцинтилляционный экран оптически связан с поверхностью активной области матрицы фоточувствительных элементов. Пучок рентгеновского излучения 2 попадает на детектор 3, сцинтилляционный экран преобразует его в видимый свет, который преобразуется сенсорами детектора в цифровую форму. Предлагаемый в настоящем изобретении фильтр осуществляет подавление шума в несколько этапов. Далее рассмотрим каждый из этих этапов более подробно. Этап 1. Покадровая оценка дисперсии шума, зависящего от интенсивности сигнала. На данном этапе вначале осуществляют оценку изображения с помощью линейной низкочастотной фильтрации текущего кадра [Hensel et. al, Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]. Для обеспечения жстких требований к скорости выполнения в приложениях реального времени целесообразно использовать простейшую линейную фильтрацию (например, биномиальным фильтром). На основе полученного сглаженного изображения строят изображение шума - разницу между исходным и фильтрованным изображениями. Оценка изображения простейшими фильтрами неидеальна - границы оказываются пересглаженны-5 017302 ми. В результате при взятии разности между исходным изображением и сглаженным изображением получают изображение шума, содержащее помимо шумовых пикселей в гладких областях определнное количество пикселей, соответствующих резким изменениям (например, границы анатомических структур - далее нешумовые пиксели). Данные пиксели могут существенно исказить оцениваемое значение дисперсии и должны быть исключены из расчта. Для этого применяются различные методики [Foi et al.,Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEEMammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, September 2008], как правило, сводящиеся к пороговой обработке сглаженных производных исходного изображения, при этом величина порога определяется локальной оценкой отношения сигнал/шум. В областях изображения, содержащих большое количество деталей, такая оценка оказывается обычно неудовлетворительной. Поэтому в настоящем изобретении при отбрасывании границ предлагают более простой подход к выделению границ, не требующий вычисления производных и локальных оценок стандартного отклонения. Суть данного морфологического подхода к отбрасыванию нешумовых пикселей состоит в следующем: 1) изображение шума разбивают на две составляющие - бинарные изображения положительных и отрицательных изменений; 2) для выделения областей, соответствующих границам, полученные изображения подвергают морфологическим операциям эрозии с последующей дилатацией; 3) обработанные изображения объединяют для получения одного бинарного изображения - карты границ исходного изображения. С целью более полного сохранения тонких структур морфологические операции эрозии и дилатации выполняют масками малого размера (окно 22). При вычислении интервальных оценок дисперсии шума находят минимальную и максимальную интенсивности оценочного изображения (края диапазона интенсивностей), выбирают шаг разбиения,равный, например, 32 градациям серого цвета. Далее для каждого пикселя оценочного изображения находят интервал, которому принадлежит значение данного пикселя, и соответствующее значение пикселя изображения шума используют для вычисления оценки дисперсии шума в данном интервале (при этом исключают пиксели, соответствующие границам). При вычислении интервальной оценки дисперсии шума можно применить различные формулы, например обычную несмещнную оценку, или робастную оценку по формуле медианы абсолютных отклонений [Hensel et al., Robust and Fast Estimation of SignalDependent Noise in Medical X-Ray Image Sequences, Springer, 2006; Foi et al., Practical poissonian-gaussiannoise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754,2008]. На выходе данной процедуры получают табличную функцию, описывающую зависимость дисперсии шума от интенсивности сигнала. Неточности, возникающие при построении карты границ, могут привести к грубым ошибкам при вычислении интервальных оценок дисперсии. Поэтому данные оценки уточняют, используя технику итеративного отбрасывания выбросов [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noisein Medical X-Ray Image Sequences, Springer, 2006]: для каждого интервала итеративно исключают значения шумовых пикселей, превышающие по абсолютной величине порог, равный трм стандартным отклонениям шума, с последующим пересчтом оценки дисперсии шума в данном интервале. После того как вычислены интервальные оценки дисперсии шума, осуществляют оценку зависимости дисперсии шума от интенсивности. При параметрической оценке в принципе можно построить тем или иным способом (например, методом наименьших квадратов, минимизации функции правдоподобия,направленной оптимизации) оценку параметров модели шума (1). Возможен также учт нелинейностей сенсора [Foi et al., Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Imageand Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006], построение параметрической модели, адекватно описывающей зависимость дисперсии шума от интенсивности сигнала, в силу ряда факторов может вызвать серьезные затруднения. К данным факторам можно отнести нелинейности в работе сенсора, нелинейные предобработки исходных снимков (например, логарифмирование). Поэтому более выгодным с точки зрения удобства реализации и универсальности применения представляется подход, при котором на основе интервальных оценок дисперсии строится непараметрическая оценка зависимости шума от интенсивности. В настоящем изобретении используют непараметрический подход к построению искомой зависимости, при котором по полученным интервальным оценкам дисперсии шума создают интерполирующую табличную функцию. Данную табличную функцию формируют на основе робастной локальной линейной аппроксимации интервальных оценок дисперсии. Использование робастных методов позволяет дополнительно снизить влияние выбросов (грубых ошибок в интервальных оценках дисперсии), в то время как локальность аппроксимации обеспечивает повторение сложного хода кривой, описывающей зависимость шума от интенсивности. Таким образом, полученная табличная функция каждой интенсивности исходного изобра-6 017302 жения ставит в соответствие оценку дисперсии шума. Точками входа в таблицу могут служить, например, интенсивности оценочного изображения. На практике, в случае параметрической оценки шума, может быть использован подход, при котором применяется преобразование, стабилизирующее дисперсию шума исходного изображения [Starck et.al, Image Processing and Data Analysis: The Multiscale Approach. Cambridge University Press, 1998; Яне,Цифровая обработка изображений, М.: Техносфера, 2007]. При этом проблема фильтрации шума, зависящего от сигнала, сводится к задаче подавления аддитивного не зависимого от сигнала шума заданной дисперсии. В настоящем изобретении осуществляют непараметрическую оценку шума, поэтому при построении карты шума предлагают использовать следующий подход: на основе оценочного изображения и интерполирующей таблицы построить карту шума - изображение, каждый пиксель которого оценивает среднеквадратическое отклонение шума в соответствующем пикселе исходного изображения. Карта шума дат попиксельную оценку шума с достаточной для практического применения точностью. Возможен вариант исполнения заявляемого способа, в котором при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка (отбрасывание границ), используют разность между текущим и скомпенсированным по движению ссылочным кадром. Этап 2. Оценка и компенсация движения и фликкер-эффекта. В настоящем изобретении оценку движения осуществляют на основе поиска соответствия блоков, с учтом возможной вариативности интенсивности последовательных кадров. Для этого до применения процедуры поиска векторов движения текущий и ссылочный кадры делят на их поля средних значений. Поля средних получают фильтрацией соответствующих изображений усредняющим по окрестности фильтром, апертура которого равна размеру используемого блока оценки движения, например 1616 пикселей. При оценке и компенсации движения блоки оценки движения рассматривают с перекрытием, равным, например, четверти размера блока. При поиске соответствия блоков используют пирамидальный алгоритм оценки движения, при котором текущий кадр и ссылочный кадр, который должен быть скомпенсирован по движению по отношению к текущему кадру, раскладывают на несколько уровней в пирамиду изображений; оценку векторов перемещения блоков осуществляют на низкочастотных пирамидах изображений; при переходе с нижнего на верхний уровень пирамиды размеры блоков и их перекрытий масштабируют, а сам вектор движения данного блока уточняют поиском в области малого радиуса; при построении скомпенсированного по движению кадра области перекрытия блоков усредняют. Оценку движения от ссылочного кадра к текущему осуществляют на низкочастотных изображениях пирамиды. Используют не более трх уровней разложения, как правило два, поскольку при более глубоком разложении при оценке движения мелких объектов (тонкие сосуды) можно допустить грубые ошибки в оценке вектора движения. При переходе с нижнего на верхний уровень пирамиды размеры блока и их перекрытия удваивают и осуществляют уточняющий поиск вектора движения по области небольшого размера, например 3x3 пикселя. При расчте квадрата евклидовой метрики как меры подобия блоков используют технику раннего прекращения вычислений [Wang et al., Fast Algorithms for the Estimation ofMotion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]. Кроме того, применяют технику отсеивания неподвижных блоков из расчта оценки вектора движения: если выборочное где k - значение поростандартное отклонение блока оценки движения не превосходит величины- величина стандартного отклонения га (здесь для данного параметра используют значение от 2 до 3),шума кадра в центральном пикселе блока оценки движения. После оценки движения формируют скомпенсированный ссылочный кадр, при этом для уменьшения эффекта блочности перекрытия блоков усредняют. Усреднение перекрытий блоков осуществляют с использованием гладкого весового окна. Далее с целью компенсации фликкер-эффекта приводят локальные средние ссылочного кадра к локальным средним текущего кадра. Для этого делят скомпенсированный ссылочный кадр на его скомпенсированное по движению поле средних и умножают на поле средних значений текущего кадра согласно следующей формуле:- скомпенсированный по движению ссылочный кадр;- значение поля средних текущего кадра серии изображений в пикселе p;- скомпенсированное по движению поле средних значений ссылочного изображения;- результирующее, скомпенсированное по движению ссылочное изображение, локальная средняя интенсивность которого приведена к локальной средней интенсивности текущего изображения. Соответствующему преобразованию подвергают и карту дисперсии шума ссылочного кадра. Этап 3. Рекурсивное усреднение с коррекцией артефактов компенсации движения Грубые ошибки в компенсации движения, которые могут возникнуть из-за слишком малого радиуса области поиска, перекрытия объектов или резкой смены содержимого кадров, в настоящем изобретении уменьшают с помощью смешивания отфильтрованного кадра с усредннной в пространстве изображения-7 017302 его версией. Веса смешивания вычисляют на основе пиксельного и структурного подобий текущего изображения скомпенсированному по движению изображению. В качестве меры пиксельного подобия используют билатеральное расстояние [Tomasi et. al, Bilateral Filtering for Gray and Color Images, In Proc. 6thInt. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]. Для вычисления структурного подобия применяют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров[Buades et al., Nonlocal image and movie denoising. Int. J. comput. Vision, 76(2):123-139, 2008]. Пусть Ir(p), IC(P) обозначают значения интенсивностей ссылочного и текущего кадров в пикселе p соответственно. Тогда пиксельное подобие кадров вычисляют согласно следующей пороговой функции: где kp, Tp - параметры, определяющие характерную форму функции; 2(Ic) - величина дисперсии шума текущего кадра в пикселе p. В свою очередь, структурное подобие между этими кадрами вычисляют по формуле где ks, Ts - параметры, определяющие форму пороговой функции; P(Ic), P(Is) - квадратные блоки пикселей (патчи) соответствующих кадров с центром в пикселе с номером p; параметры P(Ic)-P(Ir)2 квадрат евклидова расстояния между блоками пикселей. Для вычисления результирующей меры подобия кадров, объединяющей пиксельное и структурное сходство кадров, используют соотношение Параметропределяет преобладание структурного веса над весом, вычисленным по интенсивности соответствующих пикселей. Данная форма вычисления смешанного веса, при значениях параметра ,близких к 1 (например, 0,9), позволяет при смешивании кадров уменьшить влияние одиночных импульсных значений интенсивности пикселей, одновременно восстанавливая правильное структурное содержание кадра (подобное содержанию исходного текущего неотфильтрованного кадра). Временное усреднение осуществляют рекурсивно с помощью одномерной Калмановской фильтрации для скомпенсированного по движению ссылочного кадра (кадр-предиктор) и текущего кадра (наблюдаемый кадр) в соответствии со следующими формулами 1. Оценка интенсивности по Калману: 3). Регулируемая сила подавления шума:- величина дисперсии шума кадра-предиктора (предыдущего кадра, ссылочного кадЗдесь ра), скомпенсированная по движению и нормированная по интенсивности, данное значение берут из карты шума; t(p) - величина дисперсии шума наблюдаемого кадра (текущего кадра);- интенсивность ссылочного кадра, скомпенсированного по движению и интенсивности (см. формулу (4; It(p) - интенсивность наблюдаемого кадра;- рекурсивная оценка интенсивности наблюдаемого кадра;- рекурсивная оценка дисперсии наблюдаемого кадра;W(p) - мера структурного и пиксельного подобия рекурсивно-усредннного кадра наблюдаемому кадру (см. формулу (5;Iet(p) - оценка интенсивности с помощью усреднения в пространстве наблюдаемого кадра биномиальным фильтром размера 33 пикселя; r - регулируемая сила подавления шума. В качестве меры ошибки наблюдаемого кадра используют рассчитанную карту дисперсии шума t(p), при этом карту дисперсии шума кадра-предиктора компенсируют по движению и интенсивности на основе найденных-8 017302 на этапе оценки движения векторов и рассчитанных полей средних значений. Лучший вариант осуществления изобретения Для оценки дисперсии шума текущего кадра вначале получают оценочное изображение и изображение шума, для чего исходное изображение I(x, y) фильтруют следующим низкочастотным линейным биномиальным фильтром размера 33: При этом получают сглаженное изображение Ie(x, y) = IН. Далее вычисляют изображение шумаNe(x, y) = I(х, у) - Ie(x, y). При отбрасывания границ из расчтов дисперсии на основе изображение шума Ne(x,y) формируют изображения положительных и отрицательных изменений Для выделения крупных областей, соответствующих границам, данные бинарные изображения последовательно подвергают морфологическим операциям эрозии и дилатации, используя маски размера 22: Затем эти изображения объединяют для получения карты границ исходного изображения Е(х, у) = В-e(x, y)В+е(x, y). Для вычисления интервальных оценок дисперсии шума находят минимальную Imin и максимальнуюImax интенсивности изображения Ie(x, y), выбирают шаг hM и разбивают диапазон интенсивностей на интервалы Mi с шагом hM (hM = 32). Для каждого пикселя изображения Ie(x, y) находят интервал, которому принадлежит значение этого пикселя, а соответствующее значение пикселя изображения Ne(x, у) используют для вычисления оценки дисперсии шума 2(i) в данном интервале Mi, при этом исключают значения пикселей, для которых значение карты границ Е(x, y) = 1. При вычислении интервальной оценки дисперсии шума применяют формулу несмещнной оценки дисперсии где Nie(j) - значение пикселя изображения шума Ne(x, y) из интервала Mi, ni - общее количество на- среднее значений шумовых копленных значений пикселей изображения шума в интервале пикселей в интервале Mi. Полученные интервальные оценки дисперсии 2(i) уточняют таким образом, что для каждого интервала Mi итеративно исключают значения шумовых пикселей, превышающие по абсолютной величине порог, равный трм стандартным отклонениям шума, с последующим пересчтом оценки дисперсии шума в данном интервале: где Nei[k+1] - набор значений шумовых пикселей в интервале Mi на итерации k. В дальнейших расчтах используются лишь те интервалы, в которых накоплено достаточное количество значений шумовых пикселей (например, не меньше 500). Кроме того, из рассмотрения исключаются те интервалы,среднее значение в которых существенно отличается от нуля (отклоняется от нуля более чем на половину шага сетки интервалов hM), поскольку в таких интервалах с высокой вероятностью доминируют остаточные значения нешумовых пикселей. Для оценки зависимости дисперсии шума от интенсивности по полученным интервальным оценкам дисперсии шума (отмечены кружками на фиг. 2) создают интерполирующую табличную функцию. Данную табличную функцию формируют на основе робастной локальной линейной аппроксимации интервальных оценок дисперсии. Для этого на сетке интервалов Mi выбирают шаг hI и радиус rI аппроксимации (которые можно сделать зависимыми от количества точек ni в интервалах Mi). Значения полученной на предыдущем этапе табличной функции дисперсии шума от интенсивности сигнала аппроксимируют по следующему правилу: где k - номер узла аппроксимации; mk - центр интервала Mi(khI) ; параметры a, b вычисляют на основе минимума суммы абсолютных отклонений [Press et al., Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press,1999] Полученную в результате такого процесса таблицу значений интерполируют на всм диапазоне интенсивностей [Imin,Imax], т.е. интерполирующую таблицу искомой зависимости дисперсии шума от интенсивности сигнала (сплошная линия фиг. 2). При построении карты шума на основе оценочного изображения и интерполирующей таблицы строят карту шума - изображение, каждый пиксель которого оценивает дисперсию шума в соответствующем пикселе исходного изображения (фиг. 3). Оценку и компенсацию движения осуществляют на основе поиска соответствия блоков, с учтом возможной вариативности интенсивности последовательных кадров. Для этого перед оценкой движения текущий и ссылочный кадры делят на их поля средних значений. Данные поля средних получают фильтрацией соответствующих изображений усредняющим по окрестности фильтром, размер которого равен размеру используемого блока оценки движения, например 1616 пикселей. При оценке движения блоки рассматривают с перекрытием, равным, например, 4 пикселям. При поиске соответствия блоков используют пирамидальную схему: текущий кадр и ссылочный кадр, который должен быть скомпенсирован по движению по отношению к текущему кадру, раскладывают на несколько уровней в пирамиду изображений; оценку векторов перемещения блоков осуществляют на низкочастотных пирамидах изображений; при переходе с нижнего на верхний уровень пирамиды размеры блоков и их перекрытий масштабируют, а сам вектор движения данного блока уточняют поиском в области малого радиуса; при построении скомпенсированного по движению кадра области перекрытия блоков усредняют. Оценку движения от ссылочного кадра к текущему осуществляют на низкочастотных изображениях пирамиды. Используют не более трх уровней разложения, как правило, два уровня, поскольку при более глубоком разложении, в особенности при оценке движения мелких объектов (сосуды), на грубых уровнях разрешения тонкие структуры пропадут, что приведт к грубым ошибкам в оценке векторов движения. При переходе с нижнего на верхний уровень пирамиды размеры блока и их перекрытия удваивают и осуществляют уточняющий поиск вектора движения по области небольшого размера, например 33 пикселя. При расчте квадрата евклидовой метрики как меры подобия блоков используют технику раннего прекращения вычислений [Wang et. al., Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans.Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]. Кроме того, применяют технику отсеивания неподвижных блоков из расчта оценки вектора движения: если выборочное стандартное отклонение блока оценки движения не превосходит величины, где k - значение порога (здесь для данного параметра используют значение от 2 до 3),- величина стандартного отклонения шума кадра в центральном пикселе блока оценки движения. После оценки движения формируют скомпенсированный ссылочный кадр, при этом для уменьшения эффекта блочности перекрытия блоков усредняют. Аналогичным образом компенсируют карту шума ссылочного кадра. Усреднение перекрытий блоков осуществляют с использованием гладкого весового окна, например косинусного. Для компенсации фликкер-эффекта приводят локальные средние ссылочного кадра к локальным средним текущего кадра. Для этого делят скомпенсированный ссылочный кадр на его скомпенсированное по движению поле средних и умножают на поле средних значений текущего кадра согласно формуле(4). Карту шума ссылочного кадра подвергают аналогичной нормировке. Рекурсивная фильтрация текущего кадра Временное усреднение скомпенсированных по движению и интенсивности кадров осуществляют с помощью рекурсивной фильтрации, используя одномерный (попиксельный) фильтр Калмана. Расчт ведут по формулам (6)-(12). Для подавления артефактов временного усреднения, обусловленных ошибками в компенсации движения, используют смешивание фильтрованного нормализованного кадра с его исходной версией (формулы (5), (9), (10. Веса смешивания (5) находят на основе вычисления яркостного и структурного подобия текущего кадра скомпенсированному по движению и интенсивности ссылочному кадру, при этом в качестве меры структурного подобия используют евклидово расстояние между блоками пикселей текущего и скомпенсированного кадров. Расчт структурного подобия вычисляют скользящим блоком, размер которого выбирают равным, например, 1111 пикселей. Расчт пиксельного подобия вычисляют как расстояние в пространстве интенсивностей между текущим кадром и ссылочным,скомпенсированным по движению кадром. Силу подавления шума r в формулах (11), (12) на реальных- 10017302 сериях берут равной, например, r = 0,8, что соответствует пятикратному подавлению уровня шума. Фиг. 4-6 демонстрируют результат применения фильтра к реальной серии (ангиографическое исследование): фиг. 4 показывает исходный кадр серии, фиг. 5 - фильтрованный кадр, фиг. 6 - разность между кадрами фиг. 4-5; для улучшения восприятия исходный и отфильтрованный кадры были подвергнуты преобразованию повышения резкости. ФОРМУЛА ИЗОБРЕТЕНИЯ 1. Способ подавления шума серий цифровых рентгенограмм, включающий получение серии цифровых рентгенограмм, покадровую оценку дисперсии шума, зависящего от сигнала, оценку и компенсацию движения, оценку и компенсацию фликкер-эффекта, рекурсивное усреднение кадров, отличающийся тем, что при покадровой оценке дисперсии шума, зависящего от интенсивности сигнала, осуществляют морфологическое удаление значений пикселей снимка шума, соответствующих границам в исходном снимке; получают табличную зависимость шума от интенсивности сигнала с помощью робастной кусочно-линейной аппроксимации интервальных оценок дисперсии шума; вычисляют карту шума в виде попиксельной оценки дисперсии шума исходного цифрового изображения; при оценке и компенсации движения и фликкер-эффекта используют пирамидальную схему поиска соответствия блоков; блоки оценки движения рассматривают с перекрытием, при компенсации движения перекрытия блоков усредняют с использованием гладкого весового окна; учт фликкер-эффекта осуществляют делением ссылочного и текущего кадров на соответствующие поля средних значений, которые получают линейной НЧфильтрацией данных кадров; для компенсации фликкер-эффекта приводят локальную среднюю яркость ссылочного кадра к локальной средней яркости текущего, для чего делят ссылочный кадр на его скомпенсированное по движению поле средних значений и умножают на поле средних значений текущего кадра; при рекурсивном усреднении с коррекцией артефактов компенсации движения осуществляют смешивание текущего отфильтрованного кадра с его исходной версией на основе вычисления пиксельного и структурного подобия данных кадров. 2. Способ по п.1, отличающийся тем, что при оценке дисперсии шума для уменьшения влияния на вычисляемую оценку пикселей снимка, соответствующих границам снимка, используют разность между текущим и скомпенсированным по движению ссылочным кадром.