Комп’ютерне моделювання методів реконструктивної томографії

Мета роботи: набуття практичних навичок з моделювання паралельних алгебраїчних методів реконструктивної томографії на послідовних ЕОМ.

1. Ознайомитись з методами відновлення (реконструкцією) зображення та вирішенням на їх основі такої прикладної задачі, як задачі реконструктивної комп”ютерної томографії, користуючись рекоме­н­до­­ваною літературою.

2. Скласти блок-схему паралельного алгоритму розв’язання системи лінійних алгебраїчних рівнянь (СЛАР) з постійною матрицею коефіцієнтів А розміром N*N і різними векторами вільних членів b_i (i=1,N) за запропонованою паралельною моделлю метода Гаусса-Жордана..

3. За розробленою в п.2 блок-схемою скласти програму обчислення розв’язку СЛАР, обраного за своїм варіантом із таблиці 4.1, використовуючи алгоритмічну мову за бажанням студента.

4. Виконати розв’язання вказаної в таблиці 4.1 СЛАР засобами Mathcad 2000 і порівняти з результатом, отриманим в попередньому пункті.

5. Роздрукувати початкові дані та результати в вигляді числових матриць, поданих в десятковій формі.

6. Зробити висновки щодо можливості комп’ютерної реалізації паралельних алгоритмів на комп”ютерах з традиційною архітектурою.

7. Оформити звіт по роботі.

Теоретичні відомості

1. Постановка задачі реконструктивної томографії та алгебраїчний підхід до її вирішення

Існує значний клас прикладних задач, в яких необхідно здійснити синтез зображення за деякими його ознаками: виміряними спектрами просторових частот, радіаційному випромінюванню, поглинанню рентгенівських променів, відбиттю, розсіюванню ультразвука тощо. Методи відновлення (реконструкції) зображення суттєво залежать від фізичного процесу і виміряного параметру.

Суттєві успіхи досягнуті в останній час в медичній реконструктивній томографії [1]. Поява скануючих реконструктивних комп’ютерних томографів здійснила повний переворот в області рентгенівських методів отримання зображення. Виник новий рентгенодіагностичний метод відновлення (реконструкції) за допомогою ЕОМ структури поперечного анатомічного перерізу тіла пацієнта – реконструктивна комп’ютерна томографія.

Як показано на рис.1.1, тіло пацієнта опромінюється значною кількістю М паралельних пучків променів, що знаходяться в одній площині. Енергія кожного з них відома як на вході джерела випромінювання, так і на виході матриці детекторів. Логарифм відношення вихідної енергії до вхідної для кожного пучка дає загальне інтегральне поглинання і-го пучка () вздовж його шляху в тілі пацієнта. Метою є обчислення локальних значень коефіцієнтів поглинання x_j () в усіх точках поперечного перерізу зображення, які відображають стан різних тканин тіла пацієнта в заданій площині. Позначимо через а_ij довжину шляху і-го променя по j-му елементу дискретного растру, яка відповідає вазі внесення j-го елемента в загальне поглинання вздовж всього шляху розповсюдження і-го променя в об’єкті дослідження (див. рис.1.1). Тоді дискретна модель процесу може бути подана в вигляді СЛАР типу .

Таким чином, дискретна модель розв’язання задачі реконструкції зображення характеризується СЛАР, яку можна подати в матричній формі так:

Ax=b, (1)

де А – невироджена матриця проекції a_ij розміром L´N;

x=(x₁,x₂,…,x_N) – вектор елементів зображення;

b=(b₁,b₂,…,b_L) – вектор вимірювань.

Для реконструкції дискретного зображення необхідно розв’язати матричне рівняння

x=A^–1b, (2)

де А^-1 – матриця, обернена до матриці А.

Розмірність системи (1.2) дуже велика, оскільки для отримання зображення з хорошою роздільною здатністю параметри N і L повинні мати порядок 10⁵. Це ставить високі системні вимоги до швидкодії обчислювальних засобів (час обробки однієї проекції £ 10 мс), що є складовою структурної схеми системи реконструктивної томографії [1,2], які важко задовільнити в межах традиційних засобів обчислювальної техніки.

В пошуках методів, орієнтованих на нетрадиційну обчислювальну техніку, а саме оптоелектронний паралельний комп”ютер, і була запропонована розглянута нижче паралельна інтарпретація метода Гаусса-Жордана для розв”язання СЛАР високих порядків та обернення великорозмірних матриць.

Рис.1.1. Схематичне зображення реконструктивної томографії

2. Паралельна інтерпретація метода Гаусса-Жордана для розв”язання СЛАР високих порядків

З традиційним методом Гаусса-Жордана для розвязання СЛАР можна ознайомитись в літературі, наприклад, [3,4].

Розглянемо нову паралельну інтерпретацію метода Гаусса-Жордана на основі обчислення зовнішнього добутку векторів, запропоновану в [5].

Для розв”язання СЛАР виду (2) та обернення матриць початкові дані подають в вигляді розширеної матриці [A|B], де А – квадратна неособлива матриця коефіцієнтів розмірністю N*N, В – квадратна матриця вільних членів розмірністю N*N, яка в випадку розв”язання СЛАР формується за вектор-стовпцями вільних членів розмірності N кожний В = [b⁽¹⁾ ,b⁽²⁾ , …b^(N) ], а в випадку знаходження оберненої матриці співпадає з одиничною матрицею В=Е .

Розширену матрицю проміжних результатів позначимо як R[1:N;1:2N], а матрицю кінцевого результату як X[1:N;1:N].

Для t=0 в розширену матрицю проміжного результату R[1:N;1:2N] будуть записані початкові дані таким чином:

R^(t=0)[1:N;1:N]=A, R^(t=0)[1:N;(N+1):2N]=B. (3)

За ведучий елемент на кожному t-му кроці обчислень обраємо елемент r[t,t]¹0 матриці R.

Розділивши одночасно коефіцієнти t-го рівняння системи, що є елементами t-го рядка матриці R, на ведучий елемент, отримаємо вектор-рядок d[1;1:2N], а саме:

(4)

Визначивши вектор-стовпець v^(t)[1:N;1] як виділений t-й стовпець матриці проміжного результату із одиничним значенням елемента, що знаходиться в t-му рядку,

v^(t)[1:N;1]= R^(t-1)[1:N;t], (5)

де v^(t)[t;1]=-1, сформуємо зовнішній добуток векторів в вигляді матриці P за формулою (6):

P^(t)[1:N;1:2N]= v^(t)[1:N;1]Äd^(t)[1;1:2N], (6)

де p^(t)[I;J]= v^(t)[I,1]´d^(t)[1,J], ().

Тоді процес формування проміжного результату R в часі описується таким матричним співвідношенням:

R^(t)=(R^(t-1)ÙMR^(t))- P^(t), (7)

де Ù – оператор логічного множення (кон’юнкція);

MR^(t) – специфічна бінарна матриця розмірності N´2N елементів, будь-який елемент mr[i,j] якої визначається таким чином:

(8)

P^(t) – зовнішній добуток векторів.

Кінцевий результат X[1:N;1:N], що являє собою або матрицю шуканих векторів x⁽¹⁾,x⁽²⁾,…,x^(L), або обернену матрицю А^-1, визначається на останньому t=N-му кроці обчислень в виді

X[1:N;1:N]=R^(t=N)[1:N;(N+1):2N]. (9)

Час виконання обчислень за вищенаведеною формою алгоритму Гаусса-Жордана для розв’язання СЛАР складає N кроків як в однозадачному, так і в потоковому режимі обчислень.

Таким чином, розроблено паралельну форму методу Гаусса-Жордана для розв’язання СЛАР, яка відтворює в просторі та часі процес обчислень за допомогою матричних рівнянь.

Таблиця 4.1.

N	A	b
1.	0.21 0.30 0.60	–0.45 0.25 -0.35	–0.20 0.43 -0.25	1.91 0.32 1.83
2.	-3.00 0.50 0.50	0.50 -6.00 0.50	0.50 0.50 -3.00	-56.50 -100.00 -210.00
3.	0.45 -0.01 -0.35	-0.94 0.34 0.05	-0.15 0.06 0.63	-0.15 0.31 0.37
4.	0.63 0.15 0.03	0.05 0.10 0.34	0.15 0.71 0.10	0.34 0.42 0.32
5.	-0.20 -0.30 1.20	1.60 0.10 -0.20	-0.10 -1.50 0.30	0.30 0.40 -0.60

N	A	b
6.	0.30 -0.10 0.05	1.20 -0.20 0.34	-0.20 1.60 0.10	-0.60 0.30 0.32
7.	0.20 0.58 0.05	0.44 -0.29 0.34	0.81 0.05 0.10	0.74 0.02 0.32
8.	6.36 7.42 5.77	11.75 19.03 7.48	10.00 11.75 6.36	-41.40 -49.49 -27.67
9.	-9.11 7.61 -4.64	1.02 6.25 1.13	-0.73 -2.32 -8.88	-1.25 2.33 -3.75
10.	1.02 6.25 1.13	-0.73 -2.32 -8.88	-9.11 7.62 4.64	-1.25 2.33 -3.75

Література

1. Бутаков Е.А. и др. Обработка изображений на ЭВМ .-М.: Радио и связь, 1987. – 240 с.

2. СБИС для распознавания образов и обработки изображений / Под ред. К.Фу.- М.: Мир, 1988- 248 с.

3. Методичні вказівки до виконання лабораторних робыт з дисципліни ‘Чисельні методи в інженерній діяльності’ для студентів бакалаврату 6.0911/ Мартинюк Т.Б., Заболотна Н.І., Шолота В.В.- Вінниця. 1997, 68 стор.

4. Лященко М.Я., Головань М.С. Чисельні методи. – Київ:Либідь, 1996 – 284 с.

5. Шолота В.В. Високопродуктивний процесор для паралельного розв’язання СЛАР та обернення матриць // ВОТТП (Хмельницький).-1999. – №1.-С.86-91.