способ оценки шума цифровых рентгенограмм

Формула изобретения

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

Описание изобретения к патенту

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

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

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

Случайные флуктуации числа поглощенных фотонов называют квантовым шумом, или шумом фотонов. В современных детекторах основным источником шума изображений является шум фотонов. К дополнительным источникам шума относят шумы детекторной системы: шум чтения, тепловой шум, шумы усилителей, шумы квантования и другие [5]. Общий эффект данных источников шума можно моделировать распределенной по Гауссу случайной величиной [1-3]. Согласно общепринятой модели, в линейных электронных схемах дисперсия суммарного шума цифровых изображений (шум фотонов плюс шумы дополнительных источников) линейно зависит от величины полезного сигнала [11]:

где I(p) - уровень интенсивности сигнала в пикселе p.

Проблеме оценки зависящего от сигнала шума цифровых изображений посвящено достаточно много публикаций [1, 3, 4, 6-8, 10]. Так, в работе [6] представлен непараметрический метод оценки шума, при этом особый акцент делается на создании алгоритма, способного работать в реальном времени.

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

Наиболее близким к заявляемому является способ, изложенный в публикации [2, 6], в соответствии с которым для построения оценки зависящего от сигнала шума по исходному цифровому изображению достаточно осуществить следующие этапы:

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

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

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

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

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

Технический результат в способе оценки шума цифровых рентгенограмм, включающем: получение исходного изображения; получение оценочного изображения с помощью низкочастотной фильтрации исходного изображения; построение изображения шума как разности между исходным и оценочным изображением; удаление значений пикселей изображения шума, соответствующих резким изменениям в исходном изображении; разбиение диапазона интенсивности оценочного изображения на интервалы, при этом каждый пиксель оценочного изображения принадлежит определенному интервалу; накопление для каждого интервала значений пикселей изображения шума, соответствующих пикселям оценочного изображения; вычисление интервальных оценок дисперсии шума по накопленным в данном интервале значениям пикселей шумового изображения; уточнение интервальных оценок дисперсии путем отбрасывания значений шумовых пикселей по критерию 3-сигма, достигается тем, что

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

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

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

Выполнение алгоритма состоит из нескольких этапов:

1) оценка полезного сигнала с помощью низкочастотной фильтрации исходного изображения и получение изображения шума с помощью вычисления разности между исходным изображением и его оценкой;

2) морфологическое отбрасывание значений пикселей изображения шума, соответствующих границам в исходном изображении;

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

4) уточнение интервальных оценок дисперсии с помощью отбрасывания значений шумовых пикселей по критерию 3-сигма;

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

6) построение карты шума по оценочному изображению и найденной табличной зависимости шума от интенсивности.

Один из возможных способов получения рентгенограмм осуществляют с помощью устройства, показанного на фиг.1. Устройство содержит рентгеновскую трубку 1, которая вырабатывает пучок 3 рентгеновского излучения. Пучок 3 рентгеновского излучения пропускают через объект 2, подвергаемый рентгеновскому обследованию. Излучение воспринимается детектором 4. Детектор 4 содержит экран, содержащий сцинтилляционный слой, и матрицу фоточувствительных сенсоров. Экран оптически связан с поверхностью активной области фоточувствительных сенсоров. Падающее рентгеновское излучение 3 сначала преобразуется в видимый свет в сцинтилляционном слое, интенсивность которого регистрируют матрицей фоточувствительных сенсоров. С матрицы фоточувствительных сенсоров считывают цифровую информацию и выводят на монитор в виде изображения.

На этапе построения оценочного изображения и получения изображения шума оценку полезного сигнала исходного изображения осуществляют, например, с помощью фильтрации исходной рентгенограммы (линейной низкочастотной свертки, медианной фильтрации и т.п.) [2, 6]. Для обеспечения жестких требований к скорости выполнения в приложениях реального времени целесообразно использовать простейшую линейную фильтрацию (например, биномиальным фильтром). На основе полученного сглаженного изображения строят изображение шума - разница между исходным и фильтрованным изображениями.

Оценка изображения простейшими фильтрами неидеальна - границы оказываются пересглаженными. В результате при взятии разности между исходным изображением и сглаженным изображением получают изображение шума, содержащее помимо шумовых пикселей в гладких областях определенное количество пикселей, соответствующих резким изменениям (границы анатомических структур - далее нешумовые пиксели). Данные пиксели должны быть исключены при расчете статистики шума, поскольку они могут существенно исказить вычисляемое значение дисперсии шума. Для этого применяются различные методики [3, 8], как правило, сводящиеся к пороговой обработке сглаженных производных исходного изображения, при этом величина порога определяется локальной оценкой отношения сигнал/шум. В областях изображения, содержащих большое количество деталей, такая оценка оказывается обычно неудовлетворительной. Поэтому в настоящем изобретении на этапе отбрасывания границ предлагают более простой подход к выделению границ, не требующий вычисления производных и локальных оценок стандартного отклонения. Суть данного морфологического подхода к отбрасыванию нешумовых пикселей состоит в следующем:

1) изображение шума разбивают на две составляющие - бинарные изображения положительных и отрицательных изменений;

2) для выделения областей, соответствующих границам, полученные изображения подвергают морфологическим операциям эрозии с последующей дилатацией;

3) обработанные изображения объединяют для получения одного бинарного изображения - карты границ исходного изображения.

С целью более полного сохранения тонких структур морфологические операции эрозии и дилатации выполняют масками малого размера (окно 2×2).

На этапе вычисления интервальных оценок дисперсии шума вычисляют минимальную и максимальную интенсивности оценочного изображения (края диапазона интенсивностей), выбирают шаг разбиения, равный, например, 32 градациям серого цвета. Далее для каждого пикселя оценочного изображения находят интервал, которому принадлежит значение данного пикселя, и соответствующее значение пикселя изображения шума используют для вычисления оценки дисперсии шума в данном интервале (при этом исключают пиксели, соответствующие границам). При вычислении интервальной оценки дисперсии шума можно применить различные формулы, например, обычную несмещенную оценку, или робастную оценку по формуле медианны абсолютных отклонений [2, 3, 6]. На выходе данной процедуры получают табличную функцию, описывающую зависимость дисперсии шума от интенсивности сигнала.

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

После того как вычислены интервальные оценки дисперсии шума, наступает этап оценки зависимости дисперсии шума от интенсивности. При параметрической оценке в принципе можно построить тем или иным способом (например, методом наименьших квадратов, минимизации функции правдоподобия, направленной оптимизации) оценку параметров модели шума (1). Возможен также учет нелинейностей сенсора [2]. Однако, как отмечено в работе [6], построение параметрической модели, адекватно описывающей зависимость дисперсии шума от интенсивности сигнала, в силу ряда факторов может вызвать серьезные затруднения. К данным факторам можно отнести нелинейности в работе сенсора, нелинейные предобработки исходных снимков (например, логарифмирование). Поэтому более выгодным с точки зрения удобства реализации и универсальности применения представляется подход, при котором на основе интервальных оценок дисперсии строится непараметрическая оценка зависимости шума от интенсивности.

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

На практике, в случае параметрической оценки шума, может быть использован подход, при котором применяется преобразование, стабилизирующее дисперсию шума исходного изображения [2, 3, 9, 11]. При этом проблема фильтрации шума, зависящего от сигнала, сводится к задаче подавления аддитивного независимого от сигнала шума постоянной (заданной) дисперсии. В настоящем изобретении осуществляют непараметрическую оценку шума, поэтому на этапе построения карты шума предлагают использовать следующий подход: на основе оценочного изображения и интерполирующей таблицы построить карту шума - изображение, каждый пиксель которого оценивает среднеквадратическое отклонение шума в соответствующем пикселе исходного изображения. Карта шума дает попиксельную оценку шума с достаточной для практического применения точностью. Использование как можно более точной информации об уровне шума на изображении позволяет существенно улучшить качество алгоритмов шумоподавления.

Заявляемый способ осуществляют следующим образом. На этапе построения оценочного изображения и получения изображения шума исходное изображение I(x, y) фильтруют следующим низкочастотным линейным биномиальным фильтром размера 3×3:

.

При этом получают сглаженное изображение I_e(x, y)=I*H. Далее вычисляют изображение шума N _e(x, y)=I(x, y)-I_e(x, y).

На этапе отбрасывания границ, используя изображение шума N_e (x, y), формируют два изображения положительных и отрицательных изменений

Для выделения крупных областей, соответствующих границам, данные бинарные изображения последовательно подвергают морфологическим операциям эрозии и дилатации, используя маски размера 2×2:

,

.

Затем эти изображения объединяют для получения карты границ исходного изображения .

На этапе вычисления интервальных оценок дисперсии шума вычисляют минимальную I_min и максимальную I_mах интенсивности изображения I_e(x, y), выбирают шаг h_M и разбивают диапазон интенсивностей на интервалы M_i, с шагом h_M (h_M =32). Для каждого пикселя изображения I_е(x, y) находят интервал, которому принадлежит значение этого пикселя, а соответствующее значение пикселя изображения N_e(x, y) используют для вычисления оценки дисперсии шума ²(i) в данном интервале M_i, при этом исключают значения пикселей, для которых значение карты границ Е(x, y)=1. При вычислении интервальной оценки дисперсии шума применяют формулу несмещенной оценки дисперсии

где - значение пикселя изображения шума N_e(x, y) из интервала M_i, n_i - общее количество накопленных значений пикселей изображения шума в интервале M _i, - среднее значений шумовых пикселей в интервале M _i.

На этапе уточнения интервальных оценок дисперсии величины ²(i) уточняют таким образом, что для каждого интервала M_i итеративно исключают значения шумовых пикселей, превышающие по абсолютной величине порог, равный трем стандартным отклонениям шума, с последующим пересчетом оценки дисперсии шума в данном интервале:

где - набор значений шумовых пикселей в интервале M_i на итерации k. В дальнейших расчетах используются лишь те интервалы, в которых накоплено достаточное количество значений шумовых пикселей (например, не меньше 500). Кроме того, из рассмотрения исключаются те интервалы, среднее значение в которых существенно отличается от нуля (отклоняется от нуля более чем на половину шага сетки интервалов h_M ), поскольку в таких интервалах с высокой вероятностью доминируют остаточные значения нешумовых пикселей.

На этапе оценки зависимости дисперсии шума от интенсивности по полученным интервальным оценкам дисперсии шума создают интерполирующую табличную функцию. Данную табличную функцию формируют на основе робастной локальной линейной аппроксимации интервальных оценок дисперсии. Для этого на сетке интервалов M_i выбирают шаг h _I и радиус r_I аппроксимации (которые можно сделать зависимыми от количества точек n_i в интервалах M _i). Значения полученной на предыдущем этапе табличной функции дисперсии шума от интенсивности сигнала аппроксимируют по следующему правилу:

,

где k - номер узла аппроксимации; m _k - центр интервала M_i(k·h_I); параметры а, b вычисляют на основе минимума суммы абсолютных отклонений [12]

Полученную в результате такого процесса таблицу значений интерполируют на всем диапазоне интенсивностей [I _min, I_max], т.е. получают интерполирующую таблицу искомой зависимости дисперсии шума от интенсивности сигнала.

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

На фиг.2 показан фрагмент реального рентгеновского изображения, иллюстрирующий нелинейный отклик сенсора, приведшего к переэкспонированию.

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

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

На фиг.5 приведена карта шума - изображение, каждый пиксель которого содержит оценку значения СКО шума.

Литература

1. Bosco A., Battiato S., Bruna A.R., Rizzo R. Noise reduction for cfa image sensors exploiting hvs behavior. Sensors 9(3), 1692-1713 (2009).

2. Foi A., Trimeche M., Katkovnik V., Egiazarian K. Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEE Transactions on 17, 1737-1754 (October 2008).

3. Foi A. Practical denoising of clipped or overexposed noisy images. Proc. 16th European Signal Process. Conf., EUSIPCO 2008, Lausanne, Switzerland, August 2008.

4. Forstner W. Image preprocessing for feature extraction in digital intensity, color and range images. In Springer Lecture Notes on Earth Sciences, 1998.

5. Gino M. Noise, Noise, Noise.

http://www.astrophys-assist.com/educate/noise/noise.htm

6. Hensel M., Lundt В., Pralow Т., Grigat R.-R. Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences. In Handels, H., et al., ed.: Bildverarbeitung fur die Medizin 2006: Algorithmen, Systeme, Anwendun-gen. Springer (2006) 46-50.

7. Liu С., Szeliski R., Kang S.B., Zitnick C.L., Freeman W.T. Automatic estimation and removal of noise from a single image. IEEE Transactions on Pattern Analysis and Machine Intelligence 30(2), 299{314 (2008).

8. Salmeri M., Mencattini A., Rabottino G., Lojacono R. Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, September 2008.

9. Starck J., Murtagh F., Bijaoui A. Image Processing and Data Analysis: The Multiscale Approach. Cambridge University Press, 1998.

10. Waegli B. Investigations into the Noise Characteristics of Digitized Aerial Images. In: Int. Arch. For Photogr. And Remote Sensing, Vol.32-2, pages 341-348, 1998.

11. Яне Б. Цифровая обработка изображений. M.: Техносфера, 2007, 583 с.

12. W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery. Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1992.