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

Реферат: Спосіб оцінки залежного від сигналу шуму цифрових медичних рентгенограм стосується галузі обробки цифрових зображень і може бути використаний для рішення задачі обробки цифрових зображень, отриманих за допомогою високоенергетичного випромінювання, у тому числі рентгенівського. Спосіб оцінки залежного від сигналу шуму цифрових медичних рентгенограм включає: одержання вихідного зображення, низькочастотну фільтрацію вихідного зображення для одержання зображення, що оцінюється, побудову зображення шуму як різниці між вихідним і оцінним зображенням; видалення за допомогою морфологічної фільтрації значень пікселів зображення шуму, що відповідають границям у вихідному зображенні, розбивку діапазону інтенсивності зображення, що оцінюється, на інтервали, при цьому кожний піксель зображення, що оцінюється, належить певному інтервалу, накопичення для кожного інтервалу значень пікселів зображення шуму, що відповідають пікселям зображення, що оцінюється, обчислення інтервальних оцінок дисперсії шуму за накопиченими у даних інтервалах значеннями пікселів шумового зображення, уточнення інтервальних оцінок дисперсії відкиданням значень шумових пікселів за критерієм 3-сигма; робастну локальну лінійну апроксимацію інтервальних оцінок UA 104164 C2 (12) UA 104164 C2 дисперсії шуму, у результаті якої одержують табличну функцію, що описує залежність шуму від інтенсивності сигналу, обчислення на основі зображення, що оцінюється, й побудованої табличної залежності дисперсії шуму від інтенсивності сигналу карти шуму у вигляді зображення, кожний піксель якого оцінює дисперсію шуму вихідного цифрового зображення. Технічним результатом способу є підвищення швидкості і якості оцінки шуму цифрових рентгенограм. UA 104164 C2 5 10 15 Винахід стосується галузі обробки цифрових зображень і може бути використаний для рішення задач обробки цифрових зображень, отриманих за допомогою високоенергетичного випромінювання, у тому числі рентгенівського. Більш конкретно, даний винахід призначений для оцінки залежного від сигналу шуму цифрових рентгенограм. У теперішній час у цифровій рентгенодіагностиці застосовуються різноманітні алгоритми обробки зображень. Такі способи як підвищення різкості або сегментація анатомічних структур можуть використовувати інформацію про рівень шуму на зображенні. Крім того, практично кожний якісний спосіб заглушення шуму використовує як параметр дисперсію шуму. Тому важливо розглянути наступну задачу: маючи у своєму розпорядженні тільки вихідне цифрове зображення, визначити рівень шуму цього зображення. Дана задача ускладнюється тим, що дисперсія шуму цифрового зображення істотно залежить від інтенсивності корисного сигналу. Під шумом звичайно розуміють випадкове відхилення обмірюваного значення деякої фізичної величини від її точного значення. У цифровій рентгенографії детектор вимірює інтенсивність рентгенівського випромінювання - кожний осередок працює як суматор фотонів, накопичуючи за час експонування в середньому N електронів шляхом поглинання фотонів. Кількість N накопичених за час експонування в осередку детектора електронів, породжених у процесі поглинання фотонів, можна моделювати розподіленою за законом Пуассона випадковою величиною:   20 25 30 35 40 45 50 55 n N PN  n  exp  N n! , n≥0. Випадкові флуктуації числа поглинених фотонів називають квантовим шумом, або шумом фотонів. У сучасних детекторах основним джерелом шуму зображень є шум фотонів. До додаткових джерел шуму відносять шуми детекторної системи: шум читання, тепловий шум, шуми підсилювачів, шуми квантування й інші [5]. Загальний ефект даних джерел шуму можна моделювати розподіленою за Гаусом випадковою величиною [1-3]. Відповідно до загальноприйнятої моделі, у лінійних електронних схемах дисперсія сумарного шуму цифрових зображень (шум фотонів плюс шуми додаткових джерел) лінійно залежить від величини корисного сигналу [11]: 2 σ (I(p))=aI(p)+b, (1) де I(р) - рівень інтенсивності сигналу в пікселі р. Проблемі оцінки залежного від сигналу шуму цифрових зображень присвячено досить багато публікацій [1, 3, 4, 6-8, 10]. Так, у роботі [6] представлений непараметричний спосіб оцінки шуму, при цьому особливий акцент ставиться на створенні алгоритму, здатного працювати в реальному часі. У публікації [3] розглядається двопараметричний підхід до оцінки шуму. Шум вихідного цифрового зображення (отриманого безпосередньо з детектора й не пройденого через нелінійні перетворення, такі як гамма-корекція й т.п.) моделюється як адитивна стосовно сигналу випадкова величина, дисперсія якої залежить від сигналу відповідно до закону (1). Цей спосіб побудови модельних кривих дисперсії шуму, ураховує нелінійності в роботі сенсора, що призводять до недоекспонування й переекспонування, тобто при порушенні закону (1) на краях динамічного діапазону. Найбільш близьким до способу, що заявляється, є спосіб, описаний у публікації [2, 6], відповідно до якого для побудови оцінки залежного від сигналу шуму за вихідним цифровим зображенням досить здійснити наступні етапи: оцінити корисний сигнал за допомогою низькочастотної фільтрації вихідного зображення й обчислити різницю між вихідним зображенням і його оцінкою, одержавши тим самим зображення шуму; відкинути тим або іншим способом значення пікселів зображення шуму, що відповідають різким змінам (границі, одиночні "гарячі" пікселі) у вихідному зображенні; розбити діапазон інтенсивностей зображення, що оцінюється, на інтервали й для кожного такого інтервалу накопичити значення пікселів зображення шуму, що відповідають пікселям зображення, що оцінюється; обчислити дисперсію шуму в кожному інтервалі інтенсивності за накопиченими уданому інтервалі значеннями пікселів зображення шуму. Даний винахід використовує викладені у вищезгаданих публікаціях принципи оцінки шуму. Задачею винаходу, що заявляється, є створення способу, більш ефективного за швидкістю і якістю оцінки шуму цифрових рентгенограм. Технічний результат у способі оцінки шуму цифрових рентгенограм, що включає: одержання вихідного зображення; одержання зображення, що оцінюється, за допомогою низькочастотної 1 UA 104164 C2 5 10 15 20 25 30 35 40 45 50 55 60 фільтрації вихідного зображення; побудову зображення шуму як різниці між вихідним і оцінним зображенням; видалення значень пікселів зображення шуму, що відповідають різким змінам у вихідному зображенні; розбивка діапазону інтенсивності зображення, що оцінюється, на інтервали, при цьому кожний піксель зображення, що оцінюється, належить певному інтервалу; накопичення для кожного інтервалу значень пікселів зображення шуму, що відповідають пікселям зображення, що оцінюється; обчислення інтервальних оцінок дисперсії шуму за накопиченим у даному інтервалі значенням пікселів шумового зображення; уточнення інтервальних оцінок дисперсії шляхом відкидання значень шумових пікселів за критерієм 3сигма, досягається за рахунок того, що: здійснюють видалення значень пікселів зображення шуму, що відповідають різким змінам у вихідному зображенні, за допомогою морфологічного виділення значень пікселів зображення шуму, що відповідають границям на вихідному зображенні; виконують робастну локальну лінійну апроксимацію інтервальних оцінок дисперсії шуму, у результаті якої одержують табличну функцію, що описує залежність шуму від інтенсивності сигналу; обчислюють на основі зображення, що оцінюється, й знайденої табличної функції карту шуму у вигляді зображення, що є попіксельною оцінкою дисперсії шуму вихідного цифрового зображення. Виконання алгоритму складається з декількох етапів: 1) оцінка корисного сигналу за допомогою низькочастотної фільтрації вихідного зображення й одержання зображення шуму за допомогою обчислення різниці між вихідним зображенням і його оцінкою; 2) морфологічне відкидання значень пікселів зображення шуму, що відповідають границям у вихідному зображенні; 3) розбивка діапазону інтенсивностей зображення, що оцінюється, на інтервали й обчислення дисперсії шуму для кожного такого інтервалу по накопиченим у даному інтервалі значенням пікселів зображення шуму; 4) уточнення інтервальних оцінок дисперсії за допомогою відкидання значень шумових пікселів за критерієм 3-сигма; 5) робастна локальна лінійна апроксимація інтервальних оцінок дисперсії шуму для побудови табличної функції, що описує залежність шуму від інтенсивності сигналу; 6) побудова карти шуму за зображенням, що оцінюється, й знайденою табличною залежністю шуму від інтенсивності. Один з можливих способів одержання рентгенограм здійснюють за допомогою пристрою, показаного на фіг. 1. Пристрій містить рентгенівську трубку 1, що виробляє пучок 3 рентгенівського випромінювання. Пучок 3 рентгенівського випромінювання пропускають крізь об'єкт 2, що піддається рентгенівському обстеженню. Випромінювання сприймається детектором 4. Детектор 4 містить екран, що містить сцинтиляційний шар, і матрицю фоточутливих сенсорів. Екран оптично пов'язаний з поверхнею активної області фоточутливих сенсорів. Падаюче рентгенівське випромінювання 3 спочатку перетвориться у видиме світло в сцинтиляційному шарі, інтенсивність якого реєструють матрицею фоточутливих сенсорів. З матриці фоточутливих сенсорів зчитують цифрову інформацію й виводять на монітор у вигляді зображення. На етапі побудови зображення, що оцінюється, й одержання зображення шуму оцінку корисного сигналу вихідного зображення здійснюють, наприклад, за допомогою фільтрації вихідної рентгенограми (лінійної низькочастотної згортки, медіанної фільтрації й т.п.) [2, 6]. Для забезпечення жорстких вимог до швидкості виконання в додатках реального часу доцільно використовувати найпростішу лінійну фільтрацію (наприклад, біноміальним фільтром). На основі отриманого згладженого зображення будують зображення шуму - різниця між вихідним і фільтрованим зображеннями. Оцінка зображення найпростішими фільтрами неідеальна - границі виявляються перезгладженими. У результаті при узятті різниці між вихідним зображенням і згладженим зображенням одержують зображення шуму, що містить крім шумових пікселів у гладких областях певну кількість пікселів, що відповідають різким змінам (границі анатомічних структур далі нешумові пікселі). Дані пікселі повинні бути виключені при розрахунку статистики шуму, оскільки вони можуть істотно спотворити значення дисперсії шуму, що обчислюється. Для цього застосовуються різні методики [3, 8], які, як правило, зводяться до граничної обробки згладжених похідних вихідного зображення, при цьому величина порогу визначається локальною оцінкою відношення сигнал/шум. В областях зображення, що містять велику кількість деталей, така оцінка виявляється звичайно незадовільною. Тому, у даному винаході на 2 UA 104164 C2 5 10 15 20 25 30 35 40 45 50 55 етапі відкидання границь пропонують простіший підхід до виділення границь, що не потребує обчислення похідних і локальних оцінок стандартного відхилення. Сутність даного морфологічного підходу до відкидання нешумових пікселів полягає в наступному: 1) зображення шуму розбивають на дві складових - бінарних зображення позитивних і негативних змін; 2) для виділення областей, що відповідають границям, отримані зображення піддають морфологічним операціям ерозії з наступною дилатацією; 3) оброблені зображення поєднують для одержання одного бінарного зображення - карти границь вихідного зображення. З метою більш повного збереження тонких структур морфологічні операції ерозії й дилатації виконують масками малого розміру (вікно 2×2). На етапі обчислення інтервальних оцінок дисперсії шуму обчислюють мінімальну й максимальну інтенсивності зображення, що оцінюється, (краї діапазону інтенсивностей), вибирають крок розбивки рівний, наприклад, 32 градаціям сірого кольору. Далі для кожного пікселя зображення, що оцінюється, знаходять інтервал, якому належить значення даного пікселя, і відповідне значення пікселя зображення шуму використовують для обчислення оцінки дисперсії шуму в даному інтервалі (при цьому виключають пікселі, що відповідають границям). При обчисленні інтервальної оцінки дисперсії шуму можна застосувати різні формули, наприклад, звичайну незміщену оцінку, або робастну оцінку за формулою медіани абсолютних відхилень [2, 3, 6]. На виході даної процедури одержують табличну функцію, що описує залежність дисперсії шуму від інтенсивності сигналу. Неточності, що виникають при побудові карти границь, є перешкодою в акуратному виключенні нешумових пікселів з розрахунку статистики, що може призвести до грубих помилок при обчисленні інтервальних оцінок дисперсії. Тому на етапі уточнення інтервальних оцінок дисперсії дані оцінки уточнюють, використовуючи техніку ітеративного відкидання викидів [6]: для кожного інтервалу ітеративно виключають значення шумових пікселів, що перевищують по абсолютній величині поріг, рівним трьом стандартним відхиленням шуму, з наступним перерахуванням оцінки дисперсії шуму в даному інтервалі. Після того як обчислено інтервальні оцінки дисперсії шуму, наступає етап оцінки залежності дисперсії шуму від інтенсивності. При параметричній оцінці в принципі можна побудувати тим або іншим способом (наприклад, методом найменших квадратів, мінімізації функції правдоподібності, спрямованої оптимізації) оцінку параметрів моделі шуму (1). Можливий також облік нелінійностей сенсора [2]. Однак, як відзначено в роботі [6], побудова параметричної моделі, що адекватно описує залежність дисперсії шуму від інтенсивності сигналу, в силу ряду факторів може викликати серйозні ускладнення. До даних факторів можна віднести нелінійності в роботі сенсора, нелінійні попередні обробки вихідних знімків (наприклад, логарифмування). Тому більш вигідним з погляду зручності реалізації й універсальності застосування представляється підхід, при якому на основі інтервальних оцінок дисперсії будується непараметрична оцінка залежності шуму від інтенсивності. У даному винаході використовується непараметричний підхід до побудови шуканої залежності, при якому за отриманими інтервальними оцінками дисперсії шуму створюється таблична функція, що інтерполює. Дана таблична функція формується на основі робастної локальної лінійної апроксимації інтервальних оцінок дисперсії. Використання робастних методів дозволяє додатково знизити вплив викидів (грубих помилок в інтервальних оцінках дисперсії), у той час як локальність апроксимації забезпечує повторення складного ходу кривої, що описує залежність шуму від інтенсивності. Таким чином, отримана таблична функція кожної інтенсивності вихідного зображення ставить у відповідність оцінку дисперсії шуму. Крапками входу в таблицю можуть служити, наприклад, інтенсивності зображення, що оцінюється. На практиці, у випадку параметричної оцінки шуму, може бути використаний підхід, при якому застосовується перетворення, що стабілізує дисперсію шуму вихідного зображення [2, 3, 9, 11]. При цьому проблема фільтрації шуму, що залежить від сигналу, зводиться до задачі пригнічення адитивного незалежного від сигналу шуму постійної (заданої) дисперсії. У цьому винаході здійснюють непараметричну оцінку шуму, тому на етапі побудови карти шуму пропонують використовувати наступний підхід: на основі зображення, що оцінюється, й таблиці, що інтерполює, побудувати карту шуму - зображення, кожний піксель якого оцінює середньоквадратичне відхилення шуму у відповідному пікселі вихідного зображення. Карта шуму дає попіксельну оцінку шуму з достатньою для практичного застосування точністю. Використання як можна більше точної інформації про рівень шуму на зображенні дозволяє істотно поліпшити якість алгоритмів шумозаглушення. 3 UA 104164 C2 Спосіб, що заявляється, здійснюють у такий спосіб. На етапі побудови зображення, що оцінюється, й одержання зображення шуму вихідне зображення І(х, у) фільтрують наступним низькочастотним лінійним біноміальним фільтром розміру 3×3: 5 10 15 T Н1=[1 2 1]/4, H  H1  H1 . При цьому одержують згладжене зображення Іе(х, у)=І*Н. Далі обчислюють зображення шуму Ne(х, у)=І(х, у)-Іе(х, у). На етапі відкидання границь, використовуючи зображення шуму Ne(x, y), формують два зображення позитивних і негативних змін: 1, Ne x, y   0  Ne x, y    0, Ne x, y   0 , 1, Ne x, y   0  Ne x, y    0, Ne x, y   0 . Для виділення великих областей, що відповідають границям, дані бінарні зображення послідовно піддають морфологічним операціям ерозії й дилатації, використовуючи маски розміру 2×2:    , erode N x, y .    Be x, y  dilate22 erode22 Ne x, y    Be x, y  dilate22  22  e Потім ці зображення поєднують для одержання карти границь вихідного зображення   Ex, y   Be x, y   Be x, y  . 20 25 На етапі обчислення інтервальних оцінок дисперсії шуму обчислюють мінімальну І min й максимальну Іmах інтенсивності зображення Iе(х, у), вибирають крок hM i розбивають діапазон інтенсивностей на інтервали Mi із кроком hM (hM=32). Для кожного пікселя зображення Iе(х, у) знаходять інтервал, якому належить значення цього пікселя, а відповідне значення пікселя 2 зображення Ne(x, y) використовують для обчислення оцінки дисперсії шуму σ (і) в даному інтервалі Мі, при цьому виключають значення пікселів, для яких значення карти границь Е(х, у)=1. При обчисленні інтервальної оцінки дисперсії шуму застосовують формулу незміщеної оцінки дисперсії: 2 i  де Nie 2 1  ni 1 i    N j  Nie    e ni  1  j1       j , - значення пікселя зображення шуму Ne(x, y) з інтервалу Mi, ni - загальна кількість i 30 35 накопичених значень пікселів зображення шуму в інтервалі Мi, Ne - середнє значень шумових пікселів в інтервалі Mi. 2 На етапі уточнення інтервальних оцінок дисперсії величини σ (і) уточнюють таким чином, що для кожного інтервалу Mi ітеративно виключають значення шумових пікселів, що перевищують за абсолютною величиною поріг, що дорівнює трьом стандартним відхиленням шуму, з наступним перерахуванням оцінки дисперсії шуму в даному інтервалі:      Nie k  1  Ne x, y  | Ne x, y   Nie k   Ne x, y   3i , де Ne k  1 - набір значень шумових пікселів в інтервалі Мi на ітерації k. У подальших розрахунках використовуються лише ті інтервали, у яких накопичена достатня кількість значень шумових пікселів (наприклад, не менше 500). Крім того, з розгляду виключаються ті інтервали, i 40 45 i середнє значення Ne  j в яких істотно відрізняється від нуля (відхиляється від нуля більш ніж на половину кроку сітки інтервалів hM), оскільки в таких інтервалах з високою ймовірністю домінують залишкові значення нешумових пікселів. На етапі оцінки залежності дисперсії шуму від інтенсивності за отриманими інтервальними оцінками дисперсії шуму створюють табличну функцію, що інтерполює. Дану табличну функцію формують на основі робастної локальної лінійної апроксимації інтервальних оцінок дисперсії. Для цього на сітці інтервалів Mi вибирають крок hI i радіус rI апроксимації (які можна зробити залежними від кількості крапок nі в інтервалах Мі). Значення отриманої на попередньому етапі 4 UA 104164 C2 табличної функції дисперсії шуму від інтенсивності сигналу апроксимують за наступним правилом: Є2 k  hI   a  mk  b , 5 де k - номер вузла апроксимації; mk - центр інтервалу Mi(k·hI); параметри а, b обчислюють на основі мінімуму суми абсолютних відхилень [12] khI rI  jkhI rI 10 15 20 25 30 35 40 45 Є2  j  a  m j  b  min a,b .   2 Отриману в результаті такого процесу таблицю значень MI k  hI , Є k  hI  інтерполюють на всьому діапазоні інтенсивностей [Imin, Imax], тобто одержують таблицю, що інтерполює, шуканої залежності дисперсії шуму від інтенсивності сигналу. На етапі побудови карти шуму на основі зображення, що оцінюється, й таблиці, що інтерполює, будують карту шуму - зображення, кожний піксель якого оцінює дисперсію шуму у відповідному пікселі вихідного зображення. На фіг. 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 C., Szeliski R., Kang S.B., Zitnick C.L., Freeman W.T. Automaticestimationan dremovalofnoisefromasingleimage. IEEETransactionsonPatternAnalysisand Machinelntelligence30(2), 299 {314(2008). 8. Salmeri M., Mencattini A., Rabottino G., Lojacono R. SignaldependentNoiseCharacterizationforMammographiclmagesDenoising. ІМЕКО ТС4 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. Яні Б. Цифрова обробка зображень. М., Техносфера, 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. ФОРМУЛА ВИНАХОДУ 50 55 Спосіб оцінки шуму цифрових рентгенограм, що включає одержання вихідного зображення, одержання зображення, що оцінюється, за допомогою низькочастотної фільтрації вихідного зображення, побудову зображення шуму як різниці між вихідним і зображенням, що оцінюється, видалення значень пікселів зображення шуму, що відповідають різким змінам у вихідному зображенні, розбивку діапазону інтенсивності зображення, що оцінюється, на інтервали, при цьому кожний піксель зображення, що оцінюється, належить певному інтервалу, накопичення для кожного інтервалу значень пікселів зображення шуму, що відповідають пікселям 5 UA 104164 C2 5 10 зображення, що оцінюється, обчислення інтервальних оцінок дисперсії шуму за накопиченими у даному інтервалі значеннями пікселів шумового зображення, уточнення інтервальних оцінок дисперсії шляхом відкидання значень шумових пікселів за критерієм 3-сигма, який відрізняється тим, що здійснюють видалення значень пікселів зображення шуму, що відповідають різким змінам у вихідному зображенні, за допомогою морфологічного виділення значень пікселів зображення шуму, що відповідають границям на вихідному зображенні, виконують робастну локальну лінійну апроксимацію інтервальних оцінок дисперсії шуму, у результаті якої одержують табличну функцію, що описує залежність шуму від інтенсивності сигналу, обчислюють на основі зображення, що оцінюється, й побудованої табличної функції, що описує залежність шуму від інтенсивності сигналу, карту шуму у вигляді попіксельної оцінки дисперсії шуму вихідного цифрового зображення. 6 UA 104164 C2 7 UA 104164 C2 Комп’ютерна верстка Л. Бурлак Державна служба інтелектуальної власності України, вул. Урицького, 45, м. Київ, МСП, 03680, Україна ДП “Український інститут промислової власності”, вул. Глазунова, 1, м. Київ – 42, 01601 8