1 Анализ на възстановяването на скоростта на разпадане в Т-претеглени изображения на мозъка 2 на магнитния резонанс Родни Харамильо Юстинико Университет Национал де Колумбия Меделин Кампус Факултет на науките Аспирант по математика юни 2014 г.

скоростта

3 Анализ на възстановяването на скоростта на разпадане в T2-претеглени мозъчни изображения на магнитния резонанс От: Родни Джарамило Юстинико Работата е представена като частично изискване за кандидатстване за званието доктор на математическите науки Директор: Марианела Лентини Гил Национален университет в Колумбия Меделин Кампус Факултет на Науки Аспирант по математика юни 2014 г.

4 Тази работа е частично подкрепена от заместник-ректората за научни изследвания чрез проекта Укрепване на групата за научни изчисления, код на Хермес 16084

5 Благодарности Благодарен съм на моите колеги и приятели, учители по математически и статистически училища, които предложиха насърчителни думи за реализирането на този проект, особено Карлос Мехия, Марко Палушни, Хуго Арбелаес и Хуан Карлос Салазар. благодарение на Беатрис Корея, без чието настояване не бих възобновил докторантурата си, имам много специално чувство на благодарност към моя съветник, професор Марианела Лентини, за нейната отдаденост и особено за нейните учения, например доверие и приятелство. И накрая, аз благодаря на безусловната привързаност на хората, които с нетърпение очакваха, а понякога и с нетърпението на любимия, завършването на тази теза: майка ми Гладис, баща ми Лудоберто, съпругата ми Олга Росио и децата ни Самуел и Хуана

9 Съдържание Резюме Резюме Съдържание i ii iii Въведение 1 1 Вариант на метода Прони 3 11 Методи тип Прони 3 12 Вариант на метода Прони 5 2 Анализ на стабилността Числени симулации 17 3 Филтри в областта на вейвлетите за намаляване на шума в ядрено-магнитен резонанс Внедряване на филтри в областта на вейвлет за ядрено-магнитен резонанс Премахване на пристрастия за данни след разпределение на ориз Внедряване на нов филтър в областта на вейвлет за ядрено-магнитен резонанс Формула за коефициенти на мащаб Филтър тип Wiener за вълнови коефициенти Двустранен филтър Алгоритъм за намаляване на шума при изображения с магнитен резонанс Валидиране на филтъра с помощта на синтетични изображения Изпълнение на филтъра върху реално изображение с магнитен резонанс 37 4 Числени резултати Приложение на метода върху претеглени изображения в Real T 2 MRI Re Числени резултати върху синтетични изображения Заключения и обсъждане на резултатите 53 iii

10 Библиография 54 iv

13 Глава 1 Вариант на метода на Прони 11 Методи от типа Прони Методите от тип Прони представляват семейство от методи, които позволяват, наред с други проблеми, експоненциалната корекция, дадена от системата от уравнения kyi = b + C je iλ jtj = 1 i = 1, n, Ако в формулировката, дадена от (1), са дефинирани b = C 0 и λ 0 = 0, тогава данните yi трябва да отговарят на модела µ (ti) = µ (it) yi, където µ е даденото функция чрез µ (t) = k C je λjt j = 0 Тези методи, известни също като полиномиални методи, се характеризират, тъй като µ (t) удовлетворява различното уравнение на формата (δk + 2 E k δ 2 E + δ 1) µ (t) = 0, (2) където операторът E е даден с (Eµ) (t) = µ (t + t) и стойностите β j = e λ jt са корените на полинома P (z ) = δ k + 2 zk δ 2 z + δ 1 = 0, (3) което е характеристичният полином, свързан с уравнението на разликата (2) Когато се изчислява (2) за ti = it, i = 1, nk 1, получаваме набора от уравнения δ k + 2 µ (t k + 2) + + δ 1 µ (t 1) = 0 δ k + 2 µ (tn) + + δ 1 µ (t n k 1) = 0 3

17 В този случай коефициентите на полинома α (z) са симетричните функции в β 1, β k, определени от w (k) 1 = β β kw (k) 2 = β l β rlrw (k) 3 = lr, ls, srw (k) k 1 = (1) k β l β r β sk (k) β lj = 1 w (k) k = (1) k + 1 ljk β j, тези коефициенти се изчисляват като решение на системата от уравнения Накрая, β j са корените на полинома j = 1 M (k) w (k) = Q (k) (12) α (k) (z) = zkkj = 1 w (k) jzkj (13) Двете теореми, които ще посочим по-долу, установяват връзката, която съществува между решението, получено по процедурата, която току-що описахме, и модифицирания метод на Прони, описан в раздел 11 Теорема 1 Нека R е матрицата от порядъка kk, определена както следва: R = 1, ако k = 1, а за k> 1 1, ако i = j, R (i, j) = 1, ако j = i + 1, 0, в противен случай нека P (z) и α (k) (z) полиномите, дефинирани съответно в (3) и (13), ако δ = [δ 1, δ k + 1, 1] е решението на задачата за оптимизация (9), тогава векторът w (k ) = R 1 [δ k, δ 1] T удовлетворява Освен това, M (k) w (k) Q (k) = XT δ и P (z) = (z 1) α (k) (z) Тест Разтворът δ = [δ 1, δ k + 1, δ k + 2] от (9) удовлетворява δ k + 2 = 1 В случая, който разглеждаме β 0 = 1 е корен от P (z), от който следва, че δ j = 1 k + 1 j = 1 7

18 Тогава M (k) w (k) Q (k) = M (k) R 1 Rw (k) Q (k) = M (k) R 1 δ ḳ y k + 1 и k + 2 δ 1 yn 1 yn = M (k) R 1 δ ḳ δ 1 [k + 1 j = 1 δ j] и k + 1 yk + 2 [k + 1 j = 1 δ j] yn 1 yn = M (k) R 1 δ ḳ + yk + 1 и k + 1 δ ḳ + δ 1 yn 1 yn 1 δ 1 δ k + 2 и k + 2 yn + δ k + 1 и k + 1 yn 1 = δ k + 2 и k + 2 yn + δ k + 1 и k + 1 yn 1 + M (k) R 1 δ ḳ + yk + 1 и k + 1 δ ḳ δ 1 yn 1 yn 1 δ 1 = δ k + 2 и k + 2 yn + δ k + 1 yk + 1 yn 1 + ykyk 1 y 1 yn 2 yn 1 ynk 1 δ ḳ δ 1 8

19 = y 1 и 2 и k + 2 и 2 и 3 и k + 3 и 3 и 4 и k + 4 ynk 1 ynkyn δ 1 δ 2 δ k + 1 δ k + 2 = W и δ От уравнение (6) следва, че сега за полинома P (z) имаме P (z) = δ k + 2 zk δ 2 z + δ 1 ((k = (z 1) zk = (z 1) = (z 1) (( j = 1 M (k) w (k) Q (k) = XT δ y (k 1 δ j) zk 1 j = 1 δ j) zk 2 (δ 1 + δ 2) z δ 1) zkw (k) 1 zk 1 w (k) 2 zk 2 w (k) k 1 zw (k) kzk = (z 1) α (k) (z) kj = 1 w (k) jzkj)) Теорема 2 Да предположим, че има само решение на задачата за оптимизация (9) Векторът δ R k + 2 е решението на задача (9) тогава и само ако R 1 [δ k, δ 1] T е решението с най-малките квадрати на линейното уравнение (12) Доказателство Нека δ R k + 2 е решението на задача (9), ζ = R 1 [δ k, δ 1] T и нека ψ е най-малкото квадратно решение на линейната система (12) От теорема 1 следва, че XT δ y = M (k) ζ Q (k) min M (k) z Q (k) z = M (k) ψ Q (k) (14) Нека разгледаме ξ R k, дадено от [ξ k, ξ 1] T = Rψ (15) и γ R k + 2, дефинирани като k γ = [ξ 1, ξ k, 1 ξ j, 1] T (16) j = 1 По теорема ( 1) имаме M (k) ψ Q (k) = Xγ T и 9

20 Тогава XT δ и XT γ и По хипотеза δ е единственото решение на задачата за оптимизация, тогава δ = γ, което предполага, че [δ k, δ 1] = [ξ k, ξ 1] Тогава имаме Rζ = Rψ, откъдето ζ = ψ По подобен начин можем да покажем, че γ в (16) е решението на задачата за оптимизация (9), при условие че ψ е решението в най-малките квадрати на линейната система (12) 10

29 Фигура 22: Средна стойност на относителните грешки, след 100 пробега, за първия модел и съответстващи на гауссов шум Фигура 23: Относителни грешки за първия модел, съответстващи на рициев шум 19

30 Фигура 24: Средна стойност на относителните грешки, след 100 пробега, за първия модел и съответстващи на шум от Рициан Фигура 25: Относителни грешки за втория модел, съответстващи на шум от Рика 20

31 Фигура 26: Средна стойност на относителните грешки, след 100 пробега, за втория модел и съответстващи на рициев шум Фигура 27: Относителни грешки за третия модел, съответстващи на рициев шум 21

32 Фигура 28: Средна стойност на относителните грешки, след 100 пробега, за третия модел и съответстващи на шум на Рициан Номер на състоянието η G η D ˆη B модел модел на модел

37 От друга страна, ако x 375 тогава () () xex I 0 (x) = xx () ɛ 1 yx () () xex I 1 (x) = xx () ɛ 2, x където ɛ 1 19 (10 ) 7 y ɛ 2 22 (10) 7 Това полиномно представяне може да се намери в текста на Абрамовиц-Стегун, [1], страница 378 Тогава ако x 15, след това xy следователно (x 2 4 ex 2 4 I0 (x2 4) = x ()) x 2 2 ex 2 x 2 4 I0 4 (= x () () xx 2 ()) ɛ 1 x 2 Тогава x 2 lim x По същия начин се доказва, че x 2 lim x 4 ex 2 4 I 0 ( x2 x 4 ex 2 4) 4 I 1 (x2 x 4) = ɛ 1 (37) 2 = ɛ 2 (38) 2 Факторът v (x) x може да бъде записан като v (x) x = [x π e I 0 (x2 x 4) + 2 x 2 4 ex 2 4 I 0 (x2 x 4) + 2 x 2 4 ex 2 4 I 1 (x2 x 4)] (39) От равенствата (36), (37 ), (38) и (39) следва, че v (x) π L: = lim xx = 2 (ɛ 1 + ɛ 2) E r, където E r π 2 () (10) (10) 7 27

43 Фигура 32: Синтетично изображение, генерирано с MATLAB 34 Проверка на филтъра с помощта на синтетични изображения В този раздел се прави сравнение между производителността на филтъра, първоначално разработен от Kazubek, и модифицираната версия. и шумни изображения се генерират чрез добавяне на оризов шум към синтетичното изображение Шумът се генерира от уравнението J e (m, n) = (J (m, n) + e 1) 2 + e 2 2, (45) където J (m, n) е стойността без шум и e 1 и e 2 са произволни числа, които съответстват на разпределение на Гаус с нулева средна стойност и стандартно отклонение σ Важно е да се отбележи, че нивата на сивото в изображението са между 0 и 88 Five нивата на σ се използват, σ = [1, 2, 5, 8, 12] Петте индекса, които се появяват най-често в литературата, бяха разгледани за сравняване на ефективността на филтрите: съотношение сигнал към шум (SNR), пиков сигнал до съотношение на шума (PSNR), средно квадратична грешка (RMSE), средна абсолютна ер ror (MAE) и структурен подобен индекс (SSIM) Дадени две изображения с размери n x n y, референтно изображение r (x, y) и второ изображение t (x, y), което може да се разглежда като нарушение на първото; количествата PSNR, RMSE и MAE са дадени с 33

45 σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки Модификация на алгоритъма на Kazubek към алгоритъма на Kazubek Таблица 31: SNR (r, t) между оригиналното изображение и филтрираното изображение σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки Модификация на алгоритъма на Казубек към алгоритъма на Казубек Таблица 32: PSNR (r, t) между оригиналното изображение и филтрираното изображение σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки Модификация на алгоритъма на Казубек към алгоритъма на Казубек Таблица 33: RMSE (r, t) между оригиналното изображение и филтрираното изображение σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки Модификация на алгоритъма на Казубек към алгоритъма на Таблица 34: MAE (r, t) между оригиналното изображение и филтрираното изображение σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки Модификация на алгоритъма на Казубек към алгоритъма на Казубек Таблица 35: SSIM (r, t) между оригиналното изображение и филтрираното изображение 35

46 σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки 326% 652% 1625% 2616% 3917% Алгоритъм на Казубек 211% 399% 938% 1441% 1947% модификация на алгоритъма на Казубек 205% 379% 857% 1305% 1731% Таблица 36: Процентна грешка в изчислението на RMSE между изображенията и оригиналното синтетично изображение σ = 1 σ = 2 σ = 5 σ = 8 σ = 12 оригинални грешки 442% 882% 2207% 3537% 5294% Алгоритъм на Казубек 237% 444% 1055% 1626% 2111% модификация на алгоритъма на Казубек 230% 423% 1000% 1526% 1948% Таблица 37: Процентна грешка при изчисляването на MAE между изображенията и оригиналното синтетично изображение Фигура 33: Над засегнатото изображение чрез шум, съответстващ на σ = 2 Под филтрираното изображение, намерено чрез модифициране на филтъра Kazubek 36

47 35 Изпълнение на филтъра върху реално MRI изображение Фигура 34: MRI мозъчно изображение и две избрани области от интерес Мозъчното изображение, показано по-долу, е реално MRI изображение, върху което ще бъде приложен филтърът, който сме приложили, обяснен в този раздел, да оценим ефективността му върху изображения, които не идват от симулации 1 Чрез количествено определяне на разликите между оригиналното изображение и филтрираното изображение се получават следните стойности за региона, обозначен в горната част RMSE = 4287 MAE = 3399 SSIM = 0938, и за района, обозначен отдолу RMSE = 4261 MAE = 3396 SSIM = благодарен съм на професор Томас М Дезерно от Uniklinik RWTH, Аахен, Германия, за предоставяне на това и други клинични изображения 37

48 Фигура 35: Областта на интерес в оригиналното изображение е показана вляво Областта на интерес в оригиналното изображение е показана вдясно Фигура 36: Областта на интерес в оригиналното изображение е показана вдясно. 38. снимка

51 Фигура 41: Отгоре, отляво и отдясно, съответно първото и третото ехо Под петото ехо Фигура 42: Графиката вдясно показва разпаданията в нивото на интензитета, които съответстват на точките, посочени в графиката отляво 41

52 Фигура 43: Областта на интерес е показана в горната част. В долната лява част е графиката, получена при прилагане на варианта на метода Прони. В долната дясна част е показано решението, получено при прилагане на филтъра, върху ехото, преди да се използва вариантът на метода Прони Фигура 44: Областта на интерес е показана в горната част В долната лява част е графиката, получена чрез прилагане на варианта на метода Прони В долната дясна част решението, получено чрез прилагане на филтъра е показано в ехото, преди да се използва вариантът на метода Prony 42

53 Фигура 45: Изображение на шестото ехо ROI1 регион отгоре и ROI2 регион под 43

54 Фигура 46: Функции на плътност на вероятността, съответстващи на ROI1 и ROI2 Отгоре: резултат, отчетен в [19] Средно: Разтвор, получен по метод на Прони без процес на филтриране Отдолу: Разтвор, получен чрез прилагане на филтъра върху изображенията преди използване на метода Prony 44

56 Параметри на последователността mri2 ехо пъти ъгъл на обръщане 30 тип изображение M inu-поле A брой ехо 1 процент inu 0 процента шум 3 фантом: msles3 произволно семе 0 референтна тъкан 0 сканираща техника SFLASH дебелина на среза 2 ti tr: 2500 параметри на последователността mri3 ехо пъти ъгъл на обръщане 30 тип изображение M inu-поле A брой ехо 1 процент inu 0 процента шум 1 фантом: msles3 произволно семе 0 референтна тъкан 0 сканираща техника SFLASH дебелина на среза 2 ti tr:

57 Параметри на последователността mri4 ехо пъти ъгъл на обръщане 30 тип изображение M inu-поле A брой ехо 1 процент inu 0 процента шум 0 фантом: msles3 произволно семе 0 референтна тъкан 0 сканираща техника SFLASH дебелина на среза 2 ti tr: 2500 mri1 mri2 mri3 mri4 процента шум: 3 процента шум: 3 процента шум: 1 процент шум: 0 процента ину: 20 процента ину: 0 процента ину: 0 процента ину: 0 Таблица 42: Различия в нивата на шума на синтетичните изображения, генерирани от Brain-Web 47

58 Фигура 47: Изображения на първото, четвъртото и петото ехо, на mri1 последователността 48

59 Фигура 48: Изображения, поискани от BrainWeb Region ROIw1 по-горе, ROIw2 в центъра и ROIw3 под 49