Узагальнений метод швидкого обчислення перетворення Радона із використанням властивостей симетрії схеми обчислення Текст научной статьи по специальности «Математика»

УДК 004.421.2:517.44

УЗАГАЛЬНЕНИЙ МЕТОД ШВИДКОГО ОБЧИСЛЕННЯ ПЕРЕТВОРЕННЯ РАДОНА ІЗ ВИКОРИСТАННЯМ ВЛАСТИВОСТЕЙ СИМЕТРІЇ СХЕМИ ОБЧИСЛЕННЯ

КОБАСЯР М.І., РУСИН Б.П.______________________

Описується метод швидкого обчислення перетворення Радона. Він використовує властивості симетрії схеми обчислення, що дозволяє обчислювати меншу кількість характеристик. Проведені дослідження показали ефективність запропонованого методу, що дозволив суттєво зменшити час обчислення перетворення Радона.

1. Вступ

Перетворення Радона (ПР) [1], як і його частковий випадок перетворення Хафа (ПХ) [2], є відомими засобами в обробці зображень, томографії, астрономії, мікроскопії тощо [4, 6].

Лінійне ПР можна визначити у загальній формі таким чином [1, 4]:

+<Х X

g ^ 0 5 Sl£ 2,) = J J g(x,y)5^ 0-^1x-S 2y)dxdy,(1)

де g(x, y) - функція густини об’ єкта,

5(£, 0 — ^x 2y) - дельта-функція Дірака.

У виразі (1) лінія описується трьома степенями свободи. Із багатьох, наведених у літературі, найбільш

поширеними є (p, т) (тангенціальна, slant stacking) та (р, 0) (нормальна) форми ПР або, відповідно, ТПР та НПР.

Звичайно дослідники використовують не ТПР, а НПР. Це пояснюється тим, що на відміну від ТПР для НПР відсутні будь-які обмеження, зокрема на орієнтацію відрізків ліній.

НПР використовує запис рівняння лінії у нормальній формі [7] та визначається як інтеграл від функції

густини зображення g(x, y) вздовж прямої інтегрування 5:

+«>

g(p, Є) = J g(x,y)ds =

= J g(p cos 0 - s • sin 0, p sin 0 + s • cos 0)ds =

+<x +<x

= J Jg(x, y)5(p- xcos 0- y • sin 0)dxdy (2)

—да —да

де p - найкоротша відстань від початку системи координат до прямої інтегрування s ; 0 - кут між прямою інтегрування s та віссю абсцис.

Дискретну версію виразу (2) називають променевою сумою або променевою проекцією. Повна система променевих сум для даного кута 0 називається проекцією. Сукупність усіх проекцій для 0 є [0; тс] формує область параметрів (ОП).

З історичної точки зору не зрозуміло, що стало причиною появи ПХ. Вірогідно, що це була нова ідея, яка з ’явилася на півстоліття пізніше від ПР. Слід відзначити про плутанину в літературі стосовно зв’язку між ПХ та ПР. Відомо, що ці два різних підходи виділення параметрів дають, в основному, подібний тип дискретного перетворення.

Якщо зображення є розрідженим, тобто бінарним із невеликою кількістю ненульових пікселів, тоді більшість машинного часу витрачається на сумування нулів, що не роблять внеску в ОП. В одному із найбільш цитованих патентів по розпізнаванню зображень [2] Хаф (P.V.C. Hough) запропонував спосіб зменшення часу обчислення. Спочатку ідея була сформульована для відрізків ліній, що описувалися нахилом і зміщенням подібно до (p, т) форми ПР. Проте багато авторів застосували цю ідею для лінії, описаної в (р, 0) формі [3].

Нормальну форму ПХ (НПХ) можна визначити із використанням нормальних параметрів прямої (р, 0) як

+<Х +<Х

g(P, Є) =J J g(x* 5y* )S(-)dx*dy* , (3)

де 5(-) = 5(р - x* cos 0 + y* sin 0). Для пояснення цих перетворень використано той факт, що ПХ відображає окремі пікселі з області зображення у криву (синусоїду) в ОП, а ПР перетворює криву з області зображення в один піксель в ОП.

Мета дослідження: створення нових ефективних алгоритмів обчислення ПР.

2. Порівняння обчислювальної складності перетворень Радона та Хафа

Згідно із [5], обчислювальна складність НПР при її реалізації за допомогою інтерполяцій найближчим сусідом, лінійної та sinc-інтерполяцій визначається відповідно як

Onn = O(RTN) * O(N3),

0 Linear = O(RTN) * O(N3), (4)

Osinc = O(RTN2) * O(N4),

де O( ) - функція порядку обчислень, що нехтує елементами меншого порядку; N - ширина/висота зображення. Для N = 100 обчислювальна складність sinc-інтерполяції у 100 раз вища за інші два методи, оскільки sinc-функція обчислюється для кожного значення у виразі (2). Лінійна інтерполяція є повільнішою за інтерполяцію найближчим сусідом [5].

104

РИ, 2007, № 4

Обчислювальна складність ПХ рівна

Oht = O(N) * O((NN)rN), (5)

де індекс r вказує, що може бути перетворено скінченну кількість пікселів. Якщо на зображенні є тільки один ненульовий піксель, тоді у виразі (5) (NN)r = 1. ПХ буде швидшим за ПР. Проте реальна величина

(NN) лежить між N і NN (залежно від структури зображення), що збільшує час обчислень за допомогою ПХ. Це також свідчить, що найкраще застосування ПХ буде мати у випадку обробки бінарних зображень.

З метою порівняння існуючих методів обчислення ПР та ПХ було проведено дослідження залежності часу обчислення ОП від розміру зображення для різних інтерполяцій НПР і НПХ, описаних у [5]. Результати дослідження наведено на рис.1. Графіки 1-4 відповідають різним методам оптимізації НПХ для інтерполяції найближчим сусідом, а графіки 5, 6 - відповідно інтерполяції найближчим сусідом і лінійній інтерполяції НПР.

2. Для методів оптимізації НПХ 3 і 4 час обчислення є приблизно однаковим. Для великої кількості ненуль-ових пікселів зображення ці методи вимагають набагато менше часу порівняно із методами 1 і 2.

3. Для усіх методів інтерполяції НПР час обчислення ОП є приблизно однаковим.

4. НПХ є швидшим від НПР для невеликої кількості ненульових пікселів на зображенні (до 8-10%) і для відносно малого розміру зображень (до 350 х 350 пікселів). Для більших зображень ПР вимагає значно менше часу обчислень у порівнянні із ПХ. Це є новим і досить цікавим фактом.

3. Характеристики системної матриці при обчисленні перетворення Радона

При формуванні ОП за допомогою ПР і ПХ піксель розглядається як точкове джерело. Згідно із виразами (2) та (3) кожний ненульовий піксель на зображенні робить однаковий вклад у проекційну суму. Це не зовсім коректно. Скажімо, у випадку застосування ОП для детектування відрізків ліній і моделювання лінії із параметрами (р*, 0*) за допомогою дельта-функції формується пік в координатах лінії та викиди у всіх інших випадках.

Розглянемо формування ОП у випадку ПР детальніше. Оскільки ПР є лінійним перетворенням по відношенню до функції g(x, y) і ми працюємо із дискретною версією ПР, тому його можна подати в матричному представленні замість інтегрального:

g(p, Є) = ЭТ g(x, y) b = W g

(6)

де W IxJ - системна матриця (СМ), яка містить

вагові коефіцієнти та у між j-ми пікселями зображення та кожною орієнтацією і прямої інтегрування s ; I -вимірний вектор b (I = RT) відповідає ОП; J -вимірний вектор g (J = MN) відповідає зображенню; R та T - відповідно, кількість вибірок по р і 0; M та N - відповідно, ширина і висота зображення.

На основі наведеного вище вираз (3), а саме проекційну суму для довільних значень р і 0, можна записати як суму добутків елементів зображення g j і СМ

Рис. 1. Залежність часу обчислень ОП від розміру зображення для різних інтерполяцій НПХ (1-4) і НПР (56). Масштаб 0..700 (а) і 75..200 (б) пікселів

N—1

g(Px ,61) = Z gjra i,j

j=0

(7)

На основі одержаних графічних залежностей на рис.1 можна зробити такі висновки щодо обчислення ОП:

1. Ініціалізація векторів і масивів НПХ вимагає дуже багато часу навіть у випадку невеликої кількості нену-льових пікселів зображення.

де gi(p, 0) = g(px, 01) - проекційна сума; та у - ваговий коефіцієнт між j -м елементом зображення та кожною і -ю орієнтацією прямої інтегрування s .

РИ, 2007, № 4

105

Звичайно визначення елементів СМ та у проводять, використовуючи один із двох варіантів обчислень:

1) ваговий, коли обчислення та у іде точно; 2) нева-говий, коли та у у виразі (7) прирівнюються до одиниці при перетині прямої інтегрування s елемента зображення та рівні нулю у всіх інших випадках. Звичайно використовують неваговий підхід, який вимагає додаткової апроксимації вкладу та у в ОП, проте має суттєві переваги у швидкості обчислень, скажімо, у випадку застосування для детектування відрізків ліній [2, 3].

СМ складається із RTMN « N х N елементів і є дуже великою. Вона не є структурованою [8]. Матриця є розрідженою і має близько O(N) ненульових елементів для кожної прямої інтегрування. Легко визначити, що кількість ненульових величин для лінії є

менше за d(N - 1) +1, де d - розмірність об’єкта (для зображень d = 2).

Наприклад, зображення розміром 100 х 100 пікселів

вимагає обчислення і збереження ~108 (~3.6 • 106 ненульових) елементів СМ, а зображення розміром

512х512 пікселів - уже ~6.8-1010 (~5.3• 108 ненульових) елементів СМ. Це вимагає значних ресурсів пам’яті навіть для сучасних ПЕОМ.

Сукупність елементів СМ довільної проекційної суми можна розглядати як дискретизований за рівнем сигнал. Згідно із [9] груба оцінка перевищення сигналу над шумом квантування визначається таким виразом:

Ddb = 10 lg 22к = 10 • 2кlg 2 « 6к, (8)

де к - кількість двійкових розрядів квантування, що пов’язана із кількістю рівнів квантуванняL виразом

L = 2к .

Згідно із [9] величина Ddb у виразі (8) рівна 6 дБ на один розряд, якщо Kpf = 2.5 (пік-фактор сигналу, тобто відношення його максимального значення до середньоквадратичного), а якщо Kpf = 3, вона зменшується до 5.5 дБ на один розряд.

Очевидно, що для невагового способу обчислення елементів СМ к = 1 і відповідне відношення сигналшум складає 5.5 дБ. Для вагового способу можна отримати довільне значення D db , що обмежується лише значенням мантиси представлення та у . У літературі [9] наводять дані, що достатньо задати к = 9 (відповідає Ddb = 54 дБ).

Іншою важливою характеристикою шуму квантування є його спектральна характеристика. Для невагово-го вибору елементів СМ спектральна характеристика

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

На рис .2 наведено результати дослідження залежності часу обчислення ОП за допомогою інтерполяції ПР найближчим сусідом та лінійної інтерполяції для невагового та вагового способів обчислення ПР. Як бачимо, ці способи мають лінійну залежність і різниця у часі порівняно незначна.

Рис.2. Залежність часу обчислень ОП від розміру зображення. Інтерполяція найближчим сусідом та лінійна інтерполяція відповідно для невагового (1, 2) та вагового (3,4) способів обчислення ПР

4. Властивості симетрії схеми обчислення перетворення Радона

Оскільки метод швидкого обчислення ПР використовує властивості симетрії схеми інтегрування, спочатку доведемо існування цих властивостей.

Для квадратного зображення розміром N х N , причому N - парне, центр системи координат O збігається із центром зображення (рис.3). Проекції вибираються із довільним рівномірним кроком по куту Д0, а зміщення проекційних сум у цих проекціях - із довільним рівномірним кроком Ар .

Центральна симетрія. Розглянемо центральну симетрію елементів СМ для деякого такого кута а, причому а є (0;450). Відповідно проекційні суми обчислюються із однаковим по модулю додатнім та від’ємним зміщенням відносно початку системи координат. Згідно із рис.3,а р = |-р| = OM = |- OM'|, де рє (0;NV2) .

Для додатного зміщення пряма, вздовж якої проводиться інтегрування, утворює прямокутний АЛБС, а для від’ємного - прямокутний ДА' Б' С'.

Слід довести, що A ABC = АА' Б' С' і відповідні довжини відрізків в елементах зображення прямої інтегрування рівні, а координати відрізків прямої інтегрування в елементах зображення визначаються через однакові матриці інцидентності по x та у координатах.

106

РИ, 2007, № 4

Рис. 3. Ілюстрація до доведень властивостей центральної (а) та поворотної і поворотно-дзеркальної (б) симетрій СМ

1. AABC є прямокутним із ZABC = 90 0 . Знайдемо його інші кути.

ZE' OK' = а . Інтегрування проводиться для рівнові-ддалених паралельних зміщень, тобто

KK'|| AC|| C'A', звідки

ZOLS = ZE' OK' = ZE' L'K' = a .

Відрізки LO і AP є перпендикулярами до осі Oy, тобто LO || AP, звідки ZPAS = a . APAS і AOLS є прямокутними. Оскільки ZASP = ZLSO (спільний кут) та ZPAS = ZOLS = a, тоді APAS і AOLS є подібними за теоремою подібності про три кути, тобто APAS<xAOLS.

Оскільки ZPAS + ZCAB = 900 , тоді

ZCAB = 900 —а . Сума кутів довільного трикутника рівна 1800, звідки РИ, 2007, № 4

ZBCA = 1800 - ZABC - ZCAB =

= 1800 -900 -(900 -а) =а .

2. Аналогічно AA'B'C' є прямокутним із ZA' B' C' = 900 . Знайдемо його інші кути.

Оскільки KK'||C'A', тоді ZE'OK' =ZE'L'A' =a . Тоді для прямокутного AA' E' L'

ZL'A'E' = 1800 - ZA' E' L'-ZE' L' A', звідки ZL'A'E' = 900 —a .

Оскільки E'L'||B'C', тоді прямокутні AA'B'C' і AA'E'L' із спільним кутом ZL'A'E = 900 -а є подібними (за згадуваною вище теоремою подібності) і ZE'L'A' = ZB'C'A' =a .

3. Отже, AABGxAA'B'C' за теоремою подібності прямокутних трикутників за гострим кутом (в AS' OL' ZOL' S' = ZEL' A' = a як взаємовертикальні кути). Оскільки OM = |OM'| і зображення квадратне, тоді AS' OL' = ASOL та відповідно AABC = AA' B' C'.

Задамо деяку точку V, що належить відрізку AC. Точку V задамо так, щоб вона лежала між двома суміжними колонками зображення. Провівши перпендикуляр до осі Oy через точку V, позначимо точкою U перетин із лівою границею зображення. Аналогічно провівши перпендикуляр до осі Oy через точку V' , позначимо точкою U' перетин із правою границею зображення. Розглянемо AAUV і AA'U'V'.

ZAVU = ZA'V'U' виходячи із подібності прямокутних трикутників за трьома кутами, тобто

AAUV<xAA'U'V'. За побудовою UV = U'V', тоді AAUV = AA' U' V'. Останній вираз є коректним за довільного паралельного переносу точки V до осі Oy по AC.

Таким же чином доводимо, що елементам СМ на рівновіддалених від центра зображення прямих інтегрування із

ає (45°;90°), ає (90°;135°) і ає (1350;1800)

по відношенню одним до одного притаманна властивість центральної симетрії.

Як і слід було довести, відповідні довжини відрізків в елементах зображення прямої інтегрування рівні і однакові матриці інцидентності відповідних координат.

Поворотна і поворотно-дзеркальна симетрії. Слід довести, що відповідні довжини відрізків в елементах зображення прямої інтегрування рівні для деякого

кута а є (0;45 0) і будуть рівні для кутів через ~ та

будуть однакові матриці інцидентності відповідних координат.

Знову відповідні проекційні суми обчислюються із однаковим по модулю додатнім та від’ ємним зміщен-

107

ням відносно початку координат. Згідно із рис.3,б р = |-р| = OM = |-OM'|, де рє (0;N/2).

Для деякого додатного зміщення пряма із кутом а , вздовж якої проводиться інтегрування, утворює прямокутний AABC (рис.3,б). Інша пряма із таким же додатнім зміщенням відносно початку координат із %

кутом — ~ а утворює прямокутний трикутник

AA' B' C'. Слід довести, що AABC = AA' B' C' і відповідні довжини відрізків прямої інтегрування рівні, а координати відрізків визначаються через однакові матриці інцидентності по х та у координатах.

1. AABC є прямокутним, що було доведено вище.

2. Аналогічно AA'B'C' є прямокутним із ZC' B'A' = 90 . Знайдемо його інші кути.

Оскільки KiKi'||C'A', тоді

ZE'OK1' = ZE'L'C' = 900 -а . Для прямокутного AC' E' L' ZL' C'E' = 1800 - ZC' E' L'-ZE' L' C',

звідки ZL'C'E' =а.

Оскільки E'L'||B'A', тоді прямокутні AA'B'C' <xAL'E'K' за теоремою про три однакові кути. Отже, ZE'L'K' =ZB'A'C' = 900 —а .

3. AABGxAA'B'C' за теоремою подібності прямокутних трикутників за гострим кутом. Оскільки OM = |OM'| і зображення квадратне, тоді AABC = AA'B'C'.

Аналогічно до доведення про центральну симетрію покажемо, що AAUV = AA'U'V'.

Точка V' є симетричною до точки V (властивість поворотно-дзеркальної симетрії). Якщо б ми використали додатне зміщення прямих інтегрування для даного кута a (A'C' є від’ємним зміщенням для деякого а є (450;900)), тоді відповідна точка V' була б симетричною до точки V згідно із властивістю поворотної симетрії.

Таким же чином доводиться, що елементам СМ на рівновіддалених від центра зображення прямих інтегрування із а є (90°;135°) і а є (1350 ;1800) притаманна властивість поворотної або поворотно-дзеркальної симетрії відповідним елементам СМ на рівно-віддалених від центра зображення прямих інтегрування із ає (00;450).

Як і слід було довести, відповідні довжини відрізків в елементах зображення прямої інтегрування рівні і однакові матриці інцидентності відповідних координат.

5. Метод швидкого обчислення вагового перетворення Радона

Було запропоновано метод швидкого обчислення вагового перетворення Радона (ВПР). Він визначає елементи СМ як довжини прямої інтегрування в елемен-

тах зображення. Для розробленого методу характерні такі особливості [10]:

1) Врахування та попереднє обчислення допоміжних величин для визначення елементів СМ

З точки зору мінімізації часових затрат доцільно здійснювати попереднє обчислення допоміжних величин елементів СМ, зумовлених схемою обчислення ПР. Ці допоміжні елементи використовують під час обчисленні ОП, проте їх пряме збереження потребує значних затр ат пам ’ яті.

Приклади даних про кількість необхідних елементів СМ залежно від розміру зображення наведено вище. Виходячи із цих даних, робимо висновок про неможливість прямого збереження усіх елементів, а також допоміжних величин, призначених для їх обчислення.

2) Використання характерних для схеми обчислення властивостей симетрії

Ретельний аналіз геометрії інтегрування показує, що за довільного вибору кроків паралельного і кутового переміщення під час інтегрування довжини відрізків прямих інтегрування s , обмежених границями елементарних елементів зображення, координати початків та кінців цих відрізків, довжини відрізків прямих інтегрування s , обмежених областю зображення, та координати початків та кінців цих відрізків підпорядковуються властивостям центральної, поворотної і поворотно-дзеркальної симетрій. Повторимося, що ці симетрії наявні для зображень N х N , причому N -парне і центр системи координат збігається із центром зображення. Це дає змогу обчислювати та зберігати значно меншу кількість характеристик, пов’язаних з геометрією інтегрування. Звертання до значень цих хар актеристик здійснюється за відпо відними індексами, нескладні формули розрахунку яких визначаються типом симетрії (див.рис.3) та будуть подані нижче.

Прямі інтегрування одного ракурсу, зміщені на однакові відстані відносно центра зображення, формують елементи СМ. Для цих елементів притаманна центральна симетрія, яка дає можливість обчислювати геометричні характеристики лише для невід’ємних зміщень. Це дозволяє зменшити кількість необхідних характеристик у 2 рази (див.рис.3,а). Використання властивостей поворотної і поворотно-дзеркальної симетрій для різноракурсних прямих інтегрування s із 9є [00 ;1800) дозволяє обчислювати геометричні характеристики прямих інтегрування s лише із

0 є [0 ;45 ] , що зменшує кількість необхідних характеристик ще у 4 рази (див.рис.3,б). У роботі [11] було описано і використано останню властивість симетрії для підвищення швидкодії обчислення ПР (але, на жаль, без наведення часових характеристик і порівнянь). У статті [10], яку було опубліковано на кілька років раніше згаданого джерела [11], описано

1 використано усі наведені вище властивості симетрії, а дана публікація є її узагальненням.

108

РИ, 2007, № 4

Дослідження показали (рис.4), що запропонований спосіб вимагає у 8 разів меншу кількість необхідних характеристик для обчислення елементів СМ та, відповідно, ОП.

Рис.4. Залежність кількості ненульових елементів СМ від розміру зображення для запропонованого (1) та існуючого (2) способів

3) Різницеве подання матриці інцидентності координат елементів зображення та прямих інтегрування s:

Під час реалізації обчислення ОП важливо за заданими параметрами р і 0 прямої інтегрування s визначити кількість та послідовність елементів зображення, через які проходить ця пряма інтегрування s .

Що б зменшити час обчислень та нео бхідну оперативну пам’ять, обрано таку схему подання інформації про інцидентність елементів зображення та прямих інтегрування s;

- для кожної з прямих інтегрування s (9є [00 ;45 0 ] , р > 0 - невід’ємні відстані від центру зображення) визначаємо координати елементів зображення входження прямих інтегрування s в область зображення;

Yj = Yj-i +ЛУ^ (10)

де Ayj = <

(-i)uf(yj,еt,Рх)Ay, 0<еt <-;

4

(-1)u+1f(xj,Єt,Рх)Ax, 7<01 <-2;

4 2

% 3

(-1)uf(xj, 01, Px )Ax, -<01 < - я;

24

3

(_1)u+1f(yj,et,px)Ay-%<Qt <%,

4

тут u = 0 - індекс для додатних зміщень рх; u = 1 -для від’ємних зміщень рх ; f(xj,01,рх), f(Уj, 61, Рх) - бінарні величини, які свідчать про наявність (або відстутність) приросту відповідно х або у координати елемента зображення, при його перетині прямою інтегрування цього елемента зображення; Ax j, Ay j - приріст відповідної координати елемента зображення. Приймемо розмір елемента зображення рівним Ax = Ay, тобто Ax х Ax .

Визначені аналітично координати точок входження прямих інтегрування подано у таблиці, де IYst - у координата точки входження прямої інтегрування у зображення для 0 є [0;450]. 1-й та 2-й октанти відповідають 0 є [0;450], 3-й і 4-й - 0 є (45°,90°), 5-й і 6-ий - 0є [900,1350], 7-й і 8-й - 0є (1350,1800). Непарні значення октантів відповідають р > 0 , а парні - р < 0 . IYst визначається із рівняння прямої для конкретного значення р і 0.

Координати точок входження прямих інтегрування у зображення

- подальшу траєкторію проходження прямих інтегрування s через область зображення подаємо значеннями приростів координат під час переходу s із одного елемента зображення до наступного.

Згідно із рис.3 координати елемента зображення розміром Ax х Ay , який перетинається ідеальною прямою інтегрування, визначаються такими виразами:

xj = xj_1 + Axj, (9)

(-1)uf(xj,01,Px)Ax, 0<01 <T;

4

де

Axj =<

(-1)u+1f(yj,01,Px)Ay, 7<01 <7;

42

3

(-1)u+1f(yj,01,Px)Ay,-<01 <7я;

24

3

(-1)uf(xj,01,Px)Ax, — n<et <n,

4

r > 0 r < 0

№ ок. Координата № ок. Координата

х0 y0 х0 y0

1 0 IYst 2 N-1 N-IYst-1

3 N-IYst-1 N-1 4 IYst 0

5 N-IYst-1 0 6 IYst N-1

7 0 N-IYst-1 8 N-1 IYst

4) Організація структури даних, яка орієнтована на мінімізацію часу та пам ’яті

На основі міркувань, найсуттєвіші з яких наведено вище, розроблено структуру даних, яка передбачає обчислення і збереження приростів координат при проходженні прямими інтегрування s елементів зображення; нормованих довжин відрізків, що утворю-

РИ, 2007, № 4

109

ються при перетині прямих інтегрування s та границь елементів зображення.

Зрозуміло, що недоцільно зберігати повну мантису значень елементів СМ, які лежать в діапазоні

(0..Дхл/2) для елемента зображення Дх хДх. Проведені дослідження розв’язку зворотної задачі (ре-конструювання зображень за їх проекціями ітерацій-ним методом ART, відомим як алгебраїчний метод корекції або попроменева корекція [8]) засвідчили доцільність зберігання 3 -4 знаків мантиси для Аг = 1. Тоді елементи СМ представляють як 2- байтове ціле, у той час як збереження повної мантиси елементів вимагає 4-10 байт залежно від реалізації. Згідно із виразом (8) використання вагового підходу обчислення елементів СМ із такою мантисою дозволяє досягнути співвідношення сигнал-шум квантування рівного 60 дБ.

5) Особливості реалізації обчислення ОП

Обчислення ОП проводять послідовно для кожного октанту. Підбір кроків дискретизації параметрів ОП має важливе значення, оскільки дозволяє мінімізувати час її обчислень без втрати інформації або навпаки без обчислення її надлишку. Розв’язуючи зворотну задачу ПР за допомогою ART, визначено, що кількість проекцій можна зменшити на 7.5%, а кількість проекційних сум у цих проекціях на 10% порівняно із теоретичними оцінками [12].

6) Застосування властивостей симетрії для зображень довільного розміру

Властивості симетрії схеми обчислення ПР, тобто його ОП, можна застосувати для зображень довільного розміру. Потрібно зображення представити у вигляді N х N , причому N - парне. Для цього доповнюємо зображення рядками (стовпцями) із нульовими значеннями елементів зображення, якщо, відповідно, ширина (висота) зображення є більшою. Якщо більша сторона зображення має непарну кількість елементів, обидві сторони зображення доповнюємо нульовими значеннями елементів зображення. Для того щоб не обчислювати ПР доповнених рядків (стовпців) зображення, можна для усіх прямих інтегрування визначити точки входження та виходу із області інтегрування (зображення).

Отриману ОП легко привести до дійсного значення. Для цього потрібно використати властивість зв’язку ПР початкового g(x,у) і зміщеного g(x-X0,y-У0) зображень. При зміщенні зображення змінюється лише р координата, тобто h(p, 0) = g(p- х0 cos0- у0 sin0,0)

[7]. Оскільки зміщення (Х0У0) відоме, ОП легко перегрупувати.

Усі додаткові обчислення у будь-якому випадку будуть вимагати значно менше часу порівняно із визначенням ПР без використання властивостей симетрії.

Проведені дослідження показали (рис.5), що використання властивостей симетрії у випадку двовимірної

інтерполяції ПР зменшує час обчислень у 2.8-5.7 разів, а у випадку одновимірної - в 2.4-3.6 разів.

Запропонований метод швидкого обчислення ВПР із використанням властивостей симетрії передбачає фактично одночасне обчислення 8 інтегралів ПР. Його реалізація в паралельних обчислювальних системах потребує мінімальних змін алгоритму. Слід очікувати підвищення швидкодії обчислення невагового та вагового ПР приблизно у 7-8 разів порівняно із класичними способами. Т еоретичного значення підвищення швидкодії у 8 разів не можна досягнути через те, що запропонований метод обчислення із використанням властивостей симетрії передбачає виконання кількох додаткових операцій (вирази (9) і (10)) для визначення координат поточних елементів зо браження і відповідних елементів СМ.

Рис.5. Залежність часу обчислень ПР від розміру зображення при лінійній двовимірній (1,2) та одно-вимірній (3,4) інтерполяції координат відповідно для існуючого (1,3) та запропонованого (2,4) способів

Висновки

Розв’язано нову наукову задачу розробки методу швидкого обчислення перетворення Радона, що використовує властивість симетрії схеми обчислення цього перетворення. Зокрема отримано такі результати:

1. Показано, що існуючі методи обчислення перетворення Радона є порівняно повільними. Тому слід усунути їх основний недолік - великі обчислювальні затрати.

2. Запропоновано метод швидкого обчислення перетворення Радона, що використовує властивості симетрії схеми обчислення. Це дозволило обчислювати і зберігати лише восьму частину характеристик, необхідних для обчислення повної системної матриці. Використання вагового підходу дало можливість збільшити відношення сигнал-шум дискретизації елементів із 5.5 до 60 дБ порівняно із неваговим підходом. Швидкість обчислення перетворення Радона завдяки використанню властивостей симетрії збільшено у 2.4-3.6 разів для одновимірної інтерполяції координат та у 2.8-5.7 разів для двовимірної.

3. Для реалізації класичних і вагових перетворень Радона та Хафа розроблено алгоритми і програмне забезпечення, які знайшли застосування для реконст-

110

РИ, 2007, № 4

рукції зображень за їх проекціями, погодження різно-ракурсних зображень об’єктів.

Література: 1. Radon J. Uber die Bestimmung von Funktionen durch ihre Integral-werte langs gewisser Mannigfaltigkeiten. Berichte Sachsische Akademie der Wissenschaften. Leipzig, Math.-Phis. Kl. 1917. Vol.69. P.262267. 2. Hough P. V. C. Methods and Means for Recognizing Complex Patterns. U.S. Patent 3069654. 1962. 3. Leavers V. Which Hough transform? A survey of Hough transform methods // CVGIP - Image Understanding. 1993. Vol.58, №2. P.250-264. 4. Наттерер Ф. Математические аспекты компьютерной томографии. М.: Мир. 1990. 386 с. 5. ToftP. The Radon Transform - Theory and Implementation: Ph.D Thesis. Department of Mathematical Modelling, Technical University of Denmark. 1996. 230 р. 6. Steeghs P. Wigner-Radon representations for 3-d seismic data analysis. Delft University of Technology, Delft, the Netherlands. 2000. 4 р.

7. Deans S. R. The Radon transform and some of its applications. Malabar, Florida: Kriegel Publishing Company. 1993. 295 р. 8. Dijke M. C. A. Iterative methods in image reconstruction: Phd thesis. University of Utrecht. 1992. 9. Го-норовский И. С. Радиотехнические цепи и сигналы: Учебник для вузов. М.: Радио и связь. 1986. 512 с. 10. Бєляєв В. П., Кобасяр М. І., Русин Б. П. Використання властивостей симетрії схеми зйому під час реконструкції зображень за їх проекціями // Відбір і обробка інформації. 1999.

Вип. 13(89). С. 143-148. 11. LeeH.-J., AhnH.-J., Song J.-H, Park R. H. Hough transform for line and plane detection based on the conjugate formulation // Machine Vision Applications in Industrical Inspection IX (Proceedings of SPIE Volume 4301). San Jose, CA. 2001. P.244-252. 12. Кобасяр М.І., Русин Б.П. Дослідження впливу параметрів схеми зйому проекційних даних на якість реконструкції зображень // Вісник Національного університету “Львівська політехніка”: “Радіоелектроніка та телекомунікації”. 2000. №399. С. 143-149.

Надійшла до редколегії 21.12.2007

Рецензент: д-р техн. наук, проф. Путятін Є.П.

Кобасяр Михайло Іванович, канд. техн. наук, м.н.с. відділу “Методи та системи аналізу, обробки та ідентифікації зображень” Фізико-механічного інституту ім. Г. В. Кар-пенка НАНУ. Наукові інтереси: обробка та розпізнавання зображень. Адреса: Україна, 79601, Львів, вул. Наукова, 5а. Email: dep3 [email protected].

Русин Богдан Павлович, д-р. техн. наук, проф., зав. відділом “Методи та системи аналізу, обробки та ідентифікації зображень” Фізико-механічного інституту ім. Г. В. Кар-пенка НАНУ. Наукові інтереси: обробка та розпізнавання зображень. Адреса: Україна, 79601, Львів, вул. Наукова, 5а. Email: [email protected].

УДК517.988.63

КРАЕВЫЕ ЗАДАЧИ ДЛЯ НЕЛИНЕЙНОГО УРАВНЕНИЯ ТРЕТЬЕГО ПОРЯДКА

ЧЕРНЫШОВ Н.Н._______________________________

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

1. Введение

Целью работы является доказательство существования и единственности решения краевых задач для нелинейного уравнения смешанного типа третьего порядка. В области D, ограниченной отрезками прямых X=0, Y=H=-X=X-L, рассмотрим задачу о сопряжении двух нелинейных уравнений [1]:

Uxxy = F(X,Y,U); (X,Y) є D1 ^

^ D n (Y > 0);

' Uxxx = -Uxyy ^ F2(X,Y,U); (X,Y) є D2 ^ (1)

^ D n (Y < 0).

Эта система является ур авнением смешанного типа в области D [2]. Прямая Y=0 является характеристикой системы уравнений (1).

lU(X) = a(X)U(X) + y(X),0 < X < L;

[UY (X) = P(X)UY (X) + 5(X), 0 < X < L, (2)

где a(X), P(X), 5(X), Y(X) - заданные гладкие функции.

Система (1) соответствует второму и третьему типу канонического вида.

Задачи исследования:

Найти в области D\(Y=0) решения (1), удовлетворяющие краевым условиям:

U(Y) = ф1 (Y),U(L,Y) = ф2 (Y), 0 < Y < H; U(X-X) = ^(X);

• ux(x,-x) = v 2(X),0 < X < l/2;

<Pi(Y) Є C(H),i = 1,2;

@1(X) є C3(0,L/2); у2(X) є C2(0,L/2).

Введем следующие обозначения [5]:

U(X) = T1(X);Uy(X) = V1(X);

• U(X) = t 2(x);Uy(X) = v 2(X);

0 < X < L.

где t i (X), v i (X) - неизвестные функции.

Задача сводится к решению следующих вспомогательных задач.

- Найти функцию U (X, Y), удовлетворяющую в области D1;2 системе (1):

U(X) = t1;2(X), 0 < X < L.

РИ, 2007, № 4

111