ˇ ˇ ´I TECHNICKE´ V PRAZE CESK E´ VYSOKE´ UCEN ´ A FYZIKALN ´ ˇ YRSK ´ ´ FAKULTA JADERNA Eˇ INZEN A
´ PRACE ´ DIPLOMOVA
2010
Ondˇrej Tich´y
ˇ ˇ ´I TECHNICKE´ V PRAZE CESK E´ VYSOKE´ UCEN ´ A FYZIKALN ´ ˇ YRSK ´ ´ FAKULTA JADERNA Eˇ INZEN A Katedra matematiky
´ PRACE ´ DIPLOMOVA
Integr´aln´ı modely dynamick´ych scintigrafick´ych studi´ı ledvin ( Integral models for dynamic renal scintigraphy )
Ondˇrej Tich´y
ˇ ˇ ıdl, Ph.D. Skolitel: Ing. V´aclav Sm´ ˇ amal, DrSc. Konzultant: Prof. MUDr. Martin S´ Akademick´y rok: 2009/2010
Sem d´at zad´an´ı DP, jeden origin´al a dvˇe kopie...
ˇ Cestn´ e prohl´asˇen´ı Prohlaˇsuji na tomto m´ıstˇe, zˇ e jsem pˇredloˇzenou pr´aci vypracoval samostatnˇe a zˇ e jsem uvedl veˇskerou pouˇzitou literaturu.
V Praze dne 7. kvˇetna 2010
.......................... Ondˇrej Tich´y
Podˇekov´an´ı Chtˇel bych na tomto m´ıstˇe podˇekovat sv´ym rodiˇcu˚ m za podporu v cel´em m´em studiu a obzvl´asˇtˇe pak tat´ınkovi Miroslavovi za korekturu textu. ˇ amalovi, za poskytnut´e scintigrafick´e studie, D´ık patˇr´ı tak´e m´emu konzultantovi, prof. S´ pomoc s kapitolou o Scintigrafii a trpˇeliv´e konzultace, kter´e usmˇernily smˇer v´yzkumu. ˇ ıdlovi, za pˇr´ıkladn´e veden´ı Nejvˇetˇs´ı podˇekov´an´ı vˇsak patˇr´ı m´emu sˇkoliteli, Dr. V´aclavu Sm´ t´eto diplomov´e pr´ace, hodiny konzultac´ı, cenn´e n´apady, pˇripom´ınky a za v´yborn´e uveden´ı do cel´e problematiky.
N´azev pr´ace: Integr´aln´ı modely dynamick´ych scintigrafick´ych studi´ı ledvin Autor: Ondˇrej Tich´y Obor: Inˇzen´yrsk´a informatika Druh pr´ace: Diplomov´a pr´ace ˇ ıdl, Ph.D., Ustav ˇ ´ Vedouc´ı pr´ace: Ing. V´aclav Sm´ teorie informace a automatizace AV CR ˇ ´ Konzultant: Prof. MUDr. Martin S´amal, DrSc., Ustav nukle´arn´ı medic´ıny 1.lf UK v Praze Abstrakt: Tato pr´ace se zab´yv´a konvoluˇcn´ı parametrizac´ı v modelu pro faktorovou anal´yzu ve scintigrafii ledvin. Souˇca´ st´ı parametr˚u tohoto nov´eho modelu jsou pˇr´ım´e diagnostick´e koeficienty, nen´ı tedy nutn´a n´asledn´a anal´yza a interakce s odborn´ym uˇzivatelem. Po pˇredstaven´ı scintigrafie a Bayesovsk´e teorie se pr´ace vˇenuje popisu a ˇreˇsen´ı standardn´ıho modelu pro faktorovou anal´yzu a n´aslednˇe na re´aln´e studii ukazuje jeho funkci. Protoˇze vˇsak standardn´ı model nen´ı zaloˇzen na re´aln´ych biologick´ych pˇredpokladech, nar´azˇ´ıme na jeho omezen´ı pˇri z´ısk´av´an´ı diagnostick´ych parametr˚u. Hlavn´ım v´ysledkem t´eto pr´ace je navrˇzen´ı integr´aln´ıho modelu, kter´y respektuje biologick´a omezen´ı a z´aroveˇn pˇr´ımo odhaduje diagnostick´e koeficienty jako sv´e parametry. Tento nov´y model je ˇreˇsen pomoc´ı Variaˇcn´ıho Bayesova teor´emu a otestov´an na re´aln´ych scintigrafick´ych studi´ıch. Zabudov´an´ı nov´ych informac´ı pˇrineslo poˇzadovanou automatizaci diagnostiky a je uk´az´ano velk´e zlepˇsen´ı oproti standardn´ımu modelu pro faktorovou anal´yzu. Kl´ıcˇ ov´a slova: scintigrafie, Bayesovsk´a statistika, nukle´arn´ı medic´ına, obrazov´a sekvence, konvoluˇcn´ı parametrizace
Title: Integral models for dynamic renal scintigraphy Author: Ondˇrej Tich´y Abstract: This work is concerned with a convolution based parametrization of model for factor analysis in renal scintigraphy. Diagnostic coefficients are used as parameter of the new integral model, hence their estimates are direct output of the estimation and no follow-up analysis and interaction with expert user is required. After introduction of the field of scintigraphy and Bayesian theory, the Variational Bayes solution of standard model for factor analysis is described in detail. Suitability of the model is tested on real data. It is shown that the estimated diagnostic coefficients are inadequate since the model do not respect properties of real biological systems. The main contribution of this work is the integral model and estimation of its parameters via the Variational Bayes approximation. The model is designed to respect biological constraints and its parameters are variables used as diagnostic coefficients. The resulting estimation algorithm is tested on real scintigraphic data. It is shown that resulting estimates are much more realistic that of the standard model. Moreover, the estimates of the diagnostic coefficient were achieved without any interaction with expert user. Key words: scintigraphy, Bayesian statistic, nuclear medicine, image sequence, convolutionbased parameterization
Obsah ´ 1 Uvod
1
ˇ ren´ı 2 Scintigraficke´ vysetˇ 2.1 Scintigrafie . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Pr˚ubˇeh vyˇsetˇren´ı . . . . . . . . . . . . . . . . . . 2.1.2 Radionuklidy . . . . . . . . . . . . . . . . . . . . 2.1.3 Scintilaˇcn´ı kamera . . . . . . . . . . . . . . . . . 2.2 Biologick´e pˇredpoklady . . . . . . . . . . . . . . . . . . . 2.2.1 Tvary faktorov´ych kˇrivek . . . . . . . . . . . . . . 2.2.2 Faktorov´e obr´azky . . . . . . . . . . . . . . . . . 2.3 Souˇcasn´e metody funkcion´aln´ı anal´yzy obrazov´e sekvence 2.3.1 Oblasti z´ajmu . . . . . . . . . . . . . . . . . . . 2.3.2 Matematick´y pohled a vyhodnocen´ı kˇrivek . . . . 2.3.3 Klinick´e vyhodnocen´ı . . . . . . . . . . . . . . . 2.4 Zhodnocen´ı . . . . . . . . . . . . . . . . . . . . . . . . . 3 Bayesovska´ teorie a aproximace 3.1 Z´aklady Bayesovsk´e teorie . . . . . . . . . . . . . . . . . 3.1.1 Volba apriorn´ıho rozdˇelen´ı . . . . . . . . . . . . . ˇ sen´ı pomoc´ı aproximace . . . . . . . . . . . . . . . . . 3.2 Reˇ 3.2.1 EM (expectation minimalization) algoritmus . . . 3.2.2 Laplaceova aproximace . . . . . . . . . . . . . . . 3.2.3 Markov Chain Monte Carlo (MCMC) aproximace 3.2.4 Variaˇcn´ı Bayesova aproximace . . . . . . . . . . . 4 Faktorova´ analyza ´ obrazovych ´ sekvenc´ı 4.1 Konstrukce standardn´ıho modelu pro sekvenci sn´ımk˚u 4.1.1 Pozn´amka k maticov´e dekompozici . . . . . . 4.2 VB metoda pro standardn´ı model (4.1) . . . . . . . . . 4.3 V´ysledky standardn´ıho modelu . . . . . . . . . . . . . 4.3.1 Faktorov´e obr´azky a jejich kˇrivky . . . . . . . 4.3.2 Tranzitn´ı cˇ as . . . . . . . . . . . . . . . . . . 4.3.3 Relativn´ı ren´aln´ı clearance . . . . . . . . . . . 4.4 Zhodnocen´ı v´ysledk˚u . . . . . . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
´ ı model pro funcionaln´ ´ ı analyzu 5 Integraln´ ´ obrazovych ´ sekvenc´ı
ii
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . . .
2 2 3 3 6 6 7 8 9 10 11 13 14
. . . . . . .
16 16 17 17 18 18 18 19
. . . . . . . .
22 22 24 24 27 28 29 30 31 32
5.1
. . . . . . . . . . . . . . . . .
32 33 34 34 35 37 37 39 40 40 41 45 48 49 52 52 52
. . . . . . .
53 53 55 55 56 56 57 58
´ er ˇ a moˇznosti dals´ ˇ ıho pokracov ˇ an´ ´ ı 7 Zav 7.1 Hlavn´ı pˇr´ınos pr´ace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Moˇznosti dalˇs´ıho pokraˇcov´an´ı . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Zhodnocen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60 60 61 62
A Matematicky´ dodatek A.1 Stopa matice, oper´ator vec(), Kronecker˚uv a Hadamard˚uv souˇcin A.2 Norm´aln´ı rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . A.2.1 V´ıcerozmˇern´e norm´aln´ı rozdˇelen´ı . . . . . . . . . . . . A.2.2 Maticov´e norm´aln´ı rozdˇelen´ı . . . . . . . . . . . . . . . A.2.3 Vektorizace v maticov´em norm´aln´ım rozdˇelen´ı . . . . . A.2.4 Oˇrezan´e norm´aln´ı rozdˇelen´ı . . . . . . . . . . . . . . . A.3 Gamma rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Exponenci´aln´ı rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . A.4.1 Oˇrezan´e exponenci´aln´ı rozdˇelen´ı . . . . . . . . . . . . . A.5 Rovnomˇern´e rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . .
63 63 64 64 64 65 65 66 66 66 67
5.2
5.3
Konstrukce modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Z´akladn´ı model dat . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Model chyb mˇeˇren´ı . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.3 Modelov´an´ı matice obr´azk˚u . . . . . . . . . . . . . . . . . . . . 5.1.4 Modelov´an´ı impulzn´ıch retenˇcn´ıch funkc´ı . . . . . . . . . . . . . 5.1.5 Model kˇrivky krve . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.6 Odhadovan´e a z´avisl´e promˇenn´e . . . . . . . . . . . . . . . . . . ˇ sen´ı integr´aln´ıho modelu VB metodou . . . . . . . . . . . . . . . . . . Reˇ 5.2.1 Sestrojen´ı modelu . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 V´ypoˇcet logaritmu sdruˇzen´eho rozdˇelen´ı . . . . . . . . . . . . . 5.2.3 V´ypoˇcet VB-margin´al . . . . . . . . . . . . . . . . . . . . . . . 5.2.4 Identifikace standardn´ıch forem . . . . . . . . . . . . . . . . . . 5.2.5 Odhad indexov´ych mnoˇzin aprioren obrazov´e a pˇr´ır˚ustkov´e matice 5.2.6 Formulace VB-moment˚u . . . . . . . . . . . . . . . . . . . . . . Shrnut´ı v´ypoˇctu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Inicializace v´ypoˇctu . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Konvergence v´ypoˇctu . . . . . . . . . . . . . . . . . . . . . . . .
´ ıho modelu 6 Vysledky ´ integraln´ 6.1 Popis odhadovan´ych parametr˚u . . 6.1.1 Tranzitn´ı cˇ as . . . . . . . 6.1.2 Relativn´ı ren´aln´ı clearance 6.2 Vyhodnocen´ı dalˇs´ıch studi´ı . . . . 6.2.1 Studie 1 . . . . . . . . . . 6.2.2 Studie 2 . . . . . . . . . . 6.3 Zhodnocen´ı v´ysledk˚u . . . . . . .
. . . . . . .
. . . . . . .
iii
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . . . . .
. . . . . . . . . .
ˇ logaritmu sdruˇzeneho ´ ˇ B Vypo ´ cet rozdelen´ ı
68
C Obrazove´ sekvence
70
Literatura
72
iv
ˇ ı Pouˇzite´ znacen´ Line´arn´ı algebra R A ∈ Rp×n A′ A−1 ai ai,j resp. (A)i,j gi diag(A) diag(a) tr(A) vec(A) A⊗B A◦B In 1p,q 0p,q
mnoˇzina re´aln´ych cˇ´ısel A je re´aln´a matice rozmˇeru p × n; matice budeme obvykle znaˇcit velk´ymi p´ısmeny transpozice matice A inverze matice A i-t´y sloupec matice A cˇ len v i-t´em sloupci a j-t´em ˇra´ dku v matici A i-t´y prvek vektoru g diagon´ala cˇ tvercov´e matice A, v´ysledkem je tedy vektor s diagon´aln´ımi prvky matice A cˇ tvercov´a matice s prvky vektoru a na diagon´ale stopa matice A (viz A.1) vektorizace matice A (viz A.1) Kronecker˚uv souˇcin matic A a B o libovoln´em rozmˇeru (viz A.1) Hadamard˚uv souˇcin matic A ∈ Rp×n a B ∈ Rp×n (viz A.1) cˇ tvercov´a jednotkov´a matice rozmˇeru n, tedy (In )i,j = 1, pokud i = j a (In )i,j = 0, pokud i 6= j matice jedniˇcek o rozmˇerech p × r matice nul o rozmˇerech p × r Matematick´a anal´yza
erf(x) Γ(x) χ(a,b) ln(..) exp(..)
error funkce (viz A.2.4) gamma funkce (viz A.3) charakteristick´a funkce intervalu (a, b) (viz A.2.4) pˇrirozen´y logaritmus argumentu; je-li argumentem matice, rozum´ı se logaritmus po prvc´ıch exponenciela argumentu; je-li argument matice, rozum´ı se exponenciela po prvc´ıch Pravdˇepodobnostn´ı poˇcet
Ef (x) (..) b A Nx (µ, σ) Nx (µ, Σ) NX (M, Σp ⊗ Σn )
oˇcek´avan´a hodnota argumentu vzhledem k funkci f (x) bodov´y odhad parametru A norm´aln´ı rozdˇelen´ı skal´aru x se stˇredn´ı hodnotou µ a varianc´ı σ v´ıcerozmˇern´e (vektorov´e) norm´aln´ı rozdˇelen´ı vektoru x se stˇredn´ı hodnotou µ a Σ (viz A.2.1) maticov´e norm´aln´ı rozdˇelen´ı matice X ∈ Rp×n se stˇredn´ıho hodnotou M a kovarianˇcn´ımi maticemi Σp a Σn (viz A.2.2)
v
tNx (µ, σ, a, b) tNx (µ, σ) Gx (α, β) Expx (λ) t Expx (λ) Ux (a, b)
oˇrezan´e norm´aln´ı rozdˇelen´ı skal´aru x se stˇredn´ı hodnotou µ a varianc´ı Σ na intervalu (a, b] (viz A.2.4) oˇrezan´e norm´aln´ı rozdˇelen´ı skal´aru x se stˇredn´ı hodnotou µ a varianc´ı Σ na intervalu (0, +∞) gamma rozdˇelen´ı s parametry α, β (viz A.3) exponenci´aln´ı rozdˇelen´ı skal´aru x s parametrem λ (viz A.4) oˇrezan´e exponenci´aln´ı rozdˇelen´ı skal´aru x s parametrem λ na intervalu (0, +∞) (viz A.4.1) rovnomˇern´e rozdˇelen´ı skal´aru x na intervalu (a, b) (viz A.5)
vi
´ 1 Uvod V´yznamnou u´ lohu v l´ekaˇrsk´e diagnostice zauj´ım´a jedna z metod nukle´arn´ı medic´ıny, scintigrafie. Ta je zaloˇzena na sn´ım´an´ı distribuce radiofarmaka, aplikovan´eho do tˇela pacienta, d´ıky cˇ emuˇz lze z´ıskat nejen tvar, ale i funkci jednotliv´ych struktur v tˇele. Z´akladn´ım u´ kolem n´asledn´e anal´yzy je rozliˇsit na sekvenci sn´ımk˚u jednotliv´e obr´azky a pr˚ubˇeh kontrastn´ı tekutiny v nich, dohromady tzv. faktory. Svoj´ı podstatou vˇsak scintigrafie d´av´a relativnˇe m´alo kvalitn´ı sn´ımky. To je zp˚usobeno n´ızk´ym rozliˇsen´ım scintigrafick´e kamery (aˇz o ˇra´ d horˇs´ı neˇz magnetick´a rezonance) a probl´em˚um se sˇ umem (z okol´ı, z krevn´ıho ˇreˇciˇstˇe, z tk´an´ı, atd.). Neexistuje proto zˇ a´ dn´a automatick´a metoda, kter´a by provedla celkovou diagnostiku ze z´ıskan´e sekvence. Souˇcasn´e metody jsou znaˇcnˇe z´avisl´e na zruˇcnosti odborn´e obsluhy, mnohdy se v´ysledky mohou znaˇcnˇe liˇsit ´ [1]. Ukolem pˇredkl´adan´e pr´ace bude tento nedostatek odstranit. Hlavn´ı u´ lohou, kterou zast´av´a matematika ve zpracov´an´ı sekvenc´ı sn´ımk˚u ve scintigrafii, je anal´yza t´eto sekvence a spr´avn´e urˇcen´ı diagnostick´ych koeficient˚u. Protoˇze se analyzuje re´aln´y biologick´y syst´em, nen´ı moˇzn´e oˇcek´avat teoreticky pˇresn´e v´ysledky. V´yhodn´e metody pro pr´aci s t´ımto syst´emem n´am pˇrin´asˇ´ı Bayesovsk´a statistika. Umoˇznˇ uje uvlivˇnovat v´ypoˇcet dle poˇzadovan´ych pˇredpoklad˚u a zav´adˇet apriorn´ı pˇredpoklady. Cenou za tyto moˇznosti je ovˇsem analytick´a neˇreˇsitelnost [2], je nutno pˇristoupit k r˚uzn´ym form´am aproximace. V t´eto pr´aci je pouˇzita aproximace pomoc´ı Variaˇcn´ıho Bayesova teor´emu. Probl´emem souˇcasn´ych pˇr´ıstup˚u v t´eto oblasti je nedostateˇcn´e respektov´an´ı biologick´ych pˇredpoklad˚u, kter´e ovˇsem nelze opomenout, chceme-li dospˇet ke spr´avn´ym re´aln´ym v´ysledk˚um. V pˇredchoz´ı pr´aci autora [3] je naznaˇceno modelov´an´ı kˇrivek faktor˚u jako konvoluce krve s konvoluˇcn´ımi j´adry jednotliv´ych faktor˚u. To je v t´eto pr´aci rozvinuto a zdokonaleno. Pr´ace se d´ale vˇenuje vytvoˇren´ı integr´aln´ıho modelu pro anal´yzu scintigrafick´ych obrazov´ych sekvenc´ı, tzn. zabudov´an´ı potˇrebn´ych diagnostick´ych koeficient˚u pˇr´ımo do v´ypoˇctu. To vede k odstranˇen´ı chyby lidsk´eho faktoru, kter´a je v t´eto oblasti nezanedbateln´a, a z´aroveˇn ke zjednoduˇsen´ı a urychlen´ı cel´e diagnostiky v praxi. Na z´avˇer je cel´y postup a nov´y model otestov´an na re´aln´ych scintigrafick´ych studi´ıch ledvin.
1
ˇ ren´ı 2 Scintigraficke´ vysetˇ Pˇri vniknut´ı nˇekter´ych farmak do organismu doch´az´ı k jejich nahromadˇen´ı na konkr´etn´ıch m´ıstech tk´anˇe nebo org´an˚u, takˇze pokud by byla dotyˇcn´a l´atka oznaˇckov´ana radionuklidem, m˚uzˇ eme toho vyuˇz´ıt k diagnostice [4]. Tˇemto l´atk´am se ˇr´ık´a radiofarmaka. Mˇeˇren´ı prob´ıh´a typicky na zˇ iv´em organismu (mˇeˇren´ı in vivo) v relativnˇe hluboko uloˇzen´ych org´anech, proto mus´ı dan´y radionuklid vyzaˇrovat dostateˇcnˇe tvrd´e z´aˇren´ı gamma. Pot´e lze tˇr´ırozmˇernou distribuci radionuklidu v organismu transformovat na dvourozmˇern´y obraz ve stupn´ıch sˇedi. Sign´alem pˇri vyˇsetˇren´ı je tedy elektromagnetick´e z´aˇren´ı gamma, ostatn´ı doprovodn´a z´aˇren´ı jsou pro n´as neˇza´ douc´ı.
2.1 Scintigrafie Scintigrafie1 je pops´ana v [5] takto: Scintigrafie je fyzik´alnˇe-elektronick´a metoda zobrazen´ı distribuce radioindik´atoru v organismu na z´akladˇe zevn´ı detekce vych´azej´ıc´ıho z´arˇen´ı gama. Do metabolismu cˇ i krevn´ıho obˇehu aplikujeme chemickou l´atku s nav´azan´ym radionuklidem. Ta se pot´e rozloˇz´ı v organismu podle farmakokinetiky dan´eho radiofarmaka. Jak moc se v konkr´etn´ım m´ıstˇe l´atka koncentruje, z´aleˇz´ı na mnoha faktorech, pˇredevˇs´ım na intenzitˇe metabolick´ych a funkˇcn´ıch dˇej˚u v org´anech a tk´an´ıch. Scintigrafii m˚uzˇ eme rozdˇelit na dva druhy [6]: • Statick´a: danou oblast z´ajmu sejmeme jednou cˇ i nˇekolikr´at z r˚uzn´ych u´ hl˚u, nez´aleˇz´ı tedy na cˇ ase (obdoba fotografie). • Dynamick´a: pomoc´ı radioindik´atoru sledujeme dˇej v organismu v cˇ ase, vznik´a tak sekvence statick´ych sn´ımk˚u (obdoba videa). To n´am umoˇznˇ uje dˇelat matematickou anal´yzu mˇeˇren´ı a sledovat funkce org´an˚u. Poznamenejme, zˇ e n´as bude v t´eto pr´aci zaj´ımat pˇredevˇs´ım plan´arn´ı scintigrafie, tedy dvojrozmˇern´e zobrazen´ı (naproti tomu tomografick´a scintigrafie, kde z´ısk´av´ame trojrozmˇern´e zobrazen´ı). C´ılem scintigrafick´eho vyˇsetˇren´ı je kvantitativn´ı zobrazen´ı distribuce radiofarmak v tˇele pacienta (viz obr´azek 2.1). Na rozd´ıl od radiodiagnostick´ych metod zobrazuj´ıc´ıch anatomii vyˇsetˇrovan´e cˇ a´ sti tˇela, scintigrafie poskytuje diagnostick´e informace pˇredevˇs´ım o funkci org´an˚u a tk´an´ı.
1
Pˇresnˇeji sp´ısˇe gamagrafie
2
Obr´azek 2.1: Sekvence sn´ımk˚u a jejich anal´yza (pˇrevzato z [2])
ˇ vysetˇ ˇ ren´ı 2.1.1 Prub ˚ eh N´azornou uk´azku, jak scintigrafie prob´ıh´a, vid´ıme na obr. 2.2. V t´eto kapitole projdeme postupnˇe jednotliv´e kroky zn´azornˇen´e na tomto n´akresu a vysvˇetl´ıme na nich z´akladn´ı fakta, principy a postupy. Na obr. 2.3 pak vid´ıme cˇ asov´y pr˚ubˇeh vyˇsetˇren´ı pomoc´ı dynamick´e scintigrafie.
2.1.2 Radionuklidy Radionuklidem naz´yv´ame prvek, jehoˇz atomy maj´ı schopnost se spont´alnˇe mˇenit na atomy jin´ych prvk˚u (radioaktivn´ı rozpad). Rychlost rozpadu r˚uzn´ych prvk˚u je r˚uzn´a, pro kaˇzd´y prvek je vˇsak konstantn´ı. Radionuklidy jsou v nukle´arn´ı medic´ınˇe nav´azany na nosnou l´atka, farmakum, dohromady tvoˇr´ı radiofarmaka. Neradioaktivn´ı cˇ a´ st radiofarmaka umoˇznˇ uje sledovat vyˇsetˇrovanou funkci tk´anˇe nebo m´a l´ecˇ ebn´y u´ cˇ inek. Radionuklid umoˇznˇ uje sledovat chov´an´ı oznaˇcen´eho farmaka v organismu, pˇr´ıpadnˇe m˚uzˇ e m´ıt tak´e vlastn´ı l´ecˇ ebn´y efekt (terapie radionuklidy). Pro znaˇcen´ı diagnostick´ych pˇr´ıpravk˚u se pouˇz´ıvaj´ı zdroje z´aˇren´ı gama.
3
Obr´azek 2.2: Sch´ema scintigrafick´eho vyˇsetˇren´ı (pˇrevzato z [5])
´ Zakon rozpadu Poˇcet rozpadl´ych jader za jednotku cˇ asu je pˇr´ımo u´ mˇern´y poˇctu dosud nerozpadl´ych jader, tedy: −
dN = λN , dt
(2.1)
ˇ sen´ım t´eto rovnice dostaneme kde N je poˇcet nerozpadl´ych jader a λ je rozpadov´a konstanta. Reˇ tvar rozpadov´eho z´akona jako: Nt = N0 e−λt ,
(2.2)
kde Nt je poˇcet zb´yvaj´ıc´ıch (nerozpadl´ych) jader radionuklidu a N0 je poˇcet jader radionuklidu v cˇ ase t = 0. Poloˇcasem rozpadu pak oznaˇc´ıme dobu, za kterou se rozpadne pr´avˇe polovina jader l´atky. Statisticky´ pohled Rozpad radionuklidu je spont´aln´ı dˇej, m´a n´ahodn´y charakter. Poˇcet vyz´aˇren´ych cˇ a´ stic kol´ıs´a kolem urˇcit´e hodnoty a m´a urˇcit´y rozptyl. Dˇej je pops´an Poissonov´ym z´akonem, pro velk´y poˇcet vyz´aˇren´ych cˇ a´ stic (v praxi N > 50) lze od Poissonova rozdˇelen´ı pˇrej´ıt ke Gaussovu norm´aln´ımu rozdˇelen´ı. Vliv pozad´ı Pˇri detekci z´aˇren´ı z radionuklidu si mus´ıme uvˇedomit dalˇs´ı vlivy, kter´e naˇse mˇeˇren´ı zkresluj´ı. Nepˇr´ıznivˇe ho ovlivˇnuje napˇr´ıklad radiace v mˇeˇr´ıc´ı m´ıstnosti, kosmick´e z´aˇren´ı, stopy radionuklid˚u na detektoru, okoln´ı tk´anˇ atd. Toho si pˇri zpracov´an´ı sn´ımk˚u vytvoˇren´ych pomoc´ı scintigrafie mus´ıme b´yt vˇedomi a neopomenout tuto skuteˇcnost zahrnout do matematick´eho modelu.
4
ˇ Obr´azek 2.3: Casov´ ych pr˚ubˇeh scintigrafick´eho vyˇsetˇren´ı (pˇrevzato z [7])
Poˇzadavky na vlastnosti radionuklidu˚ Z uveden´eho je zˇrejm´e, zˇ e pouˇziteln´e pro n´asˇ u´ cˇ el budou jen nˇekter´e radionuklidy se specifick´ymi vlastnostmi, pˇredevˇs´ım: • Poloˇcas rozpadu mus´ı b´yt dostateˇcnˇe dlouh´y, aby byla l´atka aktivn´ı bˇehem cel´e doby vyˇsetˇren´ı, na druh´e stranˇe je zbytkov´a radiace po ukonˇcen´ı vyˇsetˇren´ı v organismu znaˇcnˇe neˇza´ douc´ı. • Emise energie radionuklidu by se mˇely pohybovat pouze v uˇziteˇcn´ych hodnot´ach, obvykle v rozmez´ı 60keV aˇz ˇra´ dovˇe stovky keV (nejvyˇssˇ´ı pouˇz´ıvan´e z´aˇren´ı je 511keV ). Jak´ekoliv jin´e z´aˇren´ı si nepˇrejeme, protoˇze zbyteˇcnˇe zvyˇsuje radiaˇcn´ı z´atˇezˇ pacienta. • Org´an, kter´y chceme vyˇsetˇrit, by mˇel selektivnˇe zachycovat nosnou l´atku s radionuklidem. Toto zachycen´ı v z´ajmov´e oblasti by mˇelo probˇehnout co nejrychleji. • Vzhledem k aplikovan´emu mnoˇzstv´ı mus´ı b´yt zajiˇstˇena netoxiˇcnost. Protoˇze v´ysledn´e koncentrace radiofarmak ve tk´an´ıch jsou nano aˇz pikomol´arn´ı toxicita vˇetˇsinou nehroz´ı. ˇ Nekter e´ pouˇz´ıvane´ radionuklidy Prvn´ım radionuklidem uˇzit´ym v klinick´e medic´ınˇe byl radioj´od 131 J s poloˇcasem rozpadu 8 dn˚u a energi´ı γ 364keV . Jeho kl´ıcˇ ov´y v´yznam je pˇredevˇs´ım v diagnostice a l´ecˇ bˇe sˇt´ıtn´e zˇ l´azy. Nejd˚uleˇzitˇejˇs´ım radionuklidem pro nukle´arn´ı medic´ınu je technicium 99m Tc s poloˇcasem rozpadu 6 hodin a energi´ı γ 140keV [8]. 99m Tc splˇnuje skoro vˇsechny z´akladn´ı poˇzadovan´e vlastnosti pro scintigrafii a byl proto impulzem pro dalˇs´ı rozvoj nukle´arn´ı medic´ıny. Z 99m Tc
5
nav´ıc vznik´a 99 Tc s poloˇcasem rozpadu ≈ 2 · 105 let, takˇze jej lze povaˇzovat za prakticky stabiln´ı. Energie z´aˇren´ı je ide´aln´ı pro clonˇen´ı kolim´atorem i pro detekci a snadno se z´ısk´av´a ve formˇe aniontu technecistanu 99m TcO− y se d´ale dobˇre v´azˇ e na biologick´e l´atky. 4 , kter´ 201 Jako dalˇs´ı lze uv´est: thalium Tl (perfuze myokardu), galium 67 Ga (scintigrafie n´ador˚u a z´anˇetliv´ych loˇzisek) a krypton 81m Kr (ventilaˇcn´ı scintigrafie plic).
ˇ ı kamera 2.1.3 Scintilacn´ K popisu scintilaˇcn´ı kamery opˇet uˇzijeme [5]: Scintilaˇcn´ı kamera je pˇr´ıstroj, kter´y sn´ım´a fotony z´arˇen´ı γ souˇcasnˇe z cel´eho zorn´eho pole, pˇrev´ad´ı je na elektrick´e impulsy a pomoc´ı nich pak na displeji vytv´arˇ´ı scintigrafick´y obraz distribuce radioindik´atoru v tomto zorn´em poli. Jedn´a se tedy o znaˇcnˇe sloˇzit´e zaˇr´ızen´ı jak konstrukc´ı, tak principem. Naˇs´ım u´ kolem nen´ı scintilaˇcn´ı kameru dokonale popsat, sp´ısˇe n´am p˚ujde o pochopen´ı, jak funguje. Ve vyˇsetˇrovan´em objektu se pohybuj´ı radionuklidy, kter´e izotropnˇe vyzaˇruj´ı z´aˇren´ı gama. Abychom zachytili cˇ a´ stice let´ıc´ı pouze v jednom smˇeru a z´ıskali tak dvourozmˇernou projekci, vloˇz´ıme z´aˇren´ı do cesty olovˇenou desku (nˇekdy wolframovou) s matic´ı otvor˚u. Ta odst´ın´ı cˇ a´ stice, kter´e nejdou pˇresnˇe ve smˇeru osy otvor˚u. Za olovˇenou matic´ı je pak velkoploˇsn´y scintilaˇcn´ı krystal, kter´y pˇri dopadu fotonu vyvol´a v dan´em m´ıstˇe z´ablesk, kter´y je sn´ım´am a pˇrevedem na elektrick´y impulz ve foton´asobiˇci. Pro dalˇs´ı zpracov´an´ı je pak d˚uleˇzit´a digitalizace scintigrafick´ych sn´ımk˚u a jejich uloˇzen´ı do pamˇeti poˇc´ıtaˇce.
2.2 Biologicke´ pˇredpoklady Rozeberme z´akladn´ı vlastnosti mˇeˇren´ych a odhadovan´ych kˇrivek. V prvn´ı ˇradˇe si mus´ıme uvˇedomit, jak´e faktory sn´ım´ame scintigrafickou kamerou pˇri vyˇsetˇren´ı ledvin. Ledvina se skl´ad´a ze dvou hlavn´ıch cˇ a´ st´ı, parenchymu a p´anviˇcky (pelvis). Tyto dvˇe cˇ a´ sti se pˇrekr´yvaj´ı, a prol´ınaj´ı, pro spr´avnou anal´yzu je vˇsak d˚uleˇzit´e je umˇet oddˇelit. Z´akladn´ı sn´ıman´e cˇ a´ sti ve scintigrafii ledvin jsou: • srdce • krevn´ı pozad´ı • parenchym • p´anviˇcka • tk´anˇ ov´e pozad´ı • dalˇs´ı org´any.
6
2.2.1 Tvary faktorovych ´ kˇrivek Srdce je org´an, k jehoˇz naplnˇen´ı konstrastn´ı l´atkou dojde jako prvn´ı, je-li, jako v naˇsem pˇr´ıpadˇe, aplikov´ano radiofarmakum do krve ([9]). Z krevn´ıho obˇehu je tato l´atka cˇ iˇstˇena ledvinami, teoreticky se d´a oˇcek´avat exponenci´aln´ı tvar kˇrivky krve, jak je uk´az´ano na obr´azku 2.4. y
0
5
10
15
ˇcas [min]
Obr´azek 2.4: Typick´y pr˚ubˇeh kˇrivky krve (srdce a krevn´ı pozad´ı)
Ze srdce se pot´e l´atka krv´ı distribuuje do cel´eho tˇela, nejprve do velk´ych c´ev a pot´e do cel´eho krevn´ıho ˇreˇciˇstˇe. Radiofarmakum je vytv´aˇreno tak, aby bylo selektivnˇe zachyt´av´ano v c´ılov´em org´anu, v ledvinˇe. V poˇca´ tku cel´e sekvence vˇsak m˚uzˇ eme pozorovat i dalˇs´ı org´any, do kter´ych se radiofarmakum dostalo pˇres krev. Tento fakt ztˇezˇ uje celou anal´yzu, pomoci n´am m˚uzˇ e to, zˇ e pr˚ubˇeh kˇrivky v nich je pomˇernˇe rychl´y a tedy je sˇance tyto org´any detekovat a t´ım automaticky odeˇc´ıst z v´ysledn´ych sn´ımk˚u. Dalˇs´ı probl´em s detekc´ı krevn´ıho pozad´ı n´am vyvst´av´a v tom, zˇ e prob´ıh´a pˇrirozen´a difuze mezi krv´ı a tk´anˇemi. T´ım vnik´a dalˇs´ı typ pozad´ı, tzv. tk´anˇ ov´e pozad´ı, kter´e pˇreb´ır´a a zase vypouˇst´ı radiofarmakum z a do krve. Lok´alnˇe se tedy m˚uzˇ e st´at, zˇ e je kˇrivka krve v jedn´e chv´ıli m´ırnˇe rostouc´ı, glob´alnˇe by vˇsak mˇelo doch´azet k poklesu aktivity v krvi. V dynamick´e scintigrafii ledvin doch´az´ı k nejvˇetˇs´ı selekci radiofarmaka v ledvin´ach, zamˇeˇr´ıme se tedy nyn´ı na tento org´an. Ledvinu mus´ıme rozdˇelit na dvˇe hlavn´ı cˇ a´ sti, na parenchym a p´anviˇcku. Ty jsou navz´ajem propojeny a kaˇzd´a m´a rozd´ılnou dynamiku. Tak´e jejich oddˇelen´ı je jedn´ım z hlavn´ıch u´ kol˚u anal´yzy. Parenchym tvoˇr´ı vˇetˇs´ı cˇ a´ st ledviny a v nˇem doch´az´ı k zachyt´av´an´ı kontrastn´ı l´atky z krve. P´anviˇcka je menˇs´ı a k jej´ımu plnˇen´ı a z naˇseho pohledu tedy k jej´ı aktivaci doch´az´ı pozdˇeji, nebot’ do n´ı pˇrit´ek´a tekutina z parenchymu. Tento bod nast´av´a typicky po 2 aˇz 3 minut´ach, v z´avislosti na vˇeku a pouˇzit´e kontrastn´ı l´atce ([10]). P´anviˇcka je vˇsak jiˇz urˇcitˇe aktivovan´a ve
7
y
parenchym
pelvis
2-3
5
10
15
ˇcas [min]
Obr´azek 2.5: Typick´y pr˚ubˇeh kˇrivky parenchymu a p´anviˇcky
chv´ıli, kdy m´a kˇrivka parenchymu sv´e maximum. V tomto bodˇe je pˇr´ıtok do p´anviˇcky a pˇr´ıtok do parenchymu vyrovn´an. Typick´e pr˚ubˇehy parenchymu a p´anviˇcky u zdrav´e ledviny vid´ıme na obr´azku 2.5. Kaˇzd´a kˇrivka faktoru (x) je d´ana v´ysledkem konvoluce sv´eho konvoluˇcn´ıho j´adra (u) s kˇrivkou krve (b) ([11], [12], [13]). Dohromady tedy: xt =
t X
bt−m+1 um .
(2.3)
m=1
Tato konvoluˇcn´ı j´adra (naz´yvan´a v nukle´arn´ı medic´ınˇe impulzn´ı retenˇcn´ı funkce) maj´ı typick´y teoretick´y tvar, kter´y je zn´azornˇen na obr´azku 2.6. Pakliˇze je faktor aktivn´ı hned od poˇca´ tku sekvence, zaˇc´ın´a konstantn´ı plato od nuly, pokud tomu tak nen´ı (napˇr. p´anviˇcka), je jeho zaˇca´ tek posunut.
´ 2.2.2 Faktorove´ obrazky Faktorov´e obr´azky z´ıskan´e ve scintigrafii pˇredstavuj´ı ve vˇetˇsinˇe pˇr´ıpad˚u nˇejak´e org´any, pˇr´ıpadnˇe jejich cˇ a´ sti. Typicky jsou tedy kompaktn´ı a omezen´e. D´ıky pˇredpoklad˚um: • sn´ıman´a tk´anˇ se vzhledem ke scintigrafick´e kameˇre neh´ybe • org´any nemˇen´ı sv˚uj tvar bˇehem mˇeˇren´ı • zmˇena objemu kapaliny uvnitˇr org´anu je line´arnˇe u´ mˇern´a poˇctu vyz´aˇrovan´ych cˇ a´ stic
8
y krev
parenchym
pelvis tk´an ˇ ov´e pozad´ı 2-3
5
10
15
ˇcas [min]
Obr´azek 2.6: Teoretick´y tvar konvoluˇcn´ıch jader jednotliv´ych faktor˚u (u zdrav´eho pacienta)
m˚uzˇ eme namˇeˇren´a data modelovat jako line´arn´ı kombinaci jednotliv´ych org´an˚u (nebo jejich cˇ a´ st´ı). Naˇs´ım z´akladn´ım u´ kolem je pak zjistit tvar jednotliv´ych org´an˚u a pr˚ubˇeh pr˚utoku kontrastn´ı tekutiny v nich (viz obr´azek 2.1). Kaˇzd´y sn´ımek sekvence (index t) lze matematicky vyj´adˇrit jako dt =
r X
at xt + et ,
(2.4)
i=1
kde r je poˇcet faktor˚u, at je faktorov´y obr´azek, xt je faktorov´a kˇrivka a et pˇredstavuje sˇum. Tyto pˇredpoklady by mˇel respektovat i matematick´y model, kter´y spr´avnˇe odhadne jednotliv´e faktorov´e obr´azky a jejich kˇrivky.
ˇ ´ ı analyzy 2.3 Soucasn e´ metody funkcionaln´ ´ obrazove´ sekvence Pod´ıvejme se nyn´ı na cel´y probl´em anal´yzy obrazov´e sekvence z matematick´eho hlediska. Pokus´ıme se popsat z´akladn´ı biologick´e pˇredpoklady, kter´e na jednotliv´e v´ysledn´e i pomocn´e kˇrivky klademe. To bude m´ıt vliv na n´aslednout faktorovou anal´yzu a z´ısk´an´ı diagnostick´ych
9
koeficient˚u. Dodejme jeˇstˇe, zˇ e v souˇcasn´e dobˇe prov´ad´ı anal´yzu expert (l´ekaˇr cˇ i laborant) pomoc´ı vizu´aln´ıho dojmu, zda dan´y faktor obsahuje nebo neobsahuje konkr´etn´ı pixel je zcela na nˇem. Obdobn´e je to napˇr. s urˇcen´ım, kde dan´a kˇrivka zaˇc´ın´a a konˇc´ı (viz napˇr. obr´azky 2.10 a 2.11). Vytvoˇren´ı a vyˇreˇsen´ı matematick´eho modelu, kter´y vysvˇetluje data, sn´ıman´a scintigrafickou kamerou, jeˇstˇe nemus´ı znamenat, zˇ e jsme nalezli biologicky spr´avn´e ˇreˇsen´ı. V´ysledek m˚uzˇ e b´yt velice vzd´alen biologick´ym pˇredpoklad˚um. Jiˇz sestaven´ım modelu totiˇz vtisk´av´ame pˇr´ıpadn´emu ˇreˇsen´ı n´ami danou podobu, kter´a m˚uzˇ e b´yt v podstatˇe libovoln´a. Jako pˇr´ıklad uved’me n´asleduj´ıc´ı: mˇejme za u´ kol rozloˇzit cˇ´ıslo d na a · x. N´asˇ matematick´y model tedy bude m´ıt tvar d = ax. Ovˇsem pokud nespecifikujeme vlastnosti tohoto rozkladu, pak je moˇzn´ych ˇreˇsen´ı nekoneˇcnˇe mnoho, protoˇze napˇr. 12 = 3 · 4, ale i 2 · 6 atd. (pro dalˇs´ı viz [2], [14]). Velk´y vliv na ˇreˇsen´ı cel´eho probl´emu m´a tedy jiˇz konstrukce modelu, ve kter´em m˚uzˇ eme ˇr´ıci, jak si ”pˇrejeme” data rozloˇzit. S ohledem na to mus´ıme konstruovat naˇse modely i v n´asleduj´ıc´ıch kapitol´ach.
´ 2.3.1 Oblasti zajmu Jedn´ım ze z´akladn´ıch krok˚u pˇri anal´yze obrazov´ych sekvenc´ı vznikl´ych scintigrafi´ı je definov´an´ı tzv. oblast´ı z´ajm˚u (region of interest - ROI). Ty pot´e slouˇz´ı k mˇeˇren´ı aktivity v oznaˇcen´e oblasti, v´ysledkem je cˇ asov´a kˇrivka, kter´e se vyuˇz´ıv´a k v´ypoˇctu kvantitativn´ıch diagnostick´ych parametr˚u. Probl´em s t´ım spojen´y vysvˇetl´ıme na pˇr´ıkladu. Na obr´azku 2.7 vid´ıme sn´ımek dutiny bˇriˇsn´ı. My dok´azˇ eme pouh´ym pohledem rozliˇsit tˇri z´akladn´ı struktury, kter´e se na obr´azku nach´azej´ı, konkr´etnˇe dvˇe ledviny a v horn´ı cˇ a´ sti srdce. Poˇc´ıtaˇc ovˇsem tyto struktury tak snadno nedetekuje, v klinick´e praxi jsou proto skoro vˇzdy urˇcov´any uˇzivatelem, tedy l´ekaˇrem.
Obr´azek 2.7: Sn´ımek dutiny bˇriˇsn´ı
10
Probl´em s jak´ymkoliv implicitn´ım nastaven´ım je ilustrov´an na obr´azc´ıch 2.8 a 2.9. Okˇr´ıdlen´e tvrzen´ı, zˇ e ”neexistuj´ı dva stejn´ı lid´e” plat´ı i v tomto pˇr´ıpadˇe. Pokud by sˇlo pˇredpokl´adat, zˇ e rozm´ıstˇen´ı org´an˚u bude vˇzdy stejn´e, mohli bychom je hledat v urˇcit´ych typick´ych oblastech (viz obr. 2.8). Jak vˇsak vid´ıme, v tomto pˇr´ıpadˇe jsou ledviny rozm´ıstˇeny jinak a lev´a ledvina nen´ı ve sv´e ROI zahrnuta cel´a, coˇz by mohlo v´est ke znaˇcnˇe zkreslen´ym v´ysledk˚um. Proto je nutno interakce s l´ekaˇrem, kter´y vyznaˇc´ı ROI pˇresnˇe dle konkr´etn´ıho pacienta (obr. 2.9).
Obr´azek 2.8: Z´akladn´ı nastaven´ı
Obr´azek 2.9: Konkr´etn´ı ROI
Uvˇedomme si, zˇ e urˇcov´an´ı ROI je slab´ym cˇ l´ankem cel´eho vyˇsetˇren´ı, nebot’ interakce s uˇzivatelem zan´asˇ´ı t´emˇeˇr vˇzdy nepˇresnost. Klasick´e metody zpracov´an´ı vˇsak ROI vyuˇz´ıvaj´ı a jsou tedy do znaˇcn´e m´ıry z´avisl´e i na zkuˇsenostech a zruˇcnosti obsluhy.
2.3.2 Matematicky´ pohled a vyhodnocen´ı kˇrivek Podaˇr´ı-li se n´am z´ıskat v´ysledn´e tvary kˇrivek jednotliv´ych faktor˚u, je tˇreba je d´ale analyzovat. Uved’me nˇekter´e z´akladn´ı pˇr´ıstupy. Patlak-Rutlanduv ˚ graf Prvn´ım moˇzn´ym pˇr´ıstupem je sestrojen´ı tzv. Patlak-Rutlandova grafu [15]. Tato metoda pˇredpokl´ad´a, zˇ e kontrastn´ı l´atka do org´an˚u vstupuje a setrv´av´a v nich. Org´any jsou tedy nˇeco jako integr´atory vstupn´ı aktivity. Je vˇsak nutn´e urˇcit polohu org´anu, pˇr´ıpadnˇe oblast z´ajmu, ve kter´e org´an leˇz´ı. Necht’ poˇcet dopad˚u cˇ a´ stic v oblasti z´ajmu je d´an funkc´ı R(t) a aktivita v krvi funkc´ı B(t). Potom Z R(t) = a B(t) dt + bB(t), (2.5)
kde konstanta a urˇcuje schopnost absorbce kontrastn´ı l´atky org´anem a konstanta b vyjadˇruje vliv pozad´ı.
11
Tento model analyzujeme vydˇelen´ım cel´e rovnice funkc´ı krve B(t) a dost´av´ame R B(t) dt R(t) =a + b, B(t) B(t) R
(2.6)
B(t) dt
R(t) coˇz je graf z´avislosti B(t) na B(t) . V´ysledkem je tedy pˇr´ımka se smˇernic´ı a a vertik´aln´ım posunem b, ze kter´e vyˇcteme potˇrebn´e informace. ´ Uskal´ ı tohoto pˇr´ıstupu je v tom, jak pˇresnˇe zvl´adneme namˇeˇrit cˇ i vypoˇc´ıtat zm´ınˇen´e funkce. V´ysledek je tak´e z´avisl´y na zruˇcnosti odborn´e obsluhy [16], jak vid´ıme na obr´azc´ıch 2.10 a 2.11. Volba d´elky line´arn´ıho u´ seku je cˇ ast´a chyba lidsk´eho faktoru. Z´aroveˇn si m˚uzˇ eme vˇsimnout, zˇ e jsme cel´a data transformovali na jedinou kˇrivku, cˇ´ımˇz doch´az´ı ke ztr´atˇe cˇ a´ sti informace. Patlak−Rutland plot
Patlak−Rutland plot
0.5
3 2.5
0
−0.5
R(t)/P(t)
R(t)/P(t)
2
−1
1.5 1 0.5
−1.5
0 −2 0
1
2
3
4
5
pseudo−time ∫ P(t)/P(t)
6
−0.5 0
7
Obr´azek 2.10: Ihned zˇrejm´y tvar PatlakRutlandova grafu (z [16] )
1
2
3
4
5
pseudo−time ∫ P(t)/P(t)
6
Obr´azek 2.11: Problematick´y tvar PatlakRutlandova grafu (z [16] )
Dekonvoluce Druh´y pˇr´ıstup je zaloˇzen na pˇredpokladu, zˇ e kˇrivka kaˇzd´eho faktoru je modelov´ana jako konvoluce kˇrivky krve a konvoluˇcn´ıho j´adra kaˇzd´eho faktoru ([11], [10], [15], [12]), jak bylo diskutov´ano jiˇz v´ysˇe. Naˇs´ım c´ılem je urˇcit impulzn´ı retenˇcn´ı funkci H(t). V´yjimeˇcnˇe vˇsak m˚uzˇ eme aplikovat radiofarmakum pˇr´ımo do vyˇsetˇrovan´eho org´anu, cˇ astˇeji je aplikov´ano do krve a z nˇej teprve pˇrech´az´ı d´ale. Mnoˇzstv´ı kontrastn´ı l´atky v org´anu, Q(t) je tedy d´ano konvoluc´ı H(t) a vstupn´ı funkc´ı org´anu, I(t), coˇz m˚uzˇ eme zapsat jako Q(t) =
ZT
I(T − t)H(t) dt
(2.7)
0
a mˇeˇrit hodnoty I(t) a Q(t). Vyˇreˇsen´ım uveden´e konvoluˇcn´ı rovnice pak dostaneme hodnoty H(t), ze kter´ych opˇet vyˇcteme potˇrebn´e informace.
12
Tento postup je v principu jednoduch´y, v praxi vˇsak naˇra´ zˇ´ı na mnoh´e probl´emy, pˇredevˇs´ım urˇcen´ı spr´avn´e vstupn´ı kˇrivky a na sˇum. Poznamenejme, zˇ e porovn´an´ı Patlak-Rutlandova grafu a dekonvoluce nalezneme v [13]. Oba jsou do znaˇcn´e m´ıry z´avisl´e na spr´avn´em urˇcen´ı a odeˇcten´ı sˇumu (kter´y tvoˇr´ı samotn´y princip metody), krevn´ıho a tk´anˇ ov´eho pozad´ı atd. Moˇznost´ı, jak tento sˇum odstranit je to, zˇ e vezmeme nˇejak´e referenˇcn´ı pozad´ı a to pak od cel´eho org´anu odeˇcteme. Nar´azˇ´ıme t´ım vˇsak na probl´em, jak takov´e referenˇcn´ı pozad´ı volit a zda t´ım neodeˇc´ıt´ame moc velkou nebo malou hodnotu. Jak je uk´az´ano v [1], rozd´ıl v´ysledku oproti skuteˇcnosti m˚uzˇ e b´yt aˇz 25% v z´avislosti na volbˇe referenˇcn´ıho pozad´ı, coˇz je enormnˇe velk´e cˇ´ıslo. Analyza ´ hlavn´ıch komponent Anal´yza hlavn´ıch komponent (t´ezˇ Principal Component Analysis, PCA) je obecnˇe postup slouˇz´ıc´ı k dekorelaci dat. D˚uleˇzit´ym pˇredpokladem t´eto metody je, zˇ e pozorujeme superpozici jednotliv´ych faktor˚u. To je pˇresnˇe n´asˇ pˇr´ıpad, nebot’ ve scintigrafii sledujeme transformaci tˇr´ırozmˇern´eho prostoru na dvourozmˇern´y sn´ımek, kter´y je line´arn´ı kombinac´ı jednotliv´ych faktor˚u. Z´akladn´ı princip je n´asleduj´ıc´ı: mˇejme datovou (v´ystup experimentu) matici D a tu rozloˇzme jako D = P ′ΛP,
(2.8)
kde Λ je diagon´aln´ı matice vlastn´ıch cˇ´ısel matice D a P je matice vlastn´ıch vektor˚u pˇr´ısluˇsn´ych k matici Λ, pˇriˇcemˇz plat´ı, zˇ e P ′ P je rovno jednotkov´e matici. Prvn´ıch r nejvˇetˇs´ıch vlastn´ıch cˇ´ısel a k nim pˇr´ısluˇsn´ych vlastn´ıch vektor˚u jsou nejv´yznamˇejˇs´ı faktory v datech D. Obecnˇejˇs´ı pˇr´ıpad pak ˇreˇs´ı tzv. Singular value decomposition. Nepˇr´ıjemnost´ı pˇri tomto rozkladu je opˇet sˇum. Pouˇzit´ı t´eto techniky je demonstrov´ano napˇr. v [17]. Model faktorove´ analyzy ´ Faktorov´a anal´yza vych´az´ı z pˇredchoz´ı podkapitoly (PCA) a z kapitoly 2.2.2. Datov´a matice D zde obsahuje ve sloupc´ıch vˇsechny sn´ımky sekvence. D se pak snaˇz´ıme rozloˇzit na matice A a X, kde A obsahuje ve sloupc´ıch jednotliv´e faktorov´e obr´azky a X obsahuje ve sloupc´ıch jednotliv´e kˇrivky. Model lze zapsat jako: D = AX ′ + E,
(2.9)
kde E pˇredstavuje sˇum. Poznamenejme, zˇ e toto bude z´akladn´ı model dat pro celou tuto pr´aci.
2.3.3 Klinicke´ vyhodnocen´ı K urˇcen´ı spr´avn´e diagn´ozy pacienta by vzhledem k v´ysˇe popsan´ym pˇr´ıstup˚um bylo zapotˇreb´ı matematika - l´ekaˇre. To je v praxi n´aroˇcn´y poˇzadavek, proto je potˇreba v´ysledek anal´yzy kvantifikovat a transformovat ide´alnˇe na jedno cˇ´ıslo, jehoˇz hodnota by byla dosateˇcnˇe vypov´ıdaj´ıc´ı
13
o celkov´em stavu. Z tohoto d˚uvodu se zav´ad´ı tzv. sekund´arn´ı (diagnostick´e) parametry, kter´e z´ısk´ame pomoc´ı metod popsan´ych v pˇredchoz´ıch kapitol´ach. Clearance Clearance je imagin´arn´ı objem tekutiny, kter´a je od kontrastn´ı l´atky kompletnˇe oˇciˇstˇena za jednotku cˇ asu [18]. Typicky jsou tedy jednotky clearance ml/min. Clearance m˚uzˇ e b´yt definov´ana t´ezˇ relativnˇe jako procento tekutiny, kter´a je vyˇciˇstˇena pˇri jednom pr˚uchodu org´anem. Form´alnˇe bychom clearanci (C) vypoˇc´ıtali jako C=
cmoˇc Vmoˇc , cplasma
(2.10)
kde cmoˇc je koncentrace sledovan´e l´atky v moˇci, cplasma koncentrace sledovan´e l´atky v plasmˇe a Vmoˇc je objem vylouˇcen´e moˇci. Tento jednoduch´y postup bohuˇzel vyuˇz´ıt nem˚uzˇ eme, protoˇze sledujeme pacienta jako celek a ne jeho jednotliv´e cˇ a´ sti zvl´asˇt’. Pomoc´ı glob´aln´ıch parametr˚u tak nejsme schopni popsat fuknci napˇr. konkr´etn´ı ledviny. Z toho d˚uvodu mus´ıme sledovat sekund´arn´ı znaky a z nich odhadnout hodnotu clearence. V t´eto pr´aci se pokus´ıme namˇeˇrit tzv. relativn´ı ren´aln´ı clearance, kterou vypoˇc´ıt´ame z v´ysledn´ych intenzit faktor˚u jako relL =
L × 100, R+L
(2.11)
kde L a R jsou intenzity mˇeˇren´e v lev´em a v prav´em parenchymu. ˇ Tranzitn´ı cas Clearance nen´ı samozˇrejmˇe jedin´y pˇr´ıznak, ze kter´eho se n´aslednˇe urˇcuje diagn´oza ledvin. Jako pˇr´ıklad dalˇs´ıho parametru uved’me tranzitn´ı cˇ as (TT) [10], jehoˇz hodnota jde snadno vyˇc´ıst z impulzn´ı retenˇcn´ı funkce, podaˇr´ı-li se n´am z´ıskat napˇr. dekonvoluc´ı. TT pak m˚uzˇ eme mˇeˇrit nˇekolik typ˚u, jak je naznaˇceno na obr´azku 2.12. Metoda ovˇsem pˇredpokl´ad´a, zˇ e dok´azˇ eme pˇresnˇe urˇcit hranice pro urˇcen´ı TT, coˇz m˚uzˇ e b´yt v praxi probl´em, jak jiˇz bylo uk´az´ano.
2.4 Zhodnocen´ı V t´eto kapitole jsme popsali cel´y postup scintigrafick´e anal´yzy. Uk´azali jsme pozad´ı, v´yhody a problematick´e body t´eto diagnostick´e metody a n´aslednˇe se sousteˇredili na funcion´aln´ı anal´yzu namˇeˇren´ych sn´ımk˚u. Popsali jsme z´akladn´ı v souˇcasnosti vyuˇz´ıvan´e metody anal´yzy. Hlavn´ı nev´yhodou scintigrafie je nedokonal´e rozliˇsen´ı, aˇz o ˇra´ d horˇs´ı neˇz napˇr. u CT. Dalˇs´ım probl´emem je superpozize org´an˚u, kdy doch´az´ı k souˇctu z´aˇren´ı ze vˇsech vrstev a t´ım m˚uzˇ e doj´ıt ke zkreslen´ı informace. Diskutov´an jiˇz byl i sˇum. Na druhou stranu n´am scintigrafie umoˇznˇ uje zkoumat procesy, ke kter´ym doch´az´ı v zˇ iv´ych organismech, kter´e nejsou jin´e diagnostick´e techniky schopny rozpoznat. Ide´aln´ı je pak hybridn´ı zobrazen´ı, kter´e kombinuje
14
y
ˇcas [min]
V TT RT Tmin RT Tmean
RT Tmax Obr´azek 2.12: Tranzitn´ı cˇ as: vaskul´arn´ı TT (VTT), ren´aln´ı TT (RTT)
metodu s vysok´ym rozliˇsen´ım a scintigrafii, takˇze m´ame kvalitn´ı obraz a pˇritom m˚uzˇ eme zachytit i prob´ıhaj´ıc´ı procesy. N´aslednˇe se jako problematick´e ukazuje, zˇ e jednotliv´e metody anal´yzy sn´ımk˚u funguj´ı oddˇelenˇe, nikdy ne jako synt´eza. Doch´az´ı tak k podstatn´e ztr´atˇe informace v pr˚ubˇehu v´ypoˇctu, kter´a se pot´e jiˇz nem´a kde kompenzovat. Velk´e chyby zan´asˇ´ı do cel´e metody t´ezˇ lidsk´y faktor. Jeˇstˇe vˇetˇs´ı neurˇcitost dost´av´a do probl´emu sˇum, kter´y je mnohy velice tˇezˇ k´e oddˇelit od uˇziteˇcn´eho sign´alu. Ide´aln´ı by tedy bylo zkonstruovat model, kter´y bude kombinovat v´yhody pˇredchoz´ıch pˇr´ıstup˚u a bude fungovat automaticky, tedy bez nutnosti neust´al´e pˇr´ıtomnosti odborn´e obsluhy.
15
3 Bayesovska´ teorie a aproximace V t´eto kapitole se budeme vˇenovat matematick´emu apar´atu, kter´y pak budeme vyuˇz´ıvat v n´asleduj´ıc´ıch kapitol´ach. Bude se jednat pˇredevˇs´ım o cˇ a´ sti z teorie pravdˇepodobnosti, kter´e se vztahuj´ı k Bayesovsk´e statistice. Ta n´am v n´asleduj´ıc´ıch kapitol´ach efektivnˇe umoˇzn´ı zahrnut´ı informac´ı jako: • expertn´ı informace • pˇredpoklad zˇ a´ dn´e apriorn´ı informace do matematick´eho modelu mnohem sn´aze, neˇz pomoc´ı klasick´e statistiky. Nejprve rozebereme z´aklady Bayesovsk´e statistiky a pot´e pˇrejdeme k metod´am aproximace, pˇredevˇs´ım Bayesovu variaˇcn´ımu teor´emu.
´ 3.1 Zaklady Bayesovske´ teorie Na zaˇca´ tku poznamenejme, zˇ e se budeme zab´yvat pouze spojit´ymi a nikoliv diskr´etn´ımi hustotami pravdˇepodobnost´ı, ty bychom dostali v´ymˇenou integr´al˚u za sumy. Necht’ D jsou namˇeˇren´a data a θ jsou parametry pravdˇepodobnostn´ıho modelu s apriorn´ım rozdˇelen´ım f (θ), θ ∈ Θ, kde Θ je prostor parametr˚u θ, kter´e ovˇsem oproti klasick´e statistice budeme povaˇzovat za promˇenn´e. Hustotu pravdˇepodobnosti D podm´ınˇenou parametry θ pak zap´ısˇeme jako f (D|θ). Sdruˇzenou hustotu pravdˇepodobnosti D a θ urˇc´ıme jako R f (D, θ) = f (θ, D) = f (D|θ)f (θ). Uvˇedomme si tak´e, zˇ e f (D) = Θ f (D|θ)f (θ)dθ. Informaci o parametrech θ, kterou m˚uzˇ eme z dat D z´ıskat, vyj´adˇr´ıme jako aposteriorn´ı rozdˇelen´ı f (θ|D). T´ım se dost´av´ame k Bayesovu teor´emu, kter´y tvar tohoto aposteriorn´ıho rozdˇelen´ı urˇcuje jako: f (θ|D) = f (D) = jako:
R
Θ
f (θ, D) f (D|θ)f (θ) =R . f (D) f (D|θ)f (θ)dθ Θ
(3.1)
f (D|θ)f (θ)dθ pˇredstavuje normalizaˇcn´ı konstantu, Bayes˚uv teor´em pak lze pˇrepsat f (θ|D) ∝ f (D|θ)f (θ) ,
(3.2)
kde symbol ∝ znaˇc´ı v´yraz proporcion´alnˇe rovno, tzn. rovno aˇz na konstantu. V´ypoˇcet moment˚u aposteriorn´ıho rozdˇelen´ı, tj. pˇredpokl´adan´ych hodnot funkc´ı g(θ), pak bude vypadat jako: Z Ef (θ|D) [g(θ)] = g(θ)f (θ|D)dθ . (3.3) Θ
16
d tedy Pro Ef (θ|D) [g(θ)] budeme uˇz´ıvat znaˇcen´ı g(θ),
d = Ef (θ|D) [g(θ)] . g(θ)
(3.4)
ˇ 3.1.1 Volba apriorn´ıho rozdelen´ ı
D˚uleˇzit´ym krokem je vhodn´a volba apriorn´ıho rozdˇelen´ı, protoˇze to pak z´asadn´ım zp˚usobem ovlivˇnuje v´ypoˇcet aposteriorn´ıho rozdˇelen´ı. To je diskutov´ano napˇr. v [2],[19]. Obecnˇe plat´ı, zˇ e pokud o datech pˇredpokl´ad´ame, zˇ e nesou velkou m´ıru informace, vol´ıme vˇetˇsinou apriorn´ı rozdˇelen´ı, kter´e nese minim´aln´ı m´ıru informace, cˇ´ımˇz doc´ıl´ıme toho, zˇ e m´a mal´y vliv na aposteriorn´ı rozdˇelen´ı. To se tedy rekonstruuje pˇredevˇs´ım z dat a my do nˇej nevn´asˇ´ıme svoje pˇredpoklady, kter´e by nemusely b´yt spr´avn´e. Volba apriorna m´a t´ezˇ vliv na sloˇzitost analytick´eho v´ypoˇctu. Jako apriorn´ı rozdˇelen´ı je proto pro n´as v´yhodn´e volit distribuce z exponenci´aln´ı rodiny rozdˇelen´ı, pˇredevˇs´ım konjugovan´a apriorna ([20]). ˇ Definice 1 Rekneme, zˇ e pravdˇepodobnostn´ı funkce f je z rodiny exponenci´aln´ıch rozdˇelen´ı, pokud: fx (x|θ) = h(x) exp(η(θ)T (x) − A(θ)) ,
(3.5)
kde funkce h(x), η(θ), T (x) a A(x) jsou zn´am´e. Konjugovan´e apriorno π pro parametr η z exponenci´aln´ı rodiny je d´ano jako: π(η) ∝ exp(−η ′ α − βA(η)) ,
(3.6)
kde α a β jsou hyperparametry apriorna (parametry, kter´e rˇ´ıd´ı jin´e parametry). Konjugovan´a apriorna n´am pot´e pˇri v´ypoˇctu generuj´ı aposteriorna stejn´eho typu (rodiny), jako jsou tato apriorna, coˇz n´am umoˇzn´ı realizovat v´ypoˇcet pomoc´ı cyklick´eho pˇrepoˇc´ıt´av´an´ı distribuc´ı, jejichˇz tvar z˚ustane zachov´an.
ˇ sen´ ˇ ı pomoc´ı aproximace 3.2 Re V t´eto cˇ a´ sti se budeme vˇenovat hled´an´ı aposteriorn´ıho rozdˇelen´ı (3.1), coˇz se uk´azˇ e pˇri podrobnˇejˇs´ı u´ vaze analyticky obecnˇe prakticky nemoˇzn´e. Pro jeho v´ypoˇcet je nutno aplikovat Bayes˚uv teor´em, spoˇc´ıtat normalizaˇcn´ı konstantu a jeho momenty, coˇz lze jen pro u´ zkou tˇr´ıdu rozdˇelen´ı. Pouˇzit´ı numerick´eho ˇreˇsen´ı by zde pˇri vyˇssˇ´ıch dimenz´ıch tak´e nebylo pˇr´ıliˇs ide´aln´ı. Jako vhodn´e ˇreˇsen´ı tohoto probl´emu se ukazuje aproximace pomoc´ı vhodn´e distribuce f (θ|D) ≈ f˜(θ|D) , kde pro f˜ jsou v´ysˇe uveden´e operace zn´am´e nebo snadno dopoˇc´ıtateln´e. V n´asleduj´ıc´ıch kapitol´ach pop´ısˇeme d˚uleˇzit´e metody aproximace.
17
(3.7)
3.2.1 EM (expectation minimalization) algoritmus Tento algoritmus byl publikov´an v [21]. My pouˇzijeme jeho upravenou verzi, dle [22]. Nejlepˇs´ı aposteriorn´ı odhad (MAP - maximum aposteriory) je definov´an jako: δM AP (θ) = arg max f (θ|D) . θ∈Θ
Pro symetrick´e f stˇredn´ı hodnota parametru θ odpov´ıd´a δM AP (θ), m˚uzˇ eme tedy ps´at: Z b δM AP (θ) = θ = θf (θ|D) dθ .
(3.8)
(3.9)
θ∈Θ
V´ypoˇcet tohoto integr´alu ovˇsem nemus´ı b´yt, obzvl´asˇtˇe pro vysok´e dimenze, realizovateln´y, proto mus´ıme pˇristoupit k n´asleduj´ıc´ı aproximaci. Pˇredpokl´adejme napˇr´ıklad, zˇ e parametr θ m´a dvˇe dimenze, tedy: θ = [θ1 , θ2 ]. Potom margin´aln´ı rozdˇelen´ı θ1 aproximujeme jako f (θ1 |D) ≈ f (θ1 |D, θb2 ) , kde θb2 je pevn´y konkr´etn´ı odhad z´ıskan´y v pˇredchoz´ algoritmu. EM-algoritmus pak m´a tvar [23]: ım kroku (i−1) (i) ˜ b E-krok: f (θ1 ) ≈ f θ1 |D, θ2 R (i) M-krok: θb2 = arg max θ1 f˜(i) (θ1 ) ln f (D, θ1 , θ2 )dθ1 . θ2
3.2.2 Laplaceova aproximace
Laplaceova metoda je uˇziteˇcn´a aproximace pravdˇepodobnostn´ı funkce f (D|θ) v okol´ı MAP odhadu θb pomoc´ı norm´aln´ıho rozdˇelen´ı. Dle [24],[25] m´a f (θ|D) tvar norm´aln´ıho rozdˇelen´ı se stˇredn´ı hodnotou θb a s kovarianˇcn´ı matic´ı −H −1 , tedy: b −H −1 ) , f (θ|D) ≈ N(θ,
kde H je Hessova matice (napˇr. [26]) logaritmu f (D|θ)f (θ): ∂ 2 log f (θ|D) (H)i,j = b. ∂θi ∂θj θ=θ
(3.10)
(3.11)
3.2.3 Markov Chain Monte Carlo (MCMC) aproximace MCMC je strategie pro generov´an´ı hodnot θ(i) , kter´e pokr´yvaj´ı prostor Θ s vyuˇzit´ım Markovova ˇretˇezce, kter´y zajist´ı, zˇ e ˇretˇez tr´av´ı v´ıce cˇ asu v d˚uleˇzit´ych oblastech pravdˇepodobnostn´ıho prostoru (viz [27]). Typick´e pouˇzit´ı je pr´avˇe poˇc´ıt´an´ı v´ıcerozmˇern´ych integr´al˚u. Sekvenci {θ(1) , . . . , θ(i) } nazveme Markovov´ym ˇretˇezcem, pokud f (θ(i) |θ(i−1) , . . . , θ(1) ) = f (θ(i) |θ(i−1) ) , tedy pokud n-t´y stav zavis´ı pouze na tom pˇredchoz´ım.
18
(3.12)
Podle principu Monte Carlo pak Z N 1 X a.s. (i) g(θ ) −−−→ I(g) = g(θ)f (θ|D) dθ . IN (g) = N →∞ N i=1 Θ
(3.13)
V´ypoˇcetn´ı n´aroˇcnost je vˇsak relativnˇe vysok´a, obzvl´asˇtˇe pro vysok´e dimenze θ.
ˇ ı Bayesova aproximace 3.2.4 Variacn´ V t´eto cˇ a´ sti se budeme zab´yvat aproximac´ı pomoc´ı Variaˇcn´ıho Bayesova teor´emu. T´eto metodˇe vˇenujeme v´ıce prostoru neˇz ostatn´ım, protoˇze ji budeme v n´asleduj´ıc´ıch kapitol´ach hojnˇe ˇ vyuˇz´ıvat. Cerpat pˇritom budeme pˇredevˇs´ım z [2], m´enˇe pak z [22] a [20]. Naˇs´ı snahou bude vybrat optim´aln´ı funkci f˜(θ|D) ∈ F, kde F je prostor pravdˇepodobnostn´ıch funkc´ı, takovou, aby byla z konjugovan´eho syst´emu funkc´ı, tzn. f˘(θ|D) ∈ F c , a pˇritom byla co ˘ nejbl´ızˇ e aposteriorn´ı funkci f (θ|D). Chceme tedy minimalizovat m´ıru ∆ f (θ|D)kf(θ|D) , tzn. ˘ ˜ (3.14) f (θ|D) = arg min ∆ f (θ|D)kf(θ|D) . f˘∈Fc
Jako ide´aln´ı se jev´ divergence ı Kullback–Leibler (KLD) [28] a [29], tedy ˘ ˘ ∆ f (θ|D)kf(θ|D) = KLD f (θ|D)kf(θ|D) , kter´a je definov´ana jako:
Z f (θ|D) f (θ|D) ˘ KLD f (θ|D)kf(θ|D) = f (θ|D) log dθ = Ef (θ|D) log ˘ f˘(θ|D) f(θ|D) Θ
!
KL divergenci pro naˇse v´ypoˇcty pomoc´ı Variaˇcn´ıho Bayese budeme definovat jako: ˘ KLDV B = KLD f(θ|D)kf (θ|D) .
.
(3.15)
(3.16)
’ Vˇeta 2 (Variaˇcn´ı Bayesuv ˚ teor´ em) Necht f (θ|D) je aposteriorn´ı rozdˇelen´ı v´ıcerozmˇern´eho ′ ′ ′ ′ parametru θ = θ1 , θ2 , . . . , θq . Necht’ d´ale lze aproximaˇcn´ı rozdˇelen´ı f˘(θ|D) lze rozepsat jako produkt vz´ajemnˇe nez´avisl´ych rozdˇelen´ı pro θ1 , θ2 , . . . , θq , tzn.: ˘ f(θ|D) = f˘(θ1 , θ2 , . . . , θq |D) =
q Y i=1
˘ i |D) . f(θ
(3.17)
˘ ˜ (θ|D) , se nab´yv´a pro Pak minimum KLDV B , tj. f(θ|D) = arg min ∆ f(θ|D)kf f˘∈Fc
f˜(θi |D) ∝ exp Ef˜(θ/i |D) [ln(f (θ, D))] ,
Q kde θ/i znaˇc´ı komplement θi do θ a f˜(θ/i |D) = qj=1,j6=i f˜(θj |D) . Nazvˇeme f˜(θ|D) VB–aproximac´ı a f˜(θi |D) VB–margin´alou.
19
(3.18)
Algoritmus: (Iteraˇcn´ı VB (IVB) algoritmus) Tento algoritmus uk´azˇ eme pro q = 2, tedy θ = [θ1′ , θ2′ ]′ , z kontextu bude ale zˇrejm´e, zˇ e lze rozˇs´ıˇrit pro libovoln´e q. Cyklick´a iterace n´asleduj´ıc´ıch krok˚u monot´onnˇe sniˇzuje KLDV B : • V n–t´e iteraci spoˇc´ıt´ame VB–margin´alu f˜(θ2 |D) podle (3.18): Z [n−1] [n] f˜ (θ2 |D) ∝ exp f˜ (θ1 |D) ln f (θ1 , θ2 , D) dθ1
(3.19)
Θ1
• S vyuˇzit´ım pˇredchoz´ıho v´ysledku spoˇc´ıt´ame VB-margin´alu f˜(θ1 |D) dle (3.18): Z [n] [n] f˜ (θ1 |D) ∝ exp f˜ (θ2 |D) ln f (θ1 , θ2 , D) dθ2 (3.20) Θ2
Inicializaˇcn´ı hodnoty mohou b´yt voleny pomˇernˇe volnˇe, konvergence algoritmu byla dok´az´ana v [30]. VB metoda Nyn´ı zformujeme ucelen´y postup pro aplikaci VB teor´emu. Ten se bude skl´adat z osmi krok˚u, kter´e se pokus´ıme podrobnˇe popsat. Krok 1: (V´ybˇer modelu) Nejprve mus´ıme zkonstruovat Bayesovsk´y model, tzn. sdruˇzen´e rozdˇelen´ı dat a parametr˚u. Tento krok zahrnuje apriorn´ı pˇredpoklady, jako je volba f (D|θ) a apriorn´ıho rozdˇelen´ı f (θ). Pˇredpokl´adejme nyn´ı, zˇ e margin´aln´ı rozdˇelen´ı jednotliv´ych parametr˚u nelze analyticky nal´ezt. ˚ Mus´ıme zjistit, zda lze jednotliv´e parametry od sebe oddˇelit, Krok 2: (Oddˇelen´ı parametru) tedy zda je moˇzno pˇrirozen´y logaritmus sdruˇzen´eho rozdˇelen´ı rozepsat jako rozdˇelen´ı se separovan´ymi parametry. Pro q = 2 napˇr.: ln f (θ1 , θ2 , D) = f1 (θ1 , D)f2 (θ2 , D) ,
(3.21)
kde f1 (θ1 , D) a f2 (θ2 , D) jsou kompatibiln´ı vektory. Lze-li tuto operaci prov´est, pak ˇr´ık´ame, zˇ e sdruˇzen´e rozdˇelen´ı f (θ, D) patˇr´ı do rodiny funkc´ı separovateln´ych v parametrech. Zajistit tuto podm´ınku lze pr´avˇe vhodnou volbou apriorna, coˇz jsme jiˇz diskutovali v kapitole 3.1.1. Krok 3: (Zaps´an´ı VB-margin´al) Vˇse jsme pˇripravili na aplikaci VB teor´emu, konkr´etnˇe vzorce (3.18). T´ım z´ısk´ame VB-margin´aly pro θ1 , . . . , θq , tedy f˜(θ1 ), . . . , f˜(θq ). Krok 4: (Identifikace standartn´ıch forem) V tomto kroku je naˇs´ım u´ kolem identifikovat tvar VB-margin´al jako formu nˇekter´eho standartn´ıho parametrick´eho rozdˇelen´ı, napˇr: f˜(θi |D) ≡ f (θi |p1 , . . . , pkj )
20
∀i ∈ {1, . . . , p} ,
(3.22)
kde p1 , . . . , pkj jsou tvarovac´ı parametry zm´ınˇen´eho standartn´ıho rozdˇelen´ı. Tento krok m˚uzˇ e b´yt pomˇernˇe komplikovan´y a pˇri nevhodn´e volbˇe apriorna i neˇreˇsiteln´y. Pak nezb´yv´a neˇz zv´azˇ it, zda jsme zvolili vhodn´e pˇredpoklady nebo zda cel´y probl´em neˇreˇsit jinou metodou. ˚ Aby byla VB-aproximace kompletn´ı, zb´yv´a n´am jeˇstˇe Krok 5: (Formulace VB-momentu) formulovat takzvan´e VB-momenty, coˇz jsou momenty parametrick´ych rozdˇelen´ı z pˇredchoz´ıho kroku. Typicky jsou to funkce jejich tvarovac´ıch parametr˚u a pokud jsou parametrick´a rozdˇelen´ı klasick´a (zn´am´a), lze je vyhledat v tabulk´ach. V kroku 4 a 5 n´am tak vznikla soustava rovnic. Krok 6: (Redukce rovnic) Tento krok nen´ı v t´eto pr´aci vyuˇz´ıv´an a uv´ad´ıme ho s ohledem na kompletnost vzhledem k [2]. Jeho princip je takov´y, zˇ e se snaˇz´ıme o maxim´aln´ı zjednoduˇsen´ı soustavy rovnic vzhledem k poˇc´ıtaˇcov´e n´aroˇcnosti. Krok 7: (IVB algoritmus) Pokud se soustava rovnic z kroku 4 a 5 nepodaˇr´ı vyˇreˇsit analyticky (coˇz by vedlo k z´asadn´ımu urychlen´ı v´ypoˇctu), pak m´ame k dispozici IVB algoritmus, kter´y podle [30] konverguje ke skuteˇcn´ym VB-margin´al´am. Krok 8: (Vyps´an´ı VB-margin´al) IVB algoritmus n´am po ukonˇcen´ı, tzn. vhodn´em zastaven´ı, dod´a tvarovac´ı parametry hledan´ych aposteriorn´ıch rozdˇelen´ı. Tohoto postupu se budeme drˇzet i v n´asleduj´ıc´ıch kapitol´ach.
21
4 Faktorova´ analyza ´ obrazovych ´ sekvenc´ı Faktorov´a anal´yza je d˚uleˇzit´a souˇca´ st zpracov´an´ı dat v l´ekaˇrsk´e diagnostice. Jej´ım hlavn´ım pˇr´ınosem je anal´yza fyziologick´ych funkc´ı. Protoˇze je k dispozici cel´a sekvence sn´ımk˚u, m˚uzˇ eme zkoumat chov´an´ı org´an˚u v pr˚ubˇehu cˇ asu a t´ım je l´epe diagnostikovat. K tomu se obvykle pouˇz´ıv´a kontrastn´ı l´atka vyzaˇruj´ıc´ı cˇ a´ stice, kter´e sn´ım´a scintigrafick´a kamera. Z´akladn´ı myˇslenka a pˇredpoklady jsou uvedeny v kapitole 2.2.2. V t´eto kapitole sestroj´ıme a vyˇreˇs´ıme standardn´ı matematick´y model pro faktorovou anal´yzu scintigrafick´e obrazov´e sekvence, v´ysledky pak otestujeme na re´aln´e studii. Hlavn´ı inspirac´ı pˇri ˇreˇsen´ı tohoto probl´emu n´am bude pˇredevˇs´ım [2].
4.1 Konstrukce standardn´ıho modelu pro sekvenci sn´ımku˚ Zkonstruujme nyn´ı standardn´ı model dle [2]. Naˇs´ım u´ kolem bude analyzovat sekvenci n sn´ımk˚u poˇr´ızen´ych v cˇ ase t = 1, ..., n. Kaˇzd´y √ z poˇr´ızen´ych obr´azk˚u se skl´ad´a z p pixel˚u (strana cˇ tvercov´eho obr´azku je p), pˇreskl´ad´an´ım tedy vznikne p-dimenzion´aln´ı vektor kaˇzd´eho sn´ımku. Protoˇze m´ame celkem n sn´ımk˚u, bude m´ıt matice pozorov´an´ı tvar D ∈ Rp×n . D´ale pˇredpokl´adejme, zˇ e m´ame r faktor˚u (org´an˚u), kter´e chceme detekovat. Kaˇzd´y sn´ımek v sekvenci je pak tvoˇren line´arn´ı kombinac´ı jednotliv´ych faktor˚u. D´a se pˇredpokl´adat, zˇ e poˇcet faktor˚u bude menˇs´ı neˇz poˇcet sn´ımk˚u a z´aroveˇn poˇcet sn´ımk˚u bude o hodnˇe menˇs´ı neˇz poˇcet pixel˚u kaˇzd´eho sn´ımku, tedy r < n ≪ p. Fyziologick´y model m˚uzˇ eme zapsat ve tvaru: D = AX ′ + E ,
(4.1)
kde A ∈ Rp×r , X ∈ Rn×r a E ∈ Rp×n , sch´ema matic je vidˇet na obr. 4.1. Rozeberme si nyn´ı v´yznam jednotliv´ych matic. Kaˇzd´y pixel matice D je poˇcet dopad˚u cˇ a´ stic na danou cˇ a´ st scintigrafick´e kamery, coˇz znamen´a, zˇ e matice D je nez´aporn´a. Sloupce matice A interpretujeme jako jednotliv´e org´any (pˇr´ıpadnˇe faktory), proto je i A nez´aporn´a. Je zˇrejm´e, zˇ e tot´ezˇ plat´ı i pro matici X, jej´ızˇ sloupce reprezentuj´ı kˇrivku kaˇzd´eho faktoru. V modelu (4.1) pˇredstavuje matice D pozorovan´a data, A a X jsou nezn´am´e parametry a matice E je matice sˇumu, kde ei,j jsou i.i.d.1 s nezn´am´ym parametrem ω −1 (variance sˇumu) s 1
independent identically distributed . . . ei,j jsou nez´avisl´e a stejnˇe rozdˇelen´e
22
Obr´azek 4.1: Tvar a symbolika matic (pˇrevzato z [2])
maticov´ym rozdˇelen´ım: f (E|ω) = tNE 0p,n , ω −1Ip ⊗ In Modely 4.1 a 4.2 tak mohou b´yt zaps´any jako:
f (D|A, X, ω) = tND AX ′ , ω −1 Ip ⊗ In .
(4.2)
(4.3)
D´ale zvol´ıme Bayesovk´y model doplnˇen´ım modelu (4.3) o apriorn´ı rozdˇelen´ı matic A a X a skal´arn´ıho parametru ω: f (A|Υ) = tNA 0p,r , Ip ⊗ Υ−1 , (4.4)
kde
Υ = diag(v), v = [v1 , . . . , vr ]′ r Y f (v) = Gvi (αi,0 , βi,0 )
(4.5)
i=1
f (X) = tNX (0n,r , In ⊗ Ir ) f (ω) = Gω (ϑ0 , ρ0 ) .
(4.6) (4.7)
α0 a β0 ∈ Rr ; Υ ∈ Rr×r je diagon´aln´ı matice s hyperparametry vi , ϑ0 a ρ0 jsou n´ami zvolen´e apriorn´ı parametry. T´ım jsme detailnˇe popsali krok 1, tedy vytvoˇren´ı modelu, VB metody.
23
´ 4.1.1 Poznamka k maticove´ dekompozici Zamysleme se jeˇstˇe nad dekompozic´ı datov´e matice D, zjednoduˇsenˇe bez sˇumu. Pro n´azornost si nejdˇr´ıve vezmˇeme dekompozici skal´arn´ı. Mˇejme cˇ´ıslo d a hledejme jeho rozklad d = a · x. Moˇznost´ı, jak tuto rovnost splnit, je pro pevn´e d a promˇenn´e a a x nekoneˇcnˇe mnoho. Obdobn´y probl´em ˇreˇs´ıme i pro matici D, kter´a m˚uzˇ e m´ıt tvar D = AX ′ ,
(4.8)
D = A(T T −1 )X ′ ,
(4.9)
ale i kde matice T je kompatibiln´ı k A a X ′ . Vloˇzen´ım jednotkov´e matice I = T T −1 totiˇz opˇet splˇnujeme rovnici, ovˇsem dost´av´ame tak jin´e tvary org´an˚u i kˇrivek, nebot’ D = A(T T −1)X ′ = (AT )(T −1X ′ ) ,
(4.10)
e = AT A e ′ = T −1 X ′ X
(4.11)
takˇze nov´e matice
splˇnuj´ı rovnost
eX e′ . D=A
(4.12)
(4.13)
Matici T naz´yv´ame rotaˇcn´ı matic´ı. Probl´em toho typu se objevuje cˇ asto a je ˇreˇsen napˇr. v [31]. Z´akladn´ım u´ kolem je vybrat z nekoneˇcnˇe mnoha jedno nejpravdˇepodobnˇejˇs´ı rozdˇelen´ı, kter´e obecnˇe splˇnuje vˇsechny poˇzadovan´e podm´ınky. Pro jednorozmˇern´y pˇr´ıpad viz [2].
4.2 VB metoda pro standardn´ı model (4.1) Krok 2 Spoˇcteme sdruˇzen´e rozdˇelen´ı (4.3) - (4.7), tzn. f (D, A, X, ω, Υ|r), a najdeme jeho pˇrirozen´y logaritmus: pn 1 ln ω − ω tr ((D − AX ′ )(D − AX ′ )′ ) + 2 2 r X 1 p 1 ln vi − tr(ΥA′ A) − tr(XX ′ )+ + 2 i=1 2 2
ln f (D, A, X, ω, Υ|r) =
+
r X i=1
(α0 − 1) ln vi +
r X i=1
β0 vi + (ϑ0 − 1) ln ω − ρ0 ω + γ , (4.14)
kde γ pˇredstavuje souˇcet v´yraz˚u, kter´e nez´avis´ı na parametrech A, X, ω ani Υ. Tento v´yraz nyn´ı splˇnuje podm´ınky Bayesova variaˇcn´ıho teor´emu pro θ1 = A, θ2 = X, θ3 = ω, θ4 = Υ a patˇr´ı tedy do rodiny funkc´ı separovateln´ych v parametrech (splˇnuje ln f (θ1 , . . . , θq , D) = g1 (θ1 , D) . . . gq (θq , D)).
24
Krok 3 Podle Bayesova variaˇcn´ıho teor´emu (vzorec (3.18)) maj´ı parametry A, X, ω a v n´asleduj´ıc´ı VB-margin´aly: 1 1 1 ′ ′ ′ ′ ′ X)A )) − ω b D ) − tr(A(b b ) [ b tr(−2AX ωX b tr(AΥA f˜(A|D, r) ∝ exp − ω 2 2 2 1 1 1 ′ ′ ′ ′ ′ A)X )) − ω b D ) − tr(X(b d f˜(X|D, r) ∝ exp − ω b tr(−2AX ωA b tr(XX ) 2 2 2 " r # r X X p 1 ′ A)) d + αi,0 − 1 ln vi − βi,0 + (diag(A f˜(v|D, r) ∝ exp i 2 2 i=1 i=1 " pn ˜ f (ω|D, r) ∝ exp + ϑ0 − 1 ln ω+ 2 # 1 1 ′ ′ ′ ′ ′ AX ′X bX b D − DX bA b + tr A d [ − ω ρ0 + tr DD − A 2 2
Krok 4 Budeme pˇredpokl´adat, zˇ e v´ysˇe uveden´e VB-margin´aly maj´ı n´asleduj´ıc´ı tvary: f˜(A|D, r) = tNA (µA , Ip ⊗ ΣA ) ˜ f(X|D, r) = tNX (µX , Ip ⊗ ΣX ) ˜ f(v|D, r) =
r Y
Gvi (αi , βi )
(4.15) (4.16) (4.17)
i=1
f˜(ω|D, r) = Gω (ϑ, ρ)
(4.18)
a postupnˇe vyj´adˇr´ıme tvarovac´ı parametry µA , ΣA , µX , ΣX , α, β, ϑ a ρ. Pro ilustraci pˇredvedeme v´ypoˇcet pro z´ısk´an´ı µA a ΣA , pro ostatn´ı parametry je postup obdobn´y. 1 −1 −1 ′ f˜(A|D, r) = tNA (µA , Ip ⊗ ΣA ) ∝ exp − tr[Ip (A − µA )(ΣA ) (A − µA )] 2 Nyn´ı porovn´ame kvadratick´e a line´arn´ı cˇ leny: ′ X)A′ + AΥA b ′ = A(Σ−1 )′ A′ [ A(b ωX A ′ b ′ ′ −2A XDb ω = −µA (Σ−1 A )A
25
a s vˇedom´ım toho, zˇ e ΣA je symetrick´a matice, snadno vyj´adˇr´ıme µA a ΣA . Celkovˇe tedy dost´av´ame soustavu rovnic: −1 ′X + Υ b [ ΣA = ω bX (4.19) −1 ′X + Υ b ω b [ µA = ω bDX bX (4.20) −1 ′A + I d ΣX = ω bA (4.21) r −1 ′A + I b ω d µX = ω bD′ A bA (4.22) r p α = α0 + 1r,1 (4.23) 2 1 ′A d (4.24) β = β0 + diag A 2 np ϑ = ϑ0 + (4.25) 2 1 ′ AX ′X bX b ′D′ − DX bA b′ + 1 tr A d [ (4.26) ρ = ρ0 + tr DD ′ − A 2 2
Krok 5
Nyn´ı n´am jeˇstˇe zb´yv´a dopoˇc´ıtat VB-momenty jednotliv´ych distribuc´ı. K tomu pouˇzijeme vzorce pro maticov´e norm´aln´ı rozdˇelen´ı a pro gamma rozdˇelen´ı z kapitol A.2 a A.3. Pomoc´ı tˇechto vzorc˚u nyn´ı vyj´adˇr´ıme z (4.15) – (4.18) momenty jako funkce tvarovac´ıch parametr˚u: b =µA A ′ A =pΣ + µ′ µ d A A A A b =µX X
Krok 7
′ X =nΣ + µ′ µ [ X X X X αi vbi = , i = 1, . . . , r βi ϑ ω b= ρ
(4.27)
(4.28)
Dost´av´ame soustavu rovnic (4.19) – (4.28), kter´a je ovˇsem analyticky prakticky neˇreˇsiteln´a, obzvl´asˇt’ pouˇzijeme-li pozdˇeji m´ısto norm´aln´ıho rozdˇelen´ı oˇrezan´e norm´aln´ı rozdˇelen´ı, abychom zajistili pozitivitu matic. Proto jsme nuceni pˇristoupit k numerick´emu ˇreˇsen´ı pomoc´ı IVB algoritmu: • nastav´ıme vhodn´e inicializaˇcn´ı hodnoty a v k-t´em kroku iterace spoˇcteme s vyuˇzit´ım hodnot z (k − 1)-kroku postupnˇe:
26
′A bA d • µA , ΣA , A,
′X b X [ • µX , ΣX , X,
• α, β, ϑ, ρ b. • ω b, Υ
T´ım z´ısk´av´ame odhad pro matice A a X.
Krok 8 Tento krok je autorem proveden a diskutov´an jiˇz v [32] na testovac´ıch datech. V n´asleduj´ıc´ı kapitole provedeme cel´y v´ypoˇcet na re´aln´ych datech a zamysl´ıme se nad v´ysledky.
4.3 Vysledky ´ standardn´ıho modelu Uvedenou metodu nyn´ı vyzkouˇs´ıme na re´aln´e obrazov´e sekvenci. Cel´y v´ypoˇcet je implementov´an v Matlabu. Nejprve se vˇenujme popisu samotn´e sekvence, jej´ım z´akladn´ım parametr˚um. Studie, na kter´e budeme pˇredv´adˇet v´ysledky standardn´ıho modelu, je z´ısk´ana ze zdrav´eho pacienta. M˚uzˇ eme na n´ı tedy pozorovat pˇribliˇznˇe stejn´e intenzity jak parenchym˚u, tak p´anviˇcek. Hlavn´ımi faktory, kter´e d´ale obsahuje, jsou srdce (a krevn´ı pozad´ı) a moˇcov´y mˇech´yˇr. Samotnou sekvenci m´ame zobrazenou na obr´azku 4.2.
10
10
10
10
10
10
20
20
20
20
20
20
30
30
30
30
30
30
40
40
40
40
40
40
50
50
50
50
50
60
60
60
60
60
20
40
60
20
40
60
20
40
60
20
40
60
50 60 20
40
60
10
10
10
10
10
10
20
20
20
20
20
20
30
30
30
30
30
30
40
40
40
40
40
40
50
50
50
50
50
50
60
60
60
60
60
20
40
60
20
40
60
20
40
60
20
40
60
40
60
10
10
10
10
10
20
20
20
20
20
20
30
30
30
30
30
30
40
40
40
40
40
40
50
50
50
50
50
60 20
40
60
60 20
40
60
60 20
40
60
40
60
60
20
40
60
20
40
60
50
60 20
40
60 20
10
60
20
60 20
40
60
Obr´azek 4.2: scintigrafick´a studie ledvin (kaˇzd´y sedm´y sn´ımek)
V naˇsem pˇr´ıpadˇe m´a jeden sn´ımek rozmˇery 64×64, z toho tedy plyne, zˇ e p = 4096. Sn´ımk˚u je celkem 120 (tedy n = 120) a jsou sn´ım´any vˇzdy po 10s.
27
Protoˇze se jedn´a o zdrav´eho pacienta se standardnˇe uloˇzen´ymi ledvinami, lze ˇr´ıci, zˇ e je tato studie pomˇernˇe lehk´a, neobjevuj´ı se tu zˇ a´ dn´e patologick´e jevy.
´ 4.3.1 Faktorove´ obrazky a jejich kˇrivky V prvn´ı cˇ a´ sti v´ysledk˚u n´am p˚ujde pouze o zobrazen´ı jednotliv´ych faktor˚u a kˇrivek pr˚utoku radiofarmaka v nich. Dalˇs´ı anal´yze se pak budeme vˇenovat v n´asleduj´ıc´ıch kapitol´ach. V´ysledek rozkladu datov´e matice na faktory a jejich kˇrivky m˚uzˇ eme vidˇet na obr´azku 4.3.
15 10 20
10
30 40
5
50 60 20
40
0
60
0
50
100
150
0
50
100
150
0
50
100
150
0
50
100
150
10 10
8
20 6
30 40
4
50
2
60 20
40
0
60
8 10 6
20 30
4
40 2
50 60 20
40
0
60
20 10 15
20 30
10
40 5
50 60 20
40
0
60
Obr´azek 4.3: v´ysledek anal´yzy studie ledvin pomoc´ı standardn´ıho modelu Prvn´ı nalezen´y faktor jsou oba parenchymy. Protoˇze maj´ı stejnou dynamiku, m´ame je urˇcen´e jako jeden faktor. Jejich odhad dopadl na prvn´ı pohled spr´avnˇe, pˇribliˇzsˇ´ım pohledu si vˇsak m˚uzˇ eme vˇsimnout, zˇ e se k tomuto obr´azku pˇriˇcetlo i nˇejak´e pozad´ı, pravdˇepodobnˇe tk´anˇ ov´e. To samozˇrejmˇe nutnˇe ovlivn´ı i v´yslednou kˇrivku.
28
Druh´ym faktorem je p´anviˇcka. Tak je nalezena spr´avnˇe a pˇri bliˇzsˇ´ım pohledu si m˚uzˇ eme vˇsimnout dokonce i odtoku do moˇcov´eho mˇech´yˇre. Jej´ı kˇrivka je na poˇca´ tku sekvence nulov´a, stoupat zaˇc´ın´a po pˇribliˇznˇe 110s, tedy dle biologick´eho pˇredpokladu. Jako tˇret´ı faktor se n´am povedlo urˇcit moˇcov´y mˇech´yˇr. Vid´ıme, zˇ e se n´am k nˇemu pˇriˇcetly i cˇ a´ sti ledvin, pˇredevˇs´ım p´anviˇcky. Protoˇze standardn´ı model nepˇredpokl´ad´a kompaktnost faktor˚u, nen´ı to chyba v´ypoˇctu, nˇekter´e pixely maj´ı zkr´atka dynamiku podobnou moˇcov´emu mˇech´yˇri, proto se k nˇemu pˇriˇcetly. Obvykle se pˇred samotnou anal´yzou m˚uzˇ e moˇcov´y mˇech´yˇr oˇr´ıznout, nebot’ jeho spr´avn´a detekce je pomˇernˇe n´aroˇcn´a a pro samotnou anal´yzu ledvin nem´a smysl. ˇ Ctvrt´ y faktor n´am tvoˇr´ı srdce a krevn´ı pozad´ı. V lev´e cˇ a´ sti obr´azku si m˚uzˇ eme tak´e vˇsimnout v´yraznˇejˇs´ı struktury, kter´a vznikla aplikac´ı radiofarmaka do pacientovi prav´e ruky. Zde m˚uzˇ eme b´yt nejm´enˇe spokojeni s odhadnutou kˇrivkou, kter´a po cca tˇrech minut´ach klesla na nulu a u n´ı se jiˇz drˇz´ı. To neodpov´ıd´a biologick´emu pˇredpokladu.
ˇ 4.3.2 Tranzitn´ı cas Samotn´e nalezen´ı jednotliv´ych faktor˚u a jejich pr˚ubˇeh˚u je jen prvn´ım krokem anal´yzy scintigrafick´e obrazov´e sekvence ledvin. Je tˇreba urˇcit parametr, kter´y pokud moˇzno co nejl´epe pop´ısˇe, jak dan´y faktor funguje. V kapitole 2.3.3 jsme dva takov´e parametry uvedli, pod´ıvejme se nejprve na tranzitn´ı cˇ as (TT) parenchymu. Jeho konstrukce a zp˚usob urˇcen´ı jsou naznaˇceny na obr´azku 2.12. 0.3
0.25
0.2
0.15
0.1
0.05
0
−0.05
0
20
40
60
80
100
120
Obr´azek 4.4: impulzn´ı retenˇcn´ı funkce parenchymu urˇcen´a standardn´ım modelem
Zkusme urˇcit napˇr´ıklad TT parenchymu. K tomu je zapotˇreb´ı z´ıskat impulzn´ı retenˇcn´ı funkci, kter´a by mˇela teoreticky vyj´ıt z dekonvoluce kˇrivky parenchymu a kˇrivky krve.
29
Prvn´ı probl´em, na kter´y nar´azˇ´ıme, je automaticky urˇcit, kter´y faktor reprezentuje kterou strukturu lidsk´eho tˇela. To n´asˇ model explicitnˇe nedˇel´a, nicm´enˇe m˚uzˇ eme s pomˇernˇe velk´ym u´ spˇechem nastavit urˇcit´e poˇca´ teˇcn´ı typick´e hodnoty kˇrivek, pouˇz´ıt je jako startovac´ı data a pˇredpokl´adat, zˇ e faktor s danou kˇrivkou se odhadne jako struktura s t´ımto pr˚ubˇehem. To samozˇrejmˇe ˇreˇs´ı pouze typick´e pˇr´ıpady, nikoliv patologick´e chov´an´ı nˇejak´eho faktoru. Pro n´asˇ experiment n´am staˇc´ı vizu´alnˇe urˇcit parenchym a krev, coˇz nen´ı probl´em, a prov´est dekonvoluci. To m˚uzˇ eme udˇelat napˇr. pomoc´ı Fourierovy transformace (FT). Provedeme FT kˇrivky parenchymu a FT kˇrivky krve, vydˇel´ıme tyto dva v´ysledky a provedeme inverzn´ı FT. V´ysledek tohoto postupu v naˇsem experimentu vid´ıme na obr´azku 4.4. Teoretick´y tvar t´eto funkce jsme jiˇz diskutovali a m˚uzˇ eme se na nˇej pod´ıvat na obr´azku 2.6. Vid´ıme, zˇ e urˇcen´ı pˇresn´eho TT z t´eto kˇrivky nen´ı jednoznaˇcn´e, protoˇze vizu´alnˇe urˇcit, kde konˇc´ı konstantn´ı cˇ a´ st a z´acˇ´ın´a kles´an´ı, nen´ı prakticky moˇzn´e. Tento postup je velmi citliv´y na pˇresn´e urˇcen´ı obou kˇrivek, pˇredevˇs´ım kˇrivky krve. V praxi tedy TT pomoc´ı tohoto modelu automaticky urˇc´ıme velmi tˇezˇ ko a nepˇresnˇe.
10
20
30
40
50
60
10
20
30
40
50
60
Obr´azek 4.5: sn´ımek parenchym˚u, ze kter´eho chceme urˇcit relativn´ı ren´aln´ı clearance
´ ı clearance 4.3.3 Relativn´ı renaln´ Urˇcen´ı relativn´ı ren´aln´ı clearance (viz kapitola 2.3.3) je dalˇs´ım u´ kolem anal´yzy. Urˇcit ji m˚uzˇ eme z v´ysledn´eho obr´azku (obr´azk˚u) parenchym˚u, opˇet za pˇredpokladu, zˇ e je dok´azˇ eme rozpoznat. K tomu staˇc´ı tentokr´at pouˇz´ıt pouze prvn´ı cˇ a´ st sekvence, napˇr. do chv´ıle, kdy nab´yv´a kˇrivka parenchymu sv´eho maxima, zjednoduˇs´ıme si tak u´ lohu o rozpozn´an´ı p´anviˇcek. V naˇsem
30
pˇr´ıpadˇe to znamen´a pouˇz´ıt pouze prvn´ıch 17 sn´ımk˚u. V´ysledn´y obr´azek parenchym˚u pak m˚uzˇ eme vidˇet na obr´azku 4.5. Nast´av´a tu ovˇsem probl´em s nedokonal´ym odstranˇen´ım sˇ umu. Pokud bychom pouze seˇcetli hodnoty pixel˚u z obou p˚ulek obr´azku, dostaneme zkreslen´y v´ysledek, tento jednoduch´y postup nelze pouˇz´ıt. Nezb´yv´a n´am tedy nic jin´eho, neˇz urˇcit oblast ledvin a seˇc´ıst pixely jen v t´eto oblasti. To se vˇsak jiˇz v klinick´e praxi dˇeje a pokud n´asˇ postup obsahuje krok ruˇcn´ıho ohraniˇcov´an´ı, nen´ı to zˇ a´ dn´y v´yznamn´y pˇr´ınos.
4.4 Zhodnocen´ı vysledk ´ u˚ Uk´azali jsme konstrukci standardn´ıho modelu pro funkcion´aln´ı anal´yzu scintigrafick´ych studi´ı ledvin a n´aslednˇe jeho ˇreˇsen´ı. V´ypoˇcet jsme pˇredvedli na studii zdrav´eho pacienta a pokusili se urˇcit dva diagnostick´e koeficienty. Samotn´y rozklad na obr´azky a kˇrivky dopadl v r´amci oˇcek´av´an´ı. Probl´emy se uk´azaly pˇri urˇcov´an´ı obou koeficient˚u, tranzitn´ıho cˇ asu parenchymu i relativn´ı ren´aln´ı clearance. Pro urˇcen´ı spr´avn´eho TT se n´am nepodaˇrilo urˇcit dostateˇcnˇe pˇresn´y tvar impulzn´ı retenˇcn´ı funkce. V´ypoˇcet spr´avn´e relativn´ı ren´aln´ı clearance by pˇredpokl´adal u´ pln´e ostranˇen´ı sˇumu, coˇz standardn´ı model nedok´azal. To by nebyl tak velk´y probl´em, pokud bychom umˇeli automaticky urˇcit polohy obou parenchym˚u. Tuto informaci vˇsak standardn´ı model neodhaduje, nen´ı to tedy zlepˇsen´ı oproti st´avaj´ıc´ı klinick´e praxi, protoˇze bychom parenchymy museli opˇet oznaˇcit ruˇcnˇe. Ne´uspˇech pˇri urˇcov´an´ı diagnostick´ych parametr˚u je cˇ a´ steˇcnˇe zp˚usoben konstrukc´ı modelu, kter´y plnˇe nerespektuje z´akladn´ı biologick´e pˇredpoklady. To je motivac´ı pro dalˇs´ı kapitolu, ve kter´e se budeme vˇenovat nov´e konstrukci modelu, kter´y se bude snaˇzit tyto nedostatky odstranit.
31
´ ı model pro funcionaln´ ´ ı 5 Integraln´ analyzu ´ obrazovych ´ sekvenc´ı V kapitole 4 jsme uk´azali, jak´ym zp˚usobem je konstruov´an standardn´ı model pro anal´yzu scintigrafick´ych obrazov´ych sekvenc´ı. Uk´azali jsme tak´e jeho ˇreˇsen´ı a n´aslednˇe v´ysledky, kter´ych pomoc´ı nˇej m˚uzˇ eme dos´ahnout. Nedostatkem tohoto modelu je, zˇ e jeho stavba neodpov´ıd´a tomu, o cˇ em pojedn´avala kapitola 2.2, nerespektuje tedy zn´am´y biologick´y model, rozkl´ad´a pouze datovou matici D na AX ′ s ohledem na sˇum. Kromˇe pozitivity neklade zˇ a´ dn´e dalˇs´ı poˇzadavky. To znamen´a, zˇ e v´ysledky splˇnuj´ı matematick´y rozklad, ovˇsem biologick´a spr´avnost kˇrivek nen´ı zaruˇcena. Pˇritom pr´avˇe na z´akladˇe nich by mˇela prob´ıhat diagnostick´a anal´yza. Nev´yhodou tak´e je, zˇ e se tento model nesnaˇz´ı doj´ıt k hotov´ym v´ysledk˚um (pˇr´ıpadnˇe tvaru, ze kter´eho lze v´ysledky pˇr´ımo vyˇc´ıst), ale je nutn´a dalˇs´ı anal´yza v´ysledn´ych kˇrivek. Z´ıskan´y v´ysledek je sice automatick´y (nen´ı nutn´e napˇr. oznaˇcovat pˇredem ROI), coˇz zmenˇsuje pravdˇepodobnost chyby lidsk´eho faktoru, ale n´asledn´a nezbytn´a anal´yza n´am chybu opˇet zan´asˇ´ı, nebot’ pˇri t´eto operaci doch´az´ı opˇet ke ztr´atˇe informace, kterou jiˇz model nem´a sˇanci nikde kompenzovat. Tyto probl´emy bychom r´adi ˇreˇsili v t´eto pr´aci. Pokus´ıme se tedy o synt´ezu informac´ı z kapitoly 2.3, cˇ´ımˇz bychom chtˇeli vybudovat integr´aln´ı model postaven´y na re´aln´ych biologick´ych pˇredpokladech. Ide´aln´ı by tak´e bylo zabudovat odhadov´an´ı d˚uleˇzit´ych parametr˚u pˇr´ımo do v´ypoˇctu, jejich hodnoty se tedy pokus´ıme dostat jako vedlejˇs´ı produkt cel´e funkcion´aln´ı anal´yzy.
5.1 Konstrukce modelu N´asˇ nov´y model bude zaloˇzen opˇet na rozkladu datov´e matice na dvˇe dalˇs´ı, kde sloupce prvn´ı matice budou reprezentovat obr´azky a sloupce druh´e matice budou tvoˇrit kˇrivky. Form´alnˇe tedy m˚uzˇ eme zapsat model opˇet ve tvaru D = AX ′ + E, kde D ∈ Rp×n , A ∈ Rp×r , X ∈ Rn×r a E ∈ Rp×n . Zachov´av´ame zde v´yznam rozmˇer˚u, p je poˇcet pixel˚u jednoho obr´azku, n poˇcet sn´ımk˚u obrazov´e sekvence a r je pˇredpokl´adan´y poˇcet faktor˚u. Z´apis matice kˇrivek jako matice X budeme i nad´ale vyuˇz´ıvat, ovˇsem nebudeme odhadovat pˇr´ımo samotnou matici X, ale parametry, pomoc´ı kter´ych budeme X jednoznaˇcnˇe reprezentovat. Tyto parametry budeme budovat na biologick´em podkladˇe, cˇ´ımˇz se pokus´ıme namodelovat poˇzadavky, kter´e na n´as klade kapitola 2.2.
32
´ 5.1.1 Zakladn´ ı model dat Pˇredevˇs´ım budeme kaˇzdou kˇrivku modelovat jako konvoluci kˇrivky krve a impulzn´ı retenˇcn´ı funkce kaˇzd´eho faktoru dle kapitoly 2.3.2. Kaˇzdou kˇrivku tak rozloˇz´ıme na tvar xt,f =
t X
bt−m+1 um,f ,
(5.1)
m=1
kde b ∈ Rn×1 je vektor tvoˇr´ıc´ı kˇrivku krve a uf je vektor obsahuj´ıc´ı impulzn´ı retenˇcn´ı funkci pro kaˇzd´y faktor. Celkovˇe tedy dostaneme matici U ∈ Rn×r , jej´ızˇ sloupce tvoˇr´ı uf jednotliv´ych faktor˚u. T´ım budeme pˇr´ımo odhadovat veliˇcinu, kterou chceme z´ıskat pro dalˇs´ı anal´yzu. Abychom s n´ı mohli d´ale pracovat a z´aroveˇn jsme zajistili jako pˇredpoklad jej´ı pozitivitu a kumulativnost, rozloˇz´ıme kaˇzd´y prvek um,f jako um,f =
n X
wφ,f ,
(5.2)
φ=m
cˇ´ımˇz dostaneme matici W ∈ Rn×r , jej´ızˇ sloupce tvoˇr´ı pˇr´ır˚ustky impulzn´ı retenˇcn´ı funkce uf . Prov´azanost prvk˚u matic U a W vid´ıme na obr´azku 5.1.
y u1,f w1,f u2,f
x Obr´azek 5.1: sch´ema prov´azanosti jednotliv´ych prvk˚u matic U a W
Obdobn´ym zp˚usobem budeme modelovat i kˇrivku krve, vektor b: bs =
n X
gς ,
ς=s
dost´av´ame tedy vektor g ∈ Rn×1 , jakoˇzto pˇr´ır˚ustkov´y vektor vektoru b.
33
(5.3)
Ztratili jsme t´ım sice moˇznost elegantn´ıho z´apisu datov´e matice D, ale zajistili jsme z´akladn´ı biologick´e pˇredpoklady. Je-li matice A tvoˇrena sloupcov´ymi vektory af ∈ Rp×1 (reprezentuj´ıc´ımi kaˇzd´y faktor), pak m˚uzˇ eme matici D zapsat jako !′ t r r X X X ′ af bt−m+1 um,f = af xf = D= f =1
=
r X f =1
af
f =1
t X
m=1
m=1
n X
ς=(t−m+1)
t=1..n
gς
n X
φ=m
′
wφ,f
t=1..n
.
(5.4)
ˇ ren´ı 5.1.2 Model chyb meˇ Pod´ıv´ame se na model matice D jako celku: vol´ıme distribuci f (D|A, X, ω) = tN(AX ′ , ω −1Ip ⊗ In ),
(5.5)
tedy pˇredpokl´ad´ame, zˇ e D m´a norm´aln´ı rozdˇelen´ı s diagon´aln´ı kovarianˇcn´ı matic´ı ovlivnˇenou parametry ω, kter´a reprezentuj´ı varianci sˇumu. Pˇresnˇeji bychom mˇeli zapsat apriorn´ı rozdˇelen´ı matice D jakoˇzto Poissonovo, ovˇsem jiˇz jsme uvedli, zˇ e to lze dobˇre aproximovat pro n > 50 ˇ jako norm´aln´ı. Sum aproximujeme parametrem ω, jemuˇz pˇriˇrad´ıme apriorno f (ω) = Gω (ϑ0 , ρ0 ),
(5.6)
s apriorn´ımi parametry ϑ0 , ρ0 . K pˇriˇrazen´ı stejn´e variance vˇsem pixel˚um jsme pˇrikroˇcili pˇredevˇs´ım z v´ypoˇcetn´ıch d˚uvod˚u, pˇresnˇejˇs´ı by bylo odhadovat kovarianˇcn´ı matici D ne jako (ω −1Ip ⊗ In ), ale napˇr´ıklad jako (ω −1 Ωp ⊗ Ωn ), kde Ωp a Ωn by byly diagon´aln´ı matice. Tˇem bychom opˇet pˇriˇradili apriorna a odhadovali je v pr˚ubˇehu v´ypoˇctu. T´ım bychom mˇeli l´epe splnit pˇredpoklady modelu, odhadov´an´ı t´eto kovariance vˇsak nen´ı pˇredmˇetem pr´ace. Podrobnˇe je tento probl´em rozebr´an v [2].
´ ı matice obrazk ´ u˚ 5.1.3 Modelovan´ Do modelu matice A bychom chtˇeli zav´est informaci o poloze nalezen´eho faktoru. To se n´am bude hodit napˇr. pro automatick´y v´ypoˇcet relativn´ı ren´aln´ı clearance. Mˇejme mnoˇzinu C, kter´a urˇcuje, zda patˇr´ı dan´y pixel matice k dan´emu faktoru, kaˇzd´y bod t´eto mnoˇziny tedy m˚uzˇ e nab´yvat hodnot 1 nebo 0. To urˇc´ıme z toho, zda je pravdˇepodobnˇejˇs´ı, aby pixel n´aleˇzel cˇ i nen´aleˇzel k dan´emu faktoru, tedy zda je jeho hodnota bl´ızˇ e k jist´e hodnotˇe cf , cˇ i k 0. Dostaneme tak vektor c ∈ Rr×1 . Poznamenejme, zˇ e mnoˇzina C tvoˇr´ı automatick´e ROI jednotliv´ych faktorov´ych obr´azk˚u. Celou tuto u´ vahu m˚uzˇ eme vyj´adˇrit jako: ( cf i ∈ C . (5.7) (zaf )i = 0 jinak
34
Pak m˚uzˇ eme ps´at apriorn´ı pˇredpoklad: (zaf )i = cf Ci,f .
(5.8)
Kaˇzd´y parametr cf bude m´ıt norm´aln´ı rozdˇelen´ı, dohromady m˚uzˇ eme ps´at f (c) = tN(0r,1 , ς0 Ir ),
(5.9)
kde apriorn´ı parametr ς0 je dostateˇcnˇe velk´e cˇ´ıslo, zajiˇst’uj´ıc´ı meze kovariance vektoru c. Apriorno jednoho sloupce af urˇc´ıme jako: f (af ) = tN(zaf , υf−1Ip )
(5.10)
a matice A, jej´ızˇ sloupce se skl´adaj´ı z vektor˚u a m´a apriorn´ı rozdˇelen´ı: f (A) = tN(ZA , Ip ⊗ Υ−1 A ).
(5.11)
Zde ZA = [za1 za2 . . . zar ] a ΥA ≡ Υ je diagon´aln´ı matice skl´adaj´ıc´ı se z hyperparametr˚u υf pˇr´ısluˇsn´ych k jednotliv´ym faktor˚um, kter´e slouˇz´ı k jejich lepˇs´ı separaci. Kaˇzd´emu hyperparametru υf pˇrisoud´ıme apriorn´ı rozdˇelen´ı gamma. To m˚uzˇ eme zapsat n´asleduj´ıc´ım zp˚usobem: Υ = diag(υ) υ = [υ1 , . . . , υr ]′ r Y f (υ) = Gυi (αi,0 , βi,0 ).
(5.12)
(5.13)
i=1
T´ım m´ame plnˇe pops´ano apriorn´ı rozdˇelen´ı kladen´e na matici A. Pozn´amka 3 Obdobn´y mechanizmus, jak´ym pouˇz´ıv´ame mnoˇzinu C jako ROI na ZA , by sˇel zav´est pˇr´ımo na matici A, v praxi to vˇsak nen´ı moˇzn´e kv˚uli v´ypoˇcetn´ı n´aroˇcnosti.
´ ı impulzn´ıch retencn´ ˇ ıch funkc´ı 5.1.4 Modelovan´ Pomoc´ı matice pˇr´ır˚ustk˚u, W , m˚uzˇ eme ovlivˇnovat a volit tvar impulzn´ıch retenˇcn´ı funkc´ı, uloˇzen´ych v matici U. Proto chceme pˇripravit takov´e apriorno, kter´e urˇc´ı pˇr´ır˚ustky jako konstantn´ı impulz o urˇcit´e d´elce. Z´aroveˇn tak zabudujeme do modelu pˇr´ım´e odhadov´an´ı tranzitn´ıch cˇ as˚u jednotliv´ych faktor˚u. Necht’ ( hf i ∈ Hf , (5.14) wi,f = 0 jinak kde Hf je spojit´a indexov´a mnoˇzina, Hf = {isf , isf +1, . . . , isf +lf }, f ∈ {1, . . . , r}, celkovˇe tedy dost´av´ame H = {H1 , . . . , Hr }. Pro popis t´eto mnoˇziny bude v´yhodn´e zav´est vektory
35
s = [s1 , . . . , sr ]′ , l = [l1 , . . . , lr ]′ a h = [h1 , . . . , hr ]′ , jejichˇz v´yznam je zˇrejm´y z konstrukce mnoˇziny H a z obr´azku 5.2. Apriorna pˇriˇrad´ıme tˇemto vektor˚um n´asledovnˇe: f (h) = tN(0r,1 , τ0 Ir ) f (si) = U(0, n) ∀i ∈ {1, . . . , r} f (li |si) = U(0, n − si ) ∀i ∈ {1, . . . , r},
(5.15) (5.16) (5.17)
kde apriorn´ı parametr τ0 je dostateˇcnˇe velk´e cˇ´ıslo, zajiˇst’uj´ıc´ı meze kovariance vektoru h. Vektor˚um s a l pˇriˇrazujeme rovnomˇern´e rozdˇelen´ı (viz A.5), nebot’ nem´ame zˇ a´ dnou preferenci na jejich hodnoty. y
uf
hf mwf
sf
ˇcas lf
Obr´azek 5.2: Zn´azornˇen´ı konstrukce apriorna matice W , pˇredstavuj´ıc´ı pˇr´ır˚ustky matice U
Apriorno wf m´a pak tvar: f (wf ) = tN(mwf (H), ξf−1In ),
(5.18)
kde mwf (H) je vznikl´e apriorno kladen´e na sloupce matice W a samotn´a matice W , jej´ızˇ sloupce se skl´adaj´ı z vektor˚u w m´a apriorn´ı rozdˇelen´ı: f (W ) = tN(MW , In ⊗ Ξ−1 W ).
36
(5.19)
Zde MW = [mw1 mw2 . . . mwr ] a ΞW je diagon´aln´ı matice skl´adaj´ıc´ı se z hyperparametr˚u ξf pˇr´ısluˇsn´ych k jednotliv´ym faktor˚um. Ty maj´ı opˇet apriorn´ı rozdˇelen´ı gamma, obdobnˇe jako u matice A. Celkovˇe tedy: Ξ = diag(ξ) ξ = [ξ1 , . . . , ξr ]′ r Y f (ξ) = Gξi (κi,0 , νi,0 ).
(5.20) (5.21) (5.22)
i=1
T´ım jsme plnˇe popsali mechanismus odhadov´an´ı impulzn´ıch retenˇcn´ıch funkc´ı.
5.1.5 Model kˇrivky krve Obdobn´ym zp˚usobem jako v pˇredchoz´ı kapitole budeme modelovat i kˇrivku krve. U n´ı pˇredpokl´ad´ame, zˇ e m´a klesaj´ıc´ı charakter (to nemus´ı b´yt re´alnˇe vˇzdy zaruˇceno) a jej´ı hodnota je nez´aporn´a. Prvky vektoru b maj´ı tvar: n X gς , (5.23) bs = ς=s
tvarujeme je tedy jako pˇr´ır˚ustky, pomoc´ı nichˇz m´ame mnohem vˇetˇs´ı vliv na tvar vektoru b, neˇz kdybychom ho modelovali samostatnˇe. Zavedli jsme pˇr´ır˚ustkov´y vektor g, kde f (g|ψ) = tN(0n , ψ −1 In ),
(5.24)
hyperparametr ψ opˇet urˇcuje velikost kovarianˇcn´ı matice vektoru g a m´a rozdˇelen´ı gamma: f (ψ) = Gψ (ζ0 , η0 )
(5.25)
s apriorn´ımi parametry ζ0 , η0 .
´ ˇ e´ 5.1.6 Odhadovane´ a zavisl e´ promenn Protoˇze v cel´em modelu vyuˇz´ıv´ame pomˇernˇe velk´y poˇcet promˇenn´ych, P pod´ıvejme se nyn´ı na jejich tvary a prov´azanost. Dobr´e by tak´e bylo odstranit znam´enka ze vˇsech v´yraz˚u a pokusit se nal´ezt jejich maticov´y z´apis. D´ıky tomu zvl´adneme cel´y v´ypoˇcet prov´est pˇrehlednˇe a efektivnˇe. ´ Pˇrevody sum na maticove´ nasoben´ ı Pod´ıvejme se nejprve na vyj´adˇren´ı vektoru xf . Sumu xt,f = 1 . . . n lze rozepsat dvˇema zp˚usoby: u1,f 0 ··· t∈1..n u2,f u1,f · · · X xf = bt−m+1 um,f = .. .. .. . . . m=1 ut,f ut−1,f · · · 37
Pt
m=1 bt−m+1 um,f ,
0 0 .. .
u1,f
b1 b2 .. . bt
kde t =
(5.26)
nebo
xf =
t∈1..n X m=1
bt−m+1 um,f
b1 b2 = .. .
0 b1 .. .
··· ··· .. .
u1,f 0 0 u2,f .. .. . . .
bt bt−1 · · · b1
Vektor xf (jednoho faktoru) se d´a vyj´adˇrit jako:
xf = [x1,f , x2,f , . . . , xn,f ] = U f b = Buf
(5.27)
ut,f
(5.28)
Matice U f tedy vznikla z vektoru uf pˇreskl´ad´an´ım jeho prvk˚u na diagon´aly t´eto matice. f je index faktoru, tˇechto matic je tedy celkem r. Obdobnˇe matice B vznikla z vektoru b pˇreskl´ad´an´ım prvk˚u na diagon´aly. Je vˇsak jen jedna, nebot’ i konvoluˇcn´ı j´adro b je jedno pro cel´y model. D´ale si pˇriprav´ıme vyj´adˇren´ı matice U pomoc´ı pˇr´ır˚ustkov´e matice W . Zd˚uraznˇeme zde rozd´ıl mezi maticemi U a U f . Matici U f jsme popsali v pˇredchoz´ım odstavci, kdeˇzto matice U = [u1 , u2 , . . . , ur ], obsahuje impulzn´ı retenˇcn´ı funkce faktor˚u a tedy existuje jednoznaˇcn´y pˇrevod mezi maticemi U a W . Ten nalezneme n´asleduj´ıc´ım zp˚usobem. w1,f 1 1 ... 1 1 u1,f u2,f 0 1 . . . 1 1 w2,f (5.29) .. = .. .. . . .. .. .. . . . . . . . wn,f 0 0 ... 0 1 un,f
1 ... 1 Matici ... . . . ... oznaˇc´ıme jako matici C. Pak tedy xf = Buf = BCwf a matici U lze 0 ... 1 dokonce vyj´adˇrit jako U = CW . Obdobnˇe vyj´adˇr´ıme i vektor b pomoc´ı pˇr´ır˚ustkov´eho vektoru g: b1 g1 1 1 ... 1 1 b2 0 1 . . . 1 1 g2 (5.30) .. = .. .. . . .. .. .. , . . . . . . . bn 0 0 ... 0 1 gn
z cˇ ehoˇz plyne b = Cg. Plat´ı tedy z´aroveˇn rovnost
X = U f Cg.
(5.31)
Celkovˇe pak m˚uzˇ eme vyj´adˇrit matici X jako X = BCW. Tento v´ysledek n´am hodnˇe pom˚uzˇ e v dalˇs´ıch v´ypoˇctech.
38
(5.32)
´ ˇ ych Zavislostn´ ı strom promenn ´ Pod´ıvejme se na z´avislosti jednotliv´ych veliˇcin, vytvoˇrme z´avislostn´ı strom, viz obr´azek 5.3. Deterministickou vazbu znaˇc´ıme dvojitou cˇ a´ rou, diskr´etn´ı promˇenn´e hranat´ym obvodem a spojit´e promˇenn´e kulat´ym obvodem. Do tohoto stromu zap´ısˇeme i odhadovan´e statistiky, kter´e budeme znaˇcit bez obvodu. Parametry, kter´e budeme pˇr´ımo odhadovat, se n´am ne vˇzdy budou hodit k popisu rovnice, ve kter´e je budeme vyuˇz´ıvat. V´yhodnˇejˇs´ı a pˇrehlednˇejˇs´ı bude z´apis pomoc´ı odvozen´ych promˇenn´ych, viz pˇredchoz´ı kapitolu. K tomu n´am bude slouˇzit a pom´ahat pr´avˇe z´avislostn´ı strom, kter´y n´am pom˚uzˇ e si uvˇedomit, jak´a vazba je mezi jednotliv´ymi promˇenn´ymi. D˚uleˇzitou u´ lohu pak m˚uzˇ e sehr´at pˇri samotn´em v´ypoˇctu na poˇc´ıtaˇci, d´ıky nˇemu m˚uzˇ eme nastavit spr´avn´e poˇrad´ı v´ypoˇctu dle kapitoly 3.2.4. D
A
ZA
Υ
α
β
c
ω
X
C
η
b
U
g
W
ψ
MW
ζ
h
s
ρ
ϑ
Ξ
l
κ
ν
Obr´azek 5.3: Z´avislostn´ı strom promˇenn´ych
ˇ sen´ ˇ ı integraln´ ´ ıho modelu VB metodou 5.2 Re V pˇredchoz´ı kapitole jsme zkonstruovali nov´y integr´aln´ı model pro anal´yzu scintigrafick´ych obrazov´ych sekvenc´ı, nyn´ı pˇredvedeme jeho ˇreˇsen´ı. Postupovat budeme opˇet dle kapitoly 3.2.4, tedy VB metodou.
39
5.2.1 Sestrojen´ı modelu Prvn´ı krok VB metody, sestrojen´ı matematick´eho modelu, jsme provedli jiˇz v pˇredchoz´ı kapitole. Bude vˇsak uˇziteˇcn´e si nyn´ı cel´y model pˇrehlednˇe sepsat. f (D|A, X, ω) = tN(AX ′ , ω −1Ip ⊗ In ) f (ω) = Gω (ϑ0 , ρ0 ) f (A|Υ) = tN(ZA , Ip ⊗ Υ−1 )
(5.33) (5.34) (5.35)
kde Υ = diag(v), v = [v1 , . . . , vr ]′ r Y f (v) = Gvi (αi,0 , βi,0 )
(5.36)
i=1
f (c) = tN(0r,1 , ς0 Ir ) f (g|ψ) = tN(0n , ψ −1 In ) f (ψ) = Gψ (ζ0 , η0 ) f (W |Ξ) = tN(MW , In ⊗ Ξ−1 W )
(5.37) (5.38) (5.39) (5.40)
kde Ξ = diag(ξ), ξ = [ξ1 , . . . , ξr ]′ r Y f (ξ) = Gξi (κi,0 , νi,0 )
(5.41) (5.42)
i=1
f (h) = tN(0r,1 , τ0 Ir ) f (si ) = U(0, n) f (li |si ) = U(0, n − si )
(5.43) (5.44) (5.45)
ˇ logaritmu sdruˇzeneho ´ ˇ 5.2.2 Vypo ´ cet rozdelen´ ı Budeme cht´ıt z´ıskat pˇrirozen´y logaritmus sdruˇzen´eho norm´aln´ıho rozdˇelen´ı a t´ım pˇr´ımo zjistit, zda dok´azˇ eme separovat jednotliv´e parametry. Mus´ıme proto rozepsat jednotliv´e distribuce naˇseho modelu (viz dodatek A) a vypoˇc´ıtat jejich logaritmus. Tento krok je technick´y a uv´ad´ıme ho v dodatku B. Celkovˇe dost´av´ame logaritmus sdruˇzen´eho rozdˇelen´ı
40
ln f (D, ω, A, Υ, c, g, ψ, W, Ξ, h, s, l|r) jako souˇcet vˇsech d´ılˇc´ıch logaritm˚u. ln f (D, ω, A, Υ, c, g, ψ, W, Ξ, h, s, l|r) =
pn ln ω+ 2
1 − ω tr [DD ′ − 2AX ′ D ′ + AX ′ XA′ ] + 2 r 1 pX ln υi − tr((A − ZA )Υ(A − ZA )′ )+ +(ϑ0 − 1) ln ω − ωρ0 + 2 i=1 2 r
r
X X 1 n 1 − ς0−1 tr(c′ c) + (αi,0 − 1) ln υi − υi βi,0 + ln ψ − ψ tr(g ′ g)+ 2 2 2 i=1 i=1 r nX 1 ′ +(ζ0 − 1) ln ψ − ψη0 + ln ξi − tr (W − MW )ΞW (W − MW ) + 2 i=1 2
+
r X i=1
(κi,0 − 1) ln ξi −
r X i=1
r
1 1 X 1 ξi νi,0 − τ0−1 tr(h′ h) + r ln + ln + γ, 2 n i=1 n − si
(5.46)
kde γ pˇredstavuje souˇcet veliˇcin, kter´a nez´avis´ı na parametrech sdruˇzen´eho rozdˇelen´ı.
ˇ VB-marginal ´ 5.2.3 Vypo ´ cet Podle Bayesova variaˇcn´ıho teor´emu vyj´adˇr´ıme VB-margin´aly. Pˇriprav´ıme si nˇekolik vztah˚u, kter´e budeme v VB-margin´al´ach vyuˇz´ıvat. Zaveden´ı pomocnych ´ vztahu˚ a rozkladu˚ Nejprve si rozep´ısˇeme matici D. Prvn´ı rozklad, jakoˇzto AX ′ vych´az´ı pˇr´ımo z konstrukce modelu. Pˇri vyj´adˇren´ı v dalˇs´ıch promˇenn´ych vypad´a n´asledovnˇe: AX ′ =
r X f =1
af x′f =
r X
af (U f b)′ =
r X f =1
f =1
af (Buf )′ =
r X
af (U f Cg)′ .
(5.47)
f =1
Nyn´ı si budeme cht´ıt pˇripravit vyj´adˇren´ı v´yrazu tr(XA′ AX ′ ) tak, abychom v nˇem byli schopni, po izolov´an´ı g, vyj´adˇrit ostatn´ı momentov´e odhady. Nedodrˇz´ıme zde znaˇcen´ı faktoru indexem f , nebot’ ho potˇrebujeme pouˇz´ıt dvakr´at, v´yhodnˇejˇs´ı tedy bude pracovat s i a j. Tuto
41
stopu jde rozepsat jako # "
"
tr XA′ AX ′ = tr ( = tr
"
"
r X
r X
i,j=1
r X
= tr = tr
"
i=1 r X i=1
= tr #
"
r X
#
xi (ai ′ aj )xj ′ =
i,j=1
r X
′
#
(ai ′ aj )b′ (U j U i )b =
i,j=1
# ′ (a′i aj )(U j U i ) Cg .
(5.48)
i,j=1
"
tr AX ′ D ′ = tr D ′ AX ′ = tr D ′ r X
#
"
′ (a′i aj )(U j U i ) b =
#
"
"
′
i,j=1
= tr g ′ C ′ D´ale si pˇripravme: # "
(ai ′ aj )U i bb′ U j
r X
#
aj xj ′ ) = tr
j=1
i=1
= tr b′ "
xi ai ′ )(
r X
U D ′ ai b′ = tr ′
#
C ′ U i D ′ ai g ′
#
"
"
ai xi ′ = tr D ′
i=1
#
i′
r X
r X
i′
r X i=1
#
ai b′ U
i′
#
=
U D ′ ai g ′ C ′ =
i=1
(5.49)
Protoˇze n´am vych´az´ı, zˇ e budeme potˇrebovat pˇri v´ypoˇctu druh´e momenty matice U f , najdˇeme jej´ı vhodn´y rozklad. K tomu u´ cˇ elu zavedeme n´asleduj´ıc´ı oznaˇcen´ı. Oznaˇcen´ı 4 ∆k ∈ Rn×n nazveme matici definovanou n´asledovnˇe: ( 1 i−j =k (∆k )i,j = 0 jinak
(5.50)
Lze rˇ´ıci, zˇ e matice ∆k je tedy jednotkov´a matice s diagon´alou posunutou o k n´ızˇe. Pro pˇr´ıklad uved’me: ∆0 = I a napˇr. pro n = 5 matice ∆2 vypad´a n´asledovnˇe: 0 0 0 0 0 0 0 0 0 0 ∆2 = 1 0 0 0 0 . 0 1 0 0 0 0 0 1 0 0 N´aslednˇe m˚uzˇeme nahl´ednout, zˇe ∆i ∆j = ∆i+j . Z´aroveˇn si uvˇedomme, zˇe pokud i + j ≧ n, pak ∆i+j≧n je nulov´a matice.
42
Potom tedy 1 0 U f = 0 .. . =
n−1 X
0 1 0 .. .
0 0 1 .. .
... 0 . . . 1 u + 0 1,f . . . .. .. . .
0 0 1 .. .
... . . . u2,f + . . . = . . . .. .
0 0 0 .. .
∆k uk+1,f .
(5.51)
k=0
′
Z tohoto poznatku vyjdeme pˇri rozepisov´an´ı v´yrazu U j U i : r X
′
Uj Ui =
i,j=1
r X
i,j=1
=
r X
i,j=1
n−1 X
∆′k uk+1,j
k=0
n−1 X
n−1 X
∆l ul+1,i
l=0
!
=
!
∆′k ∆l (uk+1,j ul+1,i) .
k,l=0
(5.52)
Pˇripravme si jeˇstˇe v´yhodnˇejˇs´ı z´apis pro v´ypoˇcet aprioren matic A a W . zaf (C) lze rozepsat jako 1 nebo 0 .. (5.53) zaf (C) = cf = jf cf . 1 nebo 0 a mwf (H) jako
0 .. . 0 1 . mwf (H) = .. hf = if hf . 1 0 . .. 0
(5.54)
Pro kaˇzd´y faktor jsme tedy zavedli vektory jf a if , se kter´ymi se n´am bude pracovat l´epe neˇz s mnoˇzinami C a H. D´elka nenulov´e cˇ a´ sti vektoru if je potom lf =
n X
(if )k .
k=1
T´ım m´ame pˇripraveny z´akladn´ı rozklady, kter´e vyuˇzijeme pˇri z´apisu VB-margin´al.
43
(5.55)
´ VB-marginaly h 1 ′ X)A′ )+ b ′ D ′ ) − 1 tr(A(b [ f˜(A|D) ∝ exp − ω b tr(−2AX ωX 2 2 i 1 b ′ ) − 1 tr(−2AΥ bZ cA ′ ) (5.56) − tr(AΥA 2 2 ! r n−1 h 1 X X 1b ′ ′ ′ b tr(−2 C′ ∆′k u[ D a c f˜(g|D) ∝ exp − ω k+1,i :,i g ) − ψ tr(g In g)+ 2 2 i=1 k=0 ! # " n−1 r i X X 1 ′ ′ \ − ω ad a ∆ ∆ (u u ) Cg (5.57) b tr g ′C ′ k l k+1,j l+1,i i j 2 i,j=1 k,l=0 h 1 1 ′ BCW A ′ AW ′ )+ b′ D BCW b d d f˜(W |D) ∝ exp − ω b tr(−2A )− ω b tr(C ′ B 2 2 i ′ 1 1 ′ d d d (5.58) − tr(−2ΞW MW W ) − tr(In W ΞW W ) 2" 2 r X p ˜ f (υ|D) ∝ exp + αi,0 − 1 ln υi+ 2 i i=1 # r X 1 1 1 ′ ′ A) + d b′ Z cA )i + diag(Z \ βi,0 + diag(A − υi diag(−2A i A Z A )i 2 2 2 i=1 f˜(ξ|D) ∝ exp
"
− ξi f˜(ω|D) ∝ exp
(5.59)
r X n
i=1 r X i=1
h pn
2
+ κi,0 − 1 ln ξi +
1 1 1 ′ ′W ) + c′M d \ νi,0 + diag(W diag(−2W diag(M\ i W )i + W MW )i 2 2 2
#
(5.60)
+ ϑ0 − 1 ln ω+ 2 i 1 ′ ′ ′ ′ ′ b b d [ − ω ρ0 + tr(DD − 2AX D + A AX X) 2 n 1 ′ ˜ f(ψ|D) ∝ exp + ζ0 − 1 ln ψ − ψ η0 + tr(gcg) 2 2 # " p X ′ 1 1 (jˆi )k ) − ς0−1 c2i f˜(ci |D) ∝ exp − υbi (−abi ′ jˆi ci − ci jˆi abi + c2i 2 2 k=1 ′ 1 1 ′ 2 −1 2 f˜(hi |D) ∝ exp − ξbi (−wbi iˆi hi − hi iˆi wbi + hi lˆi ) − τ0 hi 2 2 1 f˜(ji |D) ∝ exp ji − υbi (cb2i − 2b ci ai ) 2
44
(5.61) (5.62) (5.63) (5.64) (5.65)
5.2.4 Identifikace standardn´ıch forem V dalˇs´ım kroku identifikujeme VB-margin´aly jako hustoty pravdˇepodobnost´ı nˇekter´ych zn´am´ych rozdˇelen´ı. f˜(A|D) = tNA (µA , Ip ⊗ ΦA ) ˜ f(g|D) = tNg (µg , Σg ) f˜(vec(W )|D) = tNvec(W ) (µvec(W ) , Σvec(W ) ) r Y ˜ f (υ|D) = Gυi (αi , βi ) ˜ f(ξ|D) =
i=1 r Y
Gξi (κi , νi )
(5.66) (5.67) (5.68) (5.69) (5.70)
i=1
˜ f(ω|D) = Gω (ϑ, ρ) f˜(ψ|D) = Gψ (ζ, η)
(5.71) (5.72)
f˜(ci |D) = tNci (µci , Σci ) ˜ i |D) = tNh (µh , Σh ) f(h i
i
i
˜ i |D) = t Exp (λi ) f(j ji
(5.73) (5.74) (5.75)
N´asledn´e vyj´adˇren´ı tvarovac´ıch parametr˚u je ve vˇetˇsinˇe pˇr´ıpad˚u trivi´aln´ı, v pˇr´ıpadˇe matice W vˇsak provedeme detailn´ı odvozen´ı. Odvozen´ı tvarovac´ıch parametru˚ matice W Pro matici W se n´am nepodaˇrilo udrˇzet Kroneckerovskou formu z´apisu, coˇz n´as nut´ı k n´asleduj´ıc´ım u´ vah´am. Rozdˇel´ıme VB-margin´alu f˜(W |D) na dvˇe cˇ a´ sti takto: h ˜ f (W |D) ∝ exp − h exp −
i 1 1 ′d ′ ′ ′ ′ b b d ω b tr(−2A D BCW ) − ω b tr(C B BCW A AW ) × 2 2 i ′ 1 1 ′ d d d tr(−2ΞW MW W ) − tr(In W ΞW W ) 2 2
(5.76)
Z prvn´ıho cˇ lenu dost´av´ame rovnici:
coˇz znamen´a, zˇ e
−1 ′ ′ ′ BC)W A ′ AW ′ d d Σ−1 ωC ′ B W W (ΦW ) W = (b ′ BC ⊗ (A ′ A)′ )−1 , d d ΣW ⊗ ΦW = (b ωC ′ B ′ BC)−1 d ΣW = (b ωC ′B ′ (1) ′ A)−1 . d ΦW = (A (1)
45
(5.77) (5.78)
(5.79) (5.80)
Pro stˇredn´ı hodnotu t´ehoˇz cˇ lenu plat´ı ′ ′ −1 b′ D BCW b (Φ−1 bA W ) (µW ) ΣW W = ω b′ D BCΣ b (µW )′ = ω b (ΦW )′ A W
b ′ D ′ AΦ b µW = ω b (ΣW )′ C ′ B W . (1)
Z druh´eho cˇ lenu dost´av´ame rovnici:
coˇz znamen´a, zˇ e
(1)
(1)
−1 ′ ′ d ′ Σ−1 W W (ΦW ) W = In W ΞW W −1 d ΣW ⊗ ΦW = (In ⊗ Ξ W) ,
(5.86) −1
d . ΦW = Ξ W
(5.87)
′
′ ′ −1 dd (Φ−1 W ) (µW ) ΣW W = ΞW MW W
(5.88) ′
d d (µW )′ = (ΦW )′ Ξ W MW ΣW (2) (2) d d ′ (2) µ = (Σ )′ M W (ΞW ) Φ . W
(5.83)
(5.85)
(2)
Pro stˇredn´ı hodnotu t´ehoˇz cˇ lenu plat´ı
(5.82)
(5.84)
ΣW = In (2)
(5.81)
W
W
(5.89) (5.90)
˜ |D) tedy m˚uzˇ eme zapsat jako souˇcin dvou norm´aln´ıch rozdˇelen´ı releVB-margin´alu f(W vantn´ıch k prvn´ımu, resp. ke druh´emu, cˇ lenu: (1) (1) (1) (2) (2) (2) f˜(W |D) ∝ tN1 (µW , ΣW ⊗ ΦW )tN2 (µW , ΣW ⊗ ΦW ).
(5.91)
Abychom mohli pokraˇcovat v dalˇs´ıch u´ prav´ach, pˇrejdeme k vektorov´emu z´apisu t´eto distribuce: (2) (2) (2) (1) (1) (1) f˜(vec(W )|D) ∝ tN1 (vec(µW ), ΦW ⊗ ΣW )tN2 (vec(µW ), ΦW ⊗ ΣW ) (1) ′ A)′ ⊗ ω ′ BC)−1 )× d d bC ′B ∝ tN1 (vec(µ ), ((A W
−1 d × tN2 (vec(µW ), (Ξ W ⊗ In ) ). (2)
(5.92)
Sloˇzen´ım pak dostaneme distribuci tvaru
f˜(vec(W )|D) ∝ tN(µvec(W ) , Σvec(W ) ), jej´ızˇ tvarovac´ı parametry vyj´adˇr´ıme jako: −1 ′ A)′ ⊗ ω ′ BC) + (Ξ d d d Σvec(W ) = ((A bC ′B ⊗ I ) W n ′ A)′ ⊗ ω ′ BC)vec(µ(1) )+ d d µvec(W ) = Σvec(W ) ((A bC ′ B W (2) d + Σvec(W ) (ΞW ⊗ In )vec(µ ). W
46
(5.93)
(5.94)
(5.95)
´ Zapis tvarovac´ıch parametru˚ Nyn´ı jiˇz m˚uzˇ eme zapsat jednotliv´e tvarovac´ı parametry, pˇr´ıslusn´e k identifikovan´ym standardn´ım form´am. ′ X + ΥI b r )−1 [ ΦA = (b ωX b A b +Z cA Υ)Φ µA = (b ω DX
b n+ω Σg = (ψI bC ′ µg = ω b Σg C ′
Σvec(W )
r X
i,j=1
r X i=1
n−1 X k=0
(5.96)
′ ad i aj
n−1 X
k,l=0
!
\ ∆′k ∆l (uk+1,j ul+1,i) C)−1
(5.98)
!
(5.99)
′ ∆′k u[ bi k+1,i D a
′ A)′ ⊗ ω ′ BC) + (Ξ d d d = ((A bC ′ B W ⊗ In )
−1
′ A)′ ⊗ ω ′ BC) vec(µ(1) )+ d d µvec(W ) = Σvec(W ) ((A bC ′B W (2) + Σvec(W ) (In ⊗ d ΞW ) vec(µW ) p α = α0 + 1r,1 2 1 ′ d′ ) + 1 diag(−2A b′ Z cA ) + 1 diag(Z \ β = β0 + diag(AA A ZA ) 2 2 2 n κ = κ0 + 1r,1 2 1 1 1 ′ ′W ) + c′M d \ diag(−2W diag(M\ ν = ν0 + diag(W W) + W MW ) 2 2 2 pn ϑ = ϑ0 + 2 1 ′ AX ′ X) bX b ′ D ′ ) + 1 tr(A d [ ρ = ρ0 + tr(DD ′ − 2A 2 2 n ζ = ζ0 + 2 1 ′ g) η = η0 + tr(gc 2 !−1 p X σc = (jbi )k υbi + ς −1 0
i
µc i =
k=1 ′ˆ abi ji υbi σci
σhi = lbi ξbi + τ0−1
(5.97)
−1
(5.100)
(5.101) (5.102) (5.103) (5.104) (5.105) (5.106) (5.107) (5.108) (5.109) (5.110) (5.111) (5.112)
µhi = w b′ˆiξbi σhi 1 ci ai ) λi = − υbi (cb2i − 2b 2
(5.113) (5.114)
47
5.2.5 Odhad indexovych ´ mnoˇzin aprioren obrazove´ a pˇr´ırustkov ˚ e´ matice Odhad apriorna obrazov´e matice A a pˇr´ır˚ustkov´e matice W m˚uzˇ eme ch´apat jako samostatn´y vnitˇrn´ı cyklus: pod´ıv´ame-li se na obr´azek 5.3, vid´ıme, zˇ e matice A a matice W tvoˇr´ı kaˇzd´a samostatnou vˇetev. V kaˇzd´em kroku IVB algoritmu VB metody jsou ostatn´ı odhadovan´e veliˇciny zafixovan´e, m˚uzˇ eme tedy samostatnˇe odhadovat kaˇzdou vˇetev zvl´asˇt’, jakoˇzto vnitˇrn´ı cyklus IVB algortimu. Pakliˇze nechceme tato apriorna odhadovat, m˚uzˇ eme je nastavit jako nulov´e matice. Pod´ıvejme se nyn´ı na jejich tvary. Odhad apriorna A Mnoˇzina C, kter´a apriornˇe urˇcuje, zda pixel patˇr´ı k dan´emu faktoru cˇ i nikoliv, se skl´ad´a z prvk˚u 1 a 0. Jej´ı odhady jsme vˇsak v tomto tvaru neudrˇzeli. M´ısto toho n´am vyˇslo, zˇ e kaˇzd´y prvek odhadu t´eto mnoˇziny m´a exponenci´aln´ı rozdˇelen´ı. Pakliˇze ovˇsem pouˇzijeme oˇrezan´e exponenci´aln´ı rozdˇelen´ı (viz kap A.4.1) na intervalu [0, 1], bude n´am nyn´ı odhad mnoˇziny C urˇcovat pravdˇepodobnost, s jakou dan´y pixel patˇr´ı k dan´emu faktoru, jak je zn´azornˇeno na obr´azku 5.4.
P =1
P =0
0
cb
Obr´azek 5.4: Zn´azornˇen´ı v´ypoˇctu pravdˇepodobnosti pˇr´ısluˇsnosti pixelu k dan´emu faktoru
Odhad apriorna W Vˇenujme se kaˇzd´emu sloupci w matice W zvl´asˇt’, uˇsetˇr´ıme si tak psan´ı indexu f . Pro u´ pln´y popis apriorna vektoru w bylo naˇs´ım c´ılem odhadnout v´yrazy h a i, cˇ´ımˇz budeme m´ıt vektor mW jednoznaˇcnˇe reprezentovan´y v naˇs´ı pˇredpokl´adan´e parametrizaci.
48
Pro prvn´ı cˇ len jsme u´ spˇesˇnˇe vyuˇzili VB-teor´em, u druh´eho cˇ lenu vˇsak budeme postupovat jin´ym zp˚usobem. Na vektor i (skl´adaj´ıc´ı se z 1 a 0) totiˇz klademe poˇzadavek, aby jedniˇckov´e cˇ leny tvoˇrily kompaktn´ı ˇradu, obklopenou nulami. To se n´am povedlo parametrizovat a pˇrehlednˇe zapsat do cel´eho modelu pomoc´ı vektor˚u s a l. Pro odhad vektoru i ovˇsem nebudeme hledat celou distribuci f (i), ale pouze odhad jej´ıho extr´emu. Vyuˇzijeme tedy EM algoritmus: Z ˜ ln f (w, h, i) dh = ˆi = arg min f˜(w)f(h) ˜ f(h) ˜ i f(w) h ˜ = arg min Ef˜(h)f˜(w) (ln f (w, h, i)) − Ef˜(w) (ln f(w))+ i i ˜ − Ef˜(h) (ln f(h)) , (5.115)
kde
1 ′w − w c ′ i) d Ef˜(h)f˜(w) (ln f (w, h, i)) = 0 − ξbw (w b′ ib h − i′b hw b + hhi 2 −1 1 n ˜ Ef˜(w) (ln f(w)) = Ef˜(w) (ln |ξbw In |− 2 + 0) = Ef˜(w) ( ln ξbw ) = 2 n b = ln ξw 2 1 Ef˜(h) (ln f˜(h)) = Ef˜(h) ( ln(ξbw l + τ0−1 ) + 0) = 2 1 = ln(ξbw l + τ0−1 ). 2
(5.116)
(5.117)
(5.118)
Ve vyj´adˇren´ı stˇredn´ıch hodnot uv´ad´ıme i cˇ len ξbw , jenˇz m˚uzˇ eme pˇri v´ypoˇctu apriorna W fixovat a tedy v tomto v´yrazu ch´apat jako konstantu, ale i pˇrepoˇc´ıt´avat dle odhadovan´eho i.
5.2.6 Formulace VB-momentu˚ Ke kompletnosti cel´eho v´ypoˇctu n´am zb´yv´a jeˇstˇe zformulovat VB-momenty, vych´azej´ıc´ı z moment˚u standartn´ıch parametrick´ych distribuc´ı, uveden´ych v Matematick´em dodatku A. Momenty pro U \ Naˇs´ım hlavn´ım u´ kolem ted’ bude spoˇc´ıtat skal´ary u[ zit´e ve v´yrazech k+1,i a (uk+1,j ul+1,i ), pouˇ pro parametry vektoru g. Vych´azet m˚uzˇ eme z toho, zˇ e matice U je line´arn´ı transformac´ı matice W (U = CW ), jej´ızˇ stˇredn´ı hodnotu i kovarianˇcn´ı matici um´ıme vypoˇc´ıtat (ve vektorov´e podobˇe). V´ıme, zˇ e vec(W ) ∼ tN(µvec(W ) , Σvec(W ) ) a z´aroveˇn m˚uzˇ eme rozepsat U jako vec(U) = vec(CW ) = vec(CW Ir ) = (Ir ⊗ C) vec(W )
(5.119)
D˚uslednˇe zde mus´ıme rozliˇsovat znaˇcen´ı, totiˇz zˇ e matice U ∈ Rn×r je matice sloˇzen´a z vektor˚u uf , kdeˇzto matice U f ∈ Rn×n jsou cˇ tvercov´e matice pˇr´ısluˇsn´e k jednotliv´ym uf . vec(U) pak m´a rozdˇelen´ı vec(U) = tN((Ir ⊗ C)µvec(W ) , (Ir ⊗ C)Σvec(W ) (Ir ⊗ C)′ ).
49
(5.120)
Oznaˇcme µvec(U ) = (Ir ⊗ C)µvec(W )
(5.121)
Σvec(U ) = (Ir ⊗ C)Σvec(W ) (Ir ⊗ C)′ . \ Hledan´e skal´ary (uk+1,j ul+1,i ) z´ısk´ame z v´yrazu u1,1 u2,1 .. . un,1 u u1,1 u2,1 . . . un,1 E(vec(U) vec(U)′ ) = E( 1,2 ... u n,2 . .. un,r 2 E(u1,1) E(u1,1 u2,1) . . . E(u1,1 un,1) E(u1,1 u1,2) E(u2,1u1,1 ) . . . E(u2,1 un,1) E(u2,1 u1,2) E(u22,1 ) .. .. .. .. .. . . . . . 2 E(u u ) E(u u ) . . . E(u ) E(u n,1 1,1 n,1 2,1 n,1 u1,2 ) n,1 E(u21,2) E(u1,2u1,1 ) E(u1,2 u2,1) . . . E(u1,2 un,1) .. .. .. .. .. . . . . . E(un,r u1,1 ) E(un,r u2,1) . . . E(un,r un,1) E(un,r u1,2)
(5.122)
u1,2 . . . un,r ) =
. . . E(u1,1 un,r ) . . . E(u2,1 un,r ) .. .. . . . . . E(un,1un,r ) . . . E(u1,2 un,r ) .. .. . . ...
E(u2n,r )
= Σvec(U ) + µvec(U ) (µvec(U ) )′ ,
(5.123)
kde Σvec(U ) ∈ Rnr×nr a µvec(U ) ∈ Rnr×1 . Tedy vec(U ) vec(U ) vec(U ) \ (uk+1,j ul+1,i) = Σ(j−1)n+(k+1),(i−1)n+(l+1) + µ(j−1)n+(k+1) µ(i−1)n+(l+1) .
(5.124)
T´ım m´ame vyj´adˇreny vˇsechny stˇredn´ı hodnoty kombinac´ı prvk˚u matice U. Moment u[ adˇr´ıme jako k+1,i vyj´ vec(U )
u[ k+1,i = µ(i−1)n+(k+1) .
(5.125)
µb = Cµg Σb = CΣg C ′
(5.126) (5.127)
Momenty pro B Nejprve si vyj´adˇr´ıme moment bb.
50
Tedy bb = µb = Cµg .
(5.128)
′ B. baB d Obdobnˇe jako u matice U si pˇriprav´ıme vyj´adˇren´ı pro dalˇs´ı momenty, B Prvn´ı v´yraz vyj´adˇr´ıme snadno jako
b= B
′
n−1 X k=0
∆k bd k+1 .
(5.129)
U stˇredn´ı hodnoty B B budeme postupovat obdobnˇe jako v pˇr´ıpadˇe matic U i U j . B = ∆0 b1 + ∆1 b2 + · · · + ∆n−1 bn = Potom tedy B′B =
n−1 X k=0
∆k bk+1
!′
n−1 X l=0
∆l bl+1
!
=
n−1 X
n−1 X
∆k bk+1
(5.130)
k=0
∆k ′ ∆l (bk+1 bl+1 ).
(5.131)
k,l=0
Nyn´ı chceme vyj´adˇrit cˇ leny bd e z´ısk´ame z n´asleduj´ıc´ıho. Je-li i bj , kter´ b ∼ tN(µb , Σb ), pak E(b21 ) E(b1 b2 ) . . . E(b1 bn ) E(b2 b1 ) E(b2 ) . . . E(b2 bn ) 2 ′ E(bb′ ) = .. .. = Σb + µb µb . .. .. . . . . E(bn b1 ) E(bn b2 ) . . . E(b2n )
(5.132)
Z t´eto matice pak jiˇz snadno vybereme kombinaci dvou hledan´ych prvk˚u: \ bk+1 bl+1 = (Σb )k+1,l+1 + (µb )k+1(µb )l+1 .
(5.133)
Momenty pro X Vyj´adˇren´ı odhadu stˇredn´ı hodnoty X je trivi´aln´ı, nebot’ staˇc´ı pouˇz´ıt vzorec b = BC b W c, X
(5.134)
c z´ısk´ame zpˇetn´ym poskl´ad´an´ım v´yrazu µvec(W ) do matice o pˇr´ısluˇsn´ych rozpˇriˇcemˇz cˇ len W mˇerech. ′ B urˇ d Pro stˇredn´ı hodnotu X ′ X vyuˇzijeme moˇznost vyj´adˇrit X ′ X jako W ′C ′ B ′ BCW . B cit um´ıme z pˇredchoz´ıho, tedy ′ ′d ′X = E ′ [ X f (W |D) (W C B BCW ).
(5.135)
′ ZW (kde Z je C ′ B ′ BC), pˇriˇ d \ Probl´em se redukoval na urˇcen´ı matice tvaru W cemˇz zn´ame Σvec(W ) . Probl´em budeme ˇreˇsit algoritmicky. Tento algoritmus proch´az´ı jednotliv´e bloky (n × n) kovarianˇcn´ı matice Σvec(W ) (nr × nr), vyn´asob´ı je matic´ı Z. Pot´e spoˇc´ıt´ame stopu v´ysledn´e matice a v´ysledek uloˇz´ıme na pˇr´ısluˇsnou pozici pomocn´e matice (r×r). V´yslednou pomocnou c′Z W c. matici pak pˇriˇcteme k matici W
51
Ostatn´ı momenty Dalˇs´ı momenty jsou jiˇz trivi´aln´ı, z dodatku A.3 plyne: b = diag(b Υ υ ) = diag(α ◦ β −1 ) b = diag(q) = diag(κ ◦ ν −1 ) Ξ θ ω b= ρ ζ ψb = , η
(5.136) (5.137) (5.138) (5.139)
kde operace ”◦” znaˇc´ı Hadamard˚uv souˇcin, viz dodatek A.1. Potˇrebn´e momenty oˇrezan´ych ′ A, b bA d c′, hbi , hb2i , cbi , cb2i dopoˇcteme snadno dle dodatku A.2 a monorm´aln´ıch rozdˇelen´ı A, g , gg ment jbi dle A.4.1 ∀i ∈ 1 . . . r.
ˇ 5.3 Shrnut´ı vypo ´ ctu
Seps´an´ım VB-moment˚u jsme dokonˇcili v´ypoˇcet integr´aln´ıho modelu VB-metodou. Ostatn´ı kroky jsou technick´e a v´ypoˇcetn´ı, budeme je tedy uv´adˇet v n´asleduj´ıc´ım textu, jen bude-li to potˇreba. V´ypoˇcet se n´am pomoc´ı pomocn´ych promˇenn´ych a odvozen´ych vztah˚u podaˇrilo udrˇzet v kompaktn´ı podobˇe. Z´aroveˇn jsme dok´azali zahrnout do v´ypoˇctu oba diagnostick´e parametry, kter´ym se tato pr´ace vˇenuje. Pˇredt´ım, neˇz v dalˇs´ı kapitole pˇredvedeme v´ysledky tohoto integr´aln´ıho modelu, se pod´ıv´ame na dvˇe dalˇs´ı nastaven´ı v´ypoˇctu: inicializaci a konvergenci.
ˇ 5.3.1 Inicializace vypo ´ ctu Teoreticky m˚uzˇ eme startovat z n´ahodn´ych poˇca´ teˇcn´ıch dat. To vˇsak nen´ı pˇr´ıliˇs v´yhodn´e, nebot’ by se mohl v´ypoˇcet zastavit v nˇejak´em lok´aln´ım extr´emu a t´ım bychom nedoˇsli ke spr´avn´ym v´ysledk˚um. V´yhodnˇejˇs´ı je proto nastavit hodnoty, kter´e dok´azˇ eme dobˇre odhadnout, a iterace zah´ajit v´ypoˇctem tˇech, kter´e nezn´ame ani neodhadneme. Pod´ıv´ame-li se na obr´azek 5.3, vid´ıme, zˇ e nastav´ıme-li vˇsechny hodnoty kromˇe vˇetve matice A, tato se n´am odhadne jako prvn´ı a budeme-li m´ıt dobr´e poˇca´ teˇcn´ı odhady, v´ypoˇcet rychle dokonverguje k v´ysledku. Poˇca´ teˇcn´ı hodnoty b a W pˇritom nen´ı tˇezˇ k´e urˇcit. b m´a tvar exponenciely a W jsou pˇr´ır˚ustky impulzn´ıch retenˇcn´ıch funkc´ı, jejichˇz teoretick´e tvary jsou zn´am´e, takˇze je t´ezˇ m˚uzˇ eme dobˇre odhadnout.
ˇ 5.3.2 Konvergence vypo ´ ctu Konvergence v´ypoˇctu nen´ı v t´eto pr´aci ˇreˇsena teoreticky. Vol´ıme proto pevn´y poˇcet krok˚u IVB algoritmu (typicky 15 aˇz 50). Sledovat vˇsak m˚uzˇ eme pomocn´e ukazatele, jako napˇr. diagon´aln´ı matici Υ, jej´ızˇ hodnoty by mˇely klesat a podle nichˇz m˚uzˇ eme poznat tendenci konvergence jednotliv´ych faktor˚u. Pr´ace si vˇsak nevyˇza´ dala automatick´e odhadov´an´ı t´eto tendence.
52
´ ıho modelu 6 Vysledky ´ integraln´ Pod´ıvejme se nyn´ı na v´ysledky, jak´ych lze dos´ahnout pomoc´ı integr´aln´ıho modelu, zkonstruovan´eho a vyˇreˇsen´eho v kapitole 4. Obdobnˇe jako u standardn´ıho modelu bude c´ılem urˇcit diagnostick´e koeficienty, tranzitn´ı cˇ as (TT) a relativn´ı ren´aln´ı clearanci. V´ysledn´y algoritmus je opˇet implementov´an v Matlabu. Fungov´an´ı a chov´an´ı integr´aln´ıho modelu si uk´azˇ eme na sekvenci sn´ımk˚u z obr´azku 4.2, tedy kaˇzd´y sn´ımek m´a rozmˇery 64 × 64, z cˇ ehoˇz plyne, zˇ e p = 4096. Sn´ımk˚u je celkem 120 (n = 120) a jsou sn´ım´any vˇzdy po 10s. V´ysledek cel´e anal´yzy vid´ıme na obr´azku 6.1. Popisu v´ysledk˚u se vˇenuje n´asleduj´ıc´ı kapitola. 6
20 10
10
20
20
30
30
40
40
50
50
60 40
60
6
x 10
2
6
1.5
4
1
10 20 30
40
40
50
40
60
40
60
5
0
50
100
150
0
3
0
50
100
150
10 20
30
30
40
40
50
50
60
60 20
40
60
0
15
10 20 30
40
40
50
50
60 40
60
150
1
0
50
100
150
0
6
3
0
50
100
150
50
100
150
50
100
150
50
100
150
4
x 10
15
x 10
10
2
10
1.5
1
5
1
5
2 0.5 40
60
1
0
50
100
150
0
0.5 0
50
100
150
0
0
50
100
150
0
6
x 10
8
8
2
6
6
1.5
4
4
1
0
50
100
150
0
6
x 10
6
0 5
x 10
15
x 10
5 4
10
3 2
2
2
0.5
0
0
0
5
1 40
60
0
50
100
150
0
50
100
150
15
0
50
100
150
0
4
x 10
10
1.5 10
0
50
100
150
10
8
8
6
6
4
4
2
2
40
60
0
0 4
x 10
8
x 10
6 4
5 0.5
20
0
5
x 10
1
60 20
100
2.5
2
2
30
50
x 10
1.5
5
20
0
2.5
3
20
10
0.5
4
x 10
6
20
1
0.5
4
20
10
1.5
1.5
60 20
x 10
2 2
50
60
2
3
6
30
6
x 10
2.5
10
5
20
4
15
20
10
6
x 10
3.5
60 20
8
0
50
100
150
0
0
50
100
150
0
0
50
100
150
0
2
0
50
100
150
0
0
Obr´azek 6.1: v´ysledek anal´yzy studie ledvin pomoc´ı integr´aln´ıho modelu
6.1 Popis odhadovanych ´ parametru˚ Nyn´ı vysvˇetl´ıme, co znamenaj´ı jednotliv´e hodnoty, a n´aslednˇe uk´azˇ eme, jak z nich vyˇc´ıst poˇzadovan´e klinick´e parametry.
53
Kaˇzd´y ˇra´ dek na obr´azku 6.1 se t´yk´a jednoho nalezen´eho faktoru. Popiˇsme nyn´ı postupnˇe, jak´e hodnoty jsou v jednotliv´ych sloupc´ıch: 1. odhad mnoˇziny C, ovlivˇnuj´ıc´ı apriorno matice obr´azk˚u A; jedn´a se o pravdˇepodobnostn´ı masku, kter´a urˇcuje, s jakou pravdˇepodobnost´ı patˇr´ı dan´y pixel k tomuto faktoru, prvky t´eto masky jsou tedy v intervalu [0; 1], 2. obr´azek konkr´etn´ıho faktoru, tedy pˇreskl´adan´y dan´y sloupec matice A, 3. kˇrivka konkr´etn´ıho faktoru, tzn. dan´y sloupec matice X, 4. odhadnut´a impulzn´ı retenˇcn´ı funkce dan´eho faktoru, sloupec matice U, 5. odhadnut´e pˇr´ır˚ustky impulzn´ı retenˇcn´ı funkce dan´eho faktoru, sloupec matice W , 6. odhadovan´a apriorn´ı impulzn´ı retenˇcn´ı funkce dan´eho faktoru 7. odhadovan´e apriorn´ı pˇr´ır˚ustky impulzn´ı retenˇcn´ı funkce dan´eho faktoru, sloupec matice MW . Podaˇrilo se n´am odhadnout stejnˇe jako pomoc´ı standardn´ıho modelu cˇ tyˇri faktory: moˇcov´y mˇech´yˇr, p´anviˇcky, krevn´ı pozad´ı a parenchymy. M˚uzˇ eme si vˇsak vˇsimnout, zˇ e odhadnut´e kˇrivky se liˇs´ı. Konkr´etnˇe, u integr´aln´ıho modelu nepˇredpokl´ad´ame, zˇ e m˚uzˇ e b´yt kˇrivka faktoru na poˇca´ tku nulov´a, naˇse konstrukce impulzn´ı retenˇcn´ı funkce to neumoˇznˇ uje. D˚usledkem toho je tak´e pˇr´ıliˇsn´e odeˇcten´ı p´anviˇcek od parenchymu. To je nev´yhoda tohoto modelu a je to moˇzn´y n´amˇet pro pokraˇcov´an´ı. Pro anal´yzu poˇca´ teˇcn´ı f´aze sekvence a urˇcen´ı relativn´ı ren´aln´ı clearance to vˇsak nevad´ı, model je konstruov´an pˇredevˇs´ım pro ni. Jinak vˇsak model odhaduje parametry, kter´e jsme po nˇem poˇzadovali. M˚uzˇ eme tedy konstatovat, zˇ e funguje tak, jak jsme chtˇeli a jak jsme ho navrhli. 4
8 10
12
10
20
20
30
30
40
40
50
50
60 40
60
8
4
0
5
10
15
20
1.5 10 20
30
30
40
40
50
50
60
60 20
40
60
6.8994
6
6
1
0.5
20
40
60
0
4
6.8994 2
2
7 0
0
5
10
15
20
x 10
6.8995
9
2
60
8 6
8
40
4
x 10
6.8995
6.8994 0
5
10
15
20
0
4
20
6.8996
10
20
10
4
x 10
6
4
60 20
4
x 10
11 6
0
5
10
15
20
6.8993
4
x 10
6
0
5
10
15
20
0
4
x 10
6
6
5
5
5
5
4
4
4
4
3
3
3
3
2
2
2
2
1
1
1
0
0
0
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
4
x 10
x 10
1 0
5
10
15
20
0
−3
1.5 10
10
20
20
30
30
40
40
50
50
60 40
60
40
20
40
20
30
15
30
15
20
10
20
10
10
5
10
1
0.5
60 20
x 10
20
40
60
0
0
5
10
15
20
0
0
5
10
15
20
0
0
5
10
15
20
0
5
0
5
10
15
20
0
Obr´azek 6.2: v´ysledek anal´yzy poˇca´ teˇcn´ı f´aze studie ledvin pomoc´ı integr´aln´ıho modelu
54
ˇ 6.1.1 Tranzitn´ı cas Spr´avn´y odhad tvaru impulzn´ı retenˇcn´ı funkce, pˇresnˇe dle teoretick´eho, n´am umoˇznˇ uje pˇr´ımo z dosaˇzen´ych v´ysledk˚u vyˇc´ıst TT parenchymu. Nemus´ıme tedy tuto funkci dopoˇc´ıt´avat pomoc´ı Fourierovy transformace z odhadnut´ych kˇrivek a t´ım zan´est do v´ysledku chybu, jako to bylo u standardn´ıho modelu v kapitole 4.3.2. par Z ˇra´ dku parenchymu, z sˇest´eho sloupce, pˇr´ımo vyˇcteme, zˇ e T Tmin = 160s.
´ ı clearance 6.1.2 Relativn´ı renaln´ Pro urˇcen´ı relativn´ı ren´aln´ı clearance vyuˇzijeme opˇet jen prvn´ı cˇ a´ st cel´e sekvence, nebot’ na t´eto cˇ a´ sti n´asˇ model splˇnuje z´akladn´ı biologick´e pˇredpoklady, odpad´a probl´em s p´anviˇckou a moˇcov´ym mˇech´yˇrem, kter´e jsou na zaˇca´ tku nulov´e. Ot´azku, kolik br´at prvn´ıch sn´ımk˚u, vyˇreˇs´ıme automaticky, udˇel´ame prvn´ı anal´yzu na pln´e sekvenci a z faktoru parenchymu vezmepar me hodnotu T Tmin . V dalˇs´ım kroku analyzujeme tedy jen poˇca´ teˇcn´ı f´azi sekvence a pomoc´ı korelace pˇredpokl´adan´eho tvaru ledvin s pravdˇepodobnostn´ımi maskami (odhad mnoˇziny C), pˇr´ıpadnˇe podle TT, urˇc´ıme, kter´y faktor pˇredstavuje parenchym. V´ysledek vid´ıme na obr´azku 6.2. Anal´yzou poˇca´ teˇcn´ı f´aze se podaˇrilo rozliˇsit parenchym, krevn´ı pozad´ı a dalˇs´ı tk´anˇ ov´y faktor. Tento posledn´ı rozliˇsen´y faktor je velice slab´y, proto nen´ı vidˇet na v´ysledku celkov´e anal´yzy, kde je pod rozliˇsovac´ı schopnost´ı sˇ umu. Zde ho vˇsak jeˇstˇe nalezneme, coˇz je d˚uleˇzit´e z hlediska spr´avn´eho odeˇcten´ı pozad´ı od parenchymu, coˇz se podaˇrilo.
10
20
30
40
50
60 10
20
30
40
50
60
Obr´azek 6.3: vyˇciˇstˇen´y sn´ımek parenchym˚u, z´ıskan´y pomoc´ı integr´aln´ıho modelu
M´ame tedy velice pˇekn´y sn´ımek parenchymu a jeˇstˇe k tomu pravdˇepodobnostn´ı masku, kter´a urˇcuje, zda dan´y pixel patˇr´ı k faktoru, cˇ i nikoliv. Pro v´yslednou anal´yzu clearance je vhodn´e tuto masku prahovat hodnotou napˇr. 0, 5 na diskr´etn´ı hodnoty 1 a 0. Touto maskou
55
pak obr´azek parenchymu po prvc´ıch vyn´asob´ıme a z´ısk´ame tak oˇcek´avanou aktivitu org´anu, coˇz je prakticky vyˇciˇstˇen´y faktorov´y obr´azek, viz obr´azek 6.3. Tato maska n´am slouˇz´ı jako automatick´e ROI. Poznamenejme, zˇ e tento krok nezapad´a do celkov´eho algoritmu svoj´ı povahou, nicm´enˇe kompenzujeme tak nemoˇznost vytvoˇrit apriorno dle pozn´amky 3 kv˚uli v´ypoˇcetn´ı n´aroˇcnosti. V takto z´ıskan´em v´ysledku jiˇz jen staˇc´ı seˇc´ıst pixely prav´eho a lev´eho parenchymu a vypoˇc´ıtat relativn´ı ren´aln´ı clearance (v tomto pˇr´ıpadˇe je v´ysledek ParL ≈ 49% a ParP ≈ 51% (lev´a a prav´a strana je br´ana z naˇseho pohledu)). V t´eto oblasti medic´ıny neexistuje standard (viz [10]), nem´ame proto s cˇ´ım tyto hodnoty porovnat. M˚uzˇ eme vˇsak vizu´alnˇe odhadnout, zˇ e se bl´ızˇ´ı pravdˇe. Pro d˚ukladn´e srovn´an´ı bude tˇreba prov´est rozs´ahlou studii, kter´a pˇresahuje r´amec t´eto pr´ace. Dok´azali jsme tedy zcela automaticky analyzovat relativn´ı ren´aln´ı clearance. Nutno vˇsak dodat, zˇ e tento postup pˇredpokl´ad´a pˇr´ıtomnost obou ledvin, nepˇredpokl´ad´a jejich sr˚ust apod.
ˇ ıch studi´ı 6.2 Vyhodnocen´ı dals´ Zkusme na tomto m´ıstˇe vyhodnostit nˇekolik dalˇs´ıch studi´ı. Ty uˇz budeme volit ne tak typicky jako v pˇredchoz´ım experimentu, abychom se pod´ıvali, jak se nov´y integr´aln´ı model chov´a v tˇechto pˇr´ıpadech. Protoˇze pomoc´ı standardn´ıho modelu nelze urˇcit relativn´ı ren´aln´ı clearance, kterou nyn´ı budeme diagnostikovat, nebudeme jiˇz jeho v´ysledky uv´adˇet. 100 10
10
20
20
30
30
40
40
40
50
50
20
60
60 20
40
60
200
80 60
100
20
40
60
0
10
20
20
30
30
40
40
50
50
60
60 20
40
60
0
5
10
15
20
10
20
20
30
30
40
40
50
50
60
60 20
40
60
30
60
50
50
40
40
40
30
30
30
20
20
20
10
10
0
0
0
0
5
10
15
20
0
5
10
15
20
10 0
5
10
15
20
0
100
50
50
50
80
40
40
40
60
30
30
30
40
20
20
20
20
10
10
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
20 10
20
10
60
50
50
40 10
60
150
20
40
40
60
60
0
0
5
10
15
20
0
50
400
40
300
30
200
20
100
10
0
0
5
10
15
20
0
0
5
5
10
10
15
15
20
20
0
0
5
10
15
20
0
10
0
5
10
15
20
0
300
300
300
250
250
250
200
200
200
150
150
150
100
100
100
50
50
0
0
0
5
10
15
20
50 0
5
10
15
20
0
Obr´azek 6.4: v´ysledky studie 1
6.2.1 Studie 1 Sekvenci sn´ımk˚u prvn´ı studie m˚uzˇ eme vidˇet v doplˇnku na obr´azku C.1. Jedn´a se o znaˇcnˇe poˇskozenou ledvinu, kdy na prvn´ı pohled funguje prav´y parenchym (opˇet z naˇseho pohledu)
56
sˇpatnˇe. Kaˇzd´y sn´ımek sekvence m´a rozmˇery 64 × 64, tzn. p = 4096. Ze sekvence vezmeme prvn´ıch 20 sn´ımk˚u (n = 20), ty jsou sn´ım´any vˇzdy po 10s. V´ysledek cel´e anal´yzy vid´ıme na obr´azku 6.4. Niˇzsˇ´ı aktivita krve na poˇca´ tku m˚uzˇ e b´yt zp˚usobena t´ım, zˇ e byla aplikov´ana kontrastn´ı l´atka pˇr´ıliˇs pozdˇe, maximum aktivity krve tedy nen´ı hned na prvn´ım sn´ımku. Pro v´ypoˇcet to m˚uzˇ e pˇredstavovat probl´em a ˇreˇsen´ım je tento prvn´ı sn´ımek vyˇradit z anal´yzy, zde to vˇsak nebylo nutn´e. Pro urˇcen´ı clearance pouˇzijeme pravdˇepodobnostn´ı masku v prvn´ım sloupci. Faktor parenchym˚u s odstranˇen´ym sˇumem touto metodou vid´ıme na obr´azku 6.5. Po vyhodnocen´ı tohoto sn´ımku dost´av´ame hodnoty ParL ≈ 80, 5% a ParP ≈ 19, 5% (lev´a a prav´a strana je br´ana z naˇseho pohledu).
10
20
30
40
50
60 10
20
30
40
50
60
Obr´azek 6.5: vyˇciˇstˇen´y sn´ımek parenchym˚u, z´ıskan´y pomoc´ı integr´aln´ıho modelu ve studii 1
6.2.2 Studie 2 Sekvence druh´e studie je v doplˇnku na obr´azku C.2. M˚uzˇ eme si vˇsimnout, zˇ e zde bylo zapoˇcato sn´ım´an´ı scintigrafickou kamerou pˇr´ımo ve chv´ıli aplikace radiofarmaka. Pˇred anal´yzou takov´e sekvence proto pouˇzijeme oˇr´ıznut´ı cˇ a´ st´ı, kter´e si na sekvenci nepˇrejeme, a data uprav´ıme, aby splˇnovala naˇse pˇredpoklady. Studie zobrazuje novorozenˇe se sˇpatnˇe funguj´ıc´ı pravou ledvinou. Kaˇzd´y sn´ımek sekvence m´a rozmˇery 64 × 64, tzn. p = 4096. Ze sekvence vezmeme prvn´ıch 16 sn´ımk˚u (n = 16), jsou sn´ım´any vˇzdy po 10s. V´ysledek anal´yzy s pouˇzit´ım integr´aln´ıho modelu je na obr´azku 6.6. Pro urˇcen´ı clearance pouˇzijeme opˇet pravdˇepodobnostn´ı masku v prvn´ım sloupci. Faktor parenchym˚u s odstranˇen´ym sˇumem touto metodou vid´ıme na obr´azku 6.7. Po vyhodnocen´ı tohoto sn´ımku dost´av´ame hodnoty ParL ≈ 68% a ParP ≈ 32% (lev´a a prav´a strana je br´ana z naˇseho pohledu).
57
2500
4000 10
10
20
20
30
30
40
40
50
50
60
800
2000
3000
1500
800
600
600 1000
1500 2000
400
400
1000 1000
500 200
500
200
60 20
40
60
20
10
10
20
20
30
30
40
40
50
50
60
40
60
60 20
40
60
20
40
60
0
0
5
10
15
20
15
20
0
0
5
10
15
20
15
20
0
400
300
300
200
200
200
200
200
100
100
100
100
0
250 200
0
5
10
15
20
40
40
50
50
0
0
5
10
15
20
0
110
80
100
60
90
40
0
5
10
15
20
5
10
15
20
0
5
10
15
20
0
5
10
15
20
100
0
5
10
15
20
74.5
0
80
74
60 40
73
100 80
20
0
5
10
15
20
70
0
5
10
15
20
0
20
72.5
50 0
0
0
73.5
150
60
10
400
30
40
5
300
20
20
0
400
10
60
0
300
30
60
10
400
20
40
5
300
10
20
0
400
300
60
0
0
5
10
15
20
72
0
5
10
15
20
0
Obr´azek 6.6: v´ysledky studie 2
6.3 Zhodnocen´ı vysledk ´ u˚ V t´eto kapitole jsme na nˇekolika studi´ıch pˇredvedli fungov´an´ı integr´aln´ıho modelu pro anal´yzu scintigrafick´ych sekvenc´ı. Detailnˇe jsme popsali odhadovan´e veliˇciny a mechanismus, jak z´ıskat hledan´e diagnostick´e parametry. Z v´ysledk˚u vid´ıme, zˇ e model jeˇstˇe nefuguje pro celou sekvenci, ale jen pro jej´ı poˇca´ tek, kde jsou splnˇeny jeho pˇredpoklady. Naˇs´ım u´ kolem vˇsak bylo zkoumat a diagnostikovat parametr, kter´y z´avis´ı pr´avˇe pouze na poˇca´ tku a kter´y nen´ı ovlivnˇen n´asledn´ymi sn´ımky. Z´ıskali jsme tak automatickou metodu pro odhad relativn´ı ren´aln´ı clearance, kter´a jeˇstˇe nefunguje zcela obecnˇe, protoˇze zkoum´ame velice komplikovan´y syst´em, lidsk´e tˇelo, ale kter´a ˇreˇs´ı pˇredloˇzen´y probl´em.
58
10
20
30
40
50
60 10
20
30
40
50
60
Obr´azek 6.7: vyˇciˇstˇen´y sn´ımek parenchym˚u, z´ıskan´y pomoc´ı integr´aln´ıho modelu ve studii 2
59
´ er ˇ a moˇznosti dals´ ˇ ıho 7 Zav ˇ an´ ´ ı pokracov C´ılem pr´ace bylo prozkoumat metody zpracov´an´ı dat a odhadov´an´ı nˇekter´ych d˚uleˇzit´ych diagnostick´ych koeficient˚u v nukle´arn´ı medic´ınˇe. D´ale pak navrhnout model, kter´y by respektoval biologick´a omezen´ı a pomoc´ı nˇehoˇz bychom pˇr´ımo odhadovali potˇrebn´e hodnoty, coˇz se doposud v klinick´e praxi nedˇeje. K v´ypoˇctu tohoto modelu jsme vyuˇzili Variaˇcn´ı Bayes˚uv teor´em, zabudovan´y ve VB metodˇe. V´ysledn´y algoritmus jsme otestovali na re´aln´ych scintigrafick´ych studi´ıch.
´ 7.1 Hlavn´ı pˇr´ınos prace Kapitola 2 Sezn´amili jsme se se z´akladn´ı metodou, pouˇz´ıvanou v nukle´arn´ı medic´ınˇe, scintigrafi´ı. Podrobnˇe jsme prozkoumali fungov´an´ı t´eto metody, jej´ı v´yhody, nev´yhody a moˇzn´a u´ skal´ı. Detailnˇe jsme se pod´ıvali na radionuklidy pouˇz´ıvan´e touto metodou, s jejichˇz vlastnostmi je pot´e v matematick´em modelu nutno poˇc´ıtat. Velk´a cˇ a´ st t´eto kapitoly je pak vˇenov´ana matematick´emu pohledu na celou anal´yzu, jsou zm´ınˇeny z´akladn´ı biologick´e pˇredpoklady, matematick´e metody a vysvˇetleny diagnostick´e koeficienty, kter´ych se tato pr´ace t´yk´a. Kapitola 3 Probrali jsme z´aklady Bayesovsk´e statistiky a pouk´azali na nˇekter´e jej´ı v´yhody a nev´yhody. Mezi hlavn´ı v´yhody patˇr´ı moˇznost efektivnˇe zahrnout apriorn´ı informaci do v´ypoˇctu a tak se dostat k lepˇs´ım v´ysledk˚um, neˇz pomoc´ı klasick´e statistiky. Plat´ıme za to ovˇsem vysokou v´ypoˇcetn´ı n´aroˇcnost´ı, kter´a n´as nut´ı vyuˇz´ıvat aproximace. Tˇem je vˇenov´ana druh´a cˇ a´ st kapitoly, pˇredevˇs´ım pak aproximaci pomoc´ı Bayesova Variaˇcn´ıho teor´emu. Jeho pouˇzit´ı shrnujeme do 8 z´akladn´ıch krok˚u, tzv. VB metody. Kapitola 4 Vˇenovali jsme se konstrukci standardn´ıho modelu pro anal´yzu scintografick´ych obrazov´ych sekvenc´ı pomoc´ı Variaˇcn´ıho Bayesova teor´emu. Tento model jsme n´aslednˇe vyˇreˇsili VB metodou. V´ysledky jsme uk´azali na re´aln´e studii. Jako probl´em pˇri n´asledn´em urˇcov´an´ı diagnostick´ych parametr˚u se uk´azala nedokonal´a separace jednotliv´ych faktor˚u a nepˇresn´e urˇcen´ı kˇrivek faktor˚u, coˇz vedlo k nemoˇznosti automatick´eho urˇcen´ı n´ami poˇzadovan´ych parametr˚u.
60
Kapitola 5 P´at´a kapitole je hlavn´ım pˇr´ınosem t´eto pr´ace. Detailnˇe jsme v n´ı popsali konstrukci matematick´eho modelu, kter´y respektuje z´akladn´ı biologick´e poˇzadavky, pˇredevˇs´ım konvoluci kˇrivek. Zkonstruovali a zavedli jsme tak´e apriorna, kter´a n´am pomohla automaticky bˇehem v´ypoˇctu odhadovat poˇzadovan´e diagnostick´e koeficienty. To vedlo k znaˇcn´emu zesloˇzitˇen´ı cel´eho modelu, odvodili jsme proto z´akladn´ı metody udrˇzen´ı v´ypoˇctu v maticov´em tvaru. Model jsme vyˇreˇsili VB metodou, v pr˚ubˇehu v´ypoˇctu bylo nutno s´ahnout i k EM algoritmu. Kapitola 6 ˇ a kapitola ukazuje diagnostick´e moˇznosti nov´eho integr´aln´ıho modelu. Rozebrali jsme Sest´ jednotliv´e odhadovan´e veliˇciny a podrobnˇe uk´azali, jak pomoc´ı nich z´ısk´ame hledan´e diagnostick´e parametry. Na nˇekolika scintigrafick´ych studi´ıch jsme pot´e pˇredvedli fungov´an´ı modelu a u kaˇzd´e studie automaticky urˇcili relativn´ı ren´aln´ı clearanci parenchym˚u.
ˇ ıho pokracov ˇ an´ ´ ı 7.2 Moˇznosti dals´ Jiˇz bˇehem pr´ace jsme naznaˇcili nˇekolik smˇer˚u, kter´ymi by bylo dobr´e v budoucnu pokraˇcovat. Vypiˇsme tu na z´avˇer nˇekolik nejperspektivnˇejˇs´ıch. Prvn´ı zejdnoduˇsen´ı jsme udˇelali v urˇcov´an´ı sˇumu, kter´y modelujeme jako stejn´y ve vˇsech cˇ a´ stech obr´azku. To v praxi neplat´ı, jsou oblasti s menˇs´ım a vˇetˇs´ım sˇumem. Zm´ınili jsme se tak´e, zˇ e metoda vylouˇc´ı nˇekter´e faktory, kter´e jsou vidˇet napˇr. jen na poˇca´ tku sekvence a jsou velice m´alo v´yrazn´e, urˇc´ı je jakoˇzto sˇum. Tento probl´em jsme prozat´ım odstranili anal´yzou pouze poˇca´ teˇcn´ı f´aze sekvence, obecnˇe je ale sˇum v nukle´arn´ı medic´ınˇe a nejen v n´ı velk´y probl´em a st´alo by za to se mu vˇenovat v´ıce. Nutn´e zlepˇsen´ı bude tak´e u modelov´an´ı impulzn´ıch retenˇcn´ıch funkc´ı, kter´e jsme zjednodusˇili pˇr´ıliˇs, coˇz mˇelo za n´asledek sˇpatnou separaci p´anviˇcky od parenchymu pˇri anal´yze nezkr´acen´e sekvence. N´asˇ model nyn´ı nepoˇc´ıt´a s t´ım, zˇ e poˇca´ teˇcn´ı f´aze faktoru m˚uzˇ e b´yt nulov´a. Zaveden´ım dalˇs´ı podm´ınky do apriorna pˇr´ır˚ustkov´e matice by mˇelo j´ıt tento probl´em odstranit a dostat se tak k lepˇs´ım v´ysledk˚um. Volbou vektoru c, na jehoˇz odhadu z´avis´ı pravdˇepodobnost pˇr´ısluˇsnosti pixelu k faktoru, znaˇcnˇe penalizujeme slabˇs´ı pixely dan´eho faktoru. T´ım m˚uzˇ eme pˇrich´azet o cˇ a´ st d˚uleˇzit´e informace, kter´a se sice do znaˇcn´e m´ıry vykompenzuje v protilehl´em parenchymu, ale bylo by dobr´e zav´est pˇresnˇejˇs´ı odhadov´an´ı aprioren faktorov´ych obr´azk˚u, aby k potlaˇcov´an´ı slabˇs´ıch pixel˚u nedoch´azelo. Dalˇs´ı smˇer v´yzkumu by se mohl zaob´ırat ˇreˇsen´ım patologick´ych situac´ı. Integr´aln´ı model zat´ım nepoˇc´ıt´a s rozdˇelen´ım parenchym˚u na v´ıce faktor˚u, pˇrekryvem ledvin cˇ i dokonce jejich sr˚ustem. Protoˇze je lidsk´e tˇelo znaˇcnˇe variabiln´ı, pravdˇepodobnˇe nep˚ujde zahrnout vˇsechny probl´emy, kter´e mohou nastat, ale nˇekter´e typick´e ano.
61
7.3 Zhodnocen´ı Podaˇrilo se n´am vytvoˇrit integr´aln´ı model dle biologick´ych pˇredpoklad˚u, kter´e jsou splnˇeny pˇredevˇs´ım v prvn´ı f´azi scintigrafick´e studie ledvin. D´ıky tomu jsme schopni urˇcit automaticky diagnostick´y koeficient, relativn´ı ren´aln´ı clearanci. Aby model fungoval na cel´e sekvenci, bude muset b´yt jeˇstˇe upraven. Ve srovn´an´ı se standardn´ım modelem vˇsak dost´av´ame v´yraznˇe lepˇs´ı v´ysledky, a proto je nov´y integr´aln´ı model velk´ym posunem v anal´yze scintigrafick´ych studi´ı ledvin.
62
A Matematicky´ dodatek ´ A.1 Stopa matice, operator vec(), Kroneckeruv ˚ a ˇ Hadamarduv ˚ soucin Definice 5 Stopou cˇ tvercov´e matice A ∈ Rn×n rozum´ıme cˇ´ıslo tr(A) :=
n X
ai,i .
(A.1)
i=1
Vˇeta 6 Necht’ A ∈ Rn×n a tr(A) je jej´ı stopa. Pak plat´ı tr(A) = tr(A′ ).
(A.2)
Vˇeta 7 Necht’ A,B,C,D jsou matice takov´e, zˇe v´yraz ABCD d´av´a smysl a v´ysledkem tohoto n´asoben´ı je cˇ tvercov´a matice. Pak plat´ı, zˇe tr(ABCD) = tr(BCDA). Definice 8 Necht’ se matice A skl´ad´a ze sloupc˚u a:,1 , a:,2, ..., a:,n . Pak a:,1 a:,2 vec(A) = .. . . a:,n
(A.3)
(A.4)
Definice 9 Necht’ matice A ∈ Rn×m a matice B ∈ Rp×r . Pak Kroneckerov´ym souˇcinem matic A a B rozum´ıme v´yraz a1,1 B . . . a1,m B .. . .. A ⊗ B = ... (A.5) . . an,1 B . . . an,m B Vˇeta 10 Pro libovoln´e kompatibiln´ı matice (ve smyslu uveden´ych vzorc˚u) plat´ı vec(ABC) = (C ′ ⊗ A)vec(B) vec(ABC) = (A ⊗ C ′ )vec(B ′ ) vec(ABC) = (I ⊗ AB)vec(C) = (C ′ B ′ ⊗ I)vec(A) tr(A′ B) = vec(A)′ vec(B).
63
(A.6) (A.7) (A.8) (A.9)
Vˇeta 11 Pro Kronecker˚uv souˇcin plat´ı (A ⊗ B)′ = A′ ⊗ B ′ (A ⊗ B)−1 = A−1 ⊗ B −1 .
(A.10) (A.11)
Definice 12 Necht’ jsou d´any matice A, B stejn´eho rozmˇeru. Pak Hadamardov´ym souˇcinem matic A a B rozum´ıme matici C danou (C)i,j = (A)i,j · (B)i,j .
(A.12)
Znaˇc´ıme C = A ◦ B.
´ ı rozdelen´ ˇ A.2 Normaln´ ı ˇ e´ normaln´ ´ ı rozdelen´ ˇ A.2.1 V´ıcerozmern ı Necht’ je d´an vektor x ∈ Rn×1 . V´ıcerozmˇern´e norm´aln´ı rozdˇelen´ı vektoru x je pak: 1 1 ′ −1 − (x − µx ) Σx (x − µx ) , Nx (µx , Σx ) = 1 exp n 2 (2π) 2 |Σ| 2
(A.13)
kde Σx ∈ Rn×n je kovarianˇcn´ı matice. Pro momenty pak plat´ı: x b = µx c′ = Σx + µx µ′x xx D´ale pro C ∈ Rn×n plat´ı:
(A.14) (A.15)
′ x = tr(Σ ) + µ′ µ . xc x x x
(A.16)
Cx ∼ N(Cµx , CΣx C ′ ).
(A.17)
´ ı rozdelen´ ˇ A.2.2 Maticove´ normaln´ ı Necht’ je d´ana matice X ∈ Rn×p . Maticov´e norm´aln´ı rozdˇelen´ı matice X je pak: np
p
n
NX (µX , Σn ⊗ Φp ) = (2π)− 2 |Σn |− 2 |Φp |− 2 × 1 −1 −1 ′ ′ , (A.18) × exp − tr Σn (X − µX )(Φp ) (X − µX ) 2
kde Σn ∈ Rn×n a Φp ∈ Rp×p jsou symetrick´e pozitivnˇe definitn´ı matice, tr je stopa matice a ⊗ je Kronecker˚uv souˇcin. Pro momenty pak plat´ı: b = µX X
(A.19)
′ X = tr(Σ )Φ + µ ′ µ [ X n p X X
(A.20)
[′ = tr(Φp )Σn + µX µX ′ XX
(A.21)
64
D´ale pro C ∈ Rn×n a D ∈ Rp×p plat´ı: CXD ∼ N(CµX D, CΣn C ′ ⊗ D ′ Φp D) E(X ′ DX) = µ′X DµX + tr(Σn D)Φp E(XCX ′ ) = µX Cµ′X + Σn tr(CΦp ).
(A.22) (A.23) (A.24)
´ normaln´ ´ ım rozdelen´ ˇ A.2.3 Vektorizace v maticovem ı Nyn´ı se kr´atce pod´ıv´ame na pouˇzit´ı oper´atoru vec() na matici X ∈ Rn×p . Mˇejme opˇet X ∼ NX (µX , Σn ⊗ Φp ). Pak pro vec(X) plat´ı: vec(X) ∼ N(vec(µX ), Φp ⊗ Σn ).
(A.25)
Blokov´a struktura kovarianˇcn´ı matice vektoru vec(X) pak odpov´ıd´a vektoru vec(X) n´asledovnˇe: φ1,1 Σn ... ... ... x:,1 ... x:,2 φ2,2 Σn . . . ... (A.26) ∼ N(vec(µ ), .. .. ). .. .. X . . . . . . . ... ... . . . φp,pΣn x:,p Poznamenejme, zˇ e tato blokov´a struktura miz´ı, pokud se n´am nepodaˇr´ı udrˇzet Kroneckerovsk´a struktura kovarianˇcn´ı matice.
´ ı rozdelen´ ˇ A.2.4 Oˇrezane´ normaln´ ı Oˇrezan´e norm´aln´ı rozdˇelen´ı budeme definovat pro skal´arn´ı n´ahodnou promˇennou x na intervalu a < x ≤ b n´asledovnˇe: √ 2 exp((x − µ)2 ) tNx (x|µ, σ, a, b) = √ χ(a;b] (x), (A.27) πσ(erf (β) − erf (α)) kde α =
a−µ √ , 2σ (
β =
b−µ √ , 2σ
χ(a,b] (x) je charakteristick´a funkce intervalu (a, b] definovan´a jako
Rt 1 x ∈ (a, b] 2 a erf(t) = √2π 0 e−u du. 0 x∈ / (a, b] Momenty oˇrezan´eho norm´aln´ıho rozdˇelen´ı jsou: √ √ 2[exp(−β 2 ) − exp(−α2 )] x b=µ− σ √ π(erf (β) − erf (α)) √ √ 2[b exp(−β 2 ) − a exp(−α2 )] b 2 √ x = σ + µb x− σ . π(erf (β) − erf (α))
χ(a,b] (x) =
65
(A.28) (A.29)
ˇ A.3 Gamma rozdelen´ ı Necht’ je d´an n´ahodn´y skal´ar x. Pak definujeme gamma rozdˇelen´ı skl´aru x jako: Gx (a, b) = pro x, a, b > 0 a kde Γ(x) =
R∞
1 1 a−1 −xb x e Γ(a) b−a
(A.30)
tx−1 exp(−t) dt pro x > 0.
0
Pˇr´ısluˇsn´e momenty gamma rozdˇelen´ı jsou: a b a b 2 x = 2. b x b=
(A.31) (A.32)
´ ı rozdelen´ ˇ A.4 Exponencialn´ ı
Necht’ je d´an n´ahodn´y skal´ar x. Pak definujeme exponenci´aln´ı rozdˇelen´ı skl´aru x jako: Expx (λ) = λ exp(−λx)
(A.33)
pro x > 0. Pˇr´ısluˇsn´e momenty exponenci´aln´ıho rozdˇelen´ı jsou: 1 λ 1 xb2 = 2 . λ x b=
(A.34) (A.35)
´ ı rozdelen´ ˇ A.4.1 Oˇrezane´ exponencialn´ ı Oˇrezan´e exponenci´aln´ı rozdˇelen´ı budeme definovat pro skal´arn´ı n´ahodnou promˇennou x na intervalu a < x ≤ b n´asledovnˇe: t Expx (λ, a, b) =
1 exp(−λx)χ(a,b] (x), exp(λb) − exp(λa)
(A.36)
kde χ(a,b] (x) je charakteristick´a funkce intervalu, definovan´a v kapitole A.2.4. Prvn´ı moment oˇrezan´eho norm´aln´ıho rozdˇelen´ı je x b=
exp(λb)(1 − λb) − exp(λa)(1 − λa) . λ(exp(λa) − exp(λb))
66
(A.37)
ˇ e´ rozdelen´ ˇ A.5 Rovnomern ı Rovnomˇern´e rozdˇelen´ı pˇriˇrazuje vˇsem hodnot´am n´ahodn´e veliˇciny x na intervalu (a, b) stejnou pravdˇepodobnost. Jeho hustota pravdˇepodobnosti m´a n´asleduj´ıc´ı tvar: ( 1 x ∈ (a, b) (A.38) f (x) = b−a 0 x∈ / (a, b) Prvn´ı moment tohoto rozdˇelen´ı spoˇc´ıt´ame jako x b=
a+b . 2
67
(A.39)
ˇ logaritmu sdruˇzeneho ´ B Vypo ´ cet ˇ rozdelen´ ı Pomocn´y v´ypoˇcet logaritm˚u pravdˇepodobnostn´ıch rozdˇelen´ı z kapitoly 5: f(D): f (D|A, X, ω) = tN(AX ′ , ω −1Ip ⊗ In ) = p n 1 − pn −1 − − ′ ′ ′ = (2π) 2 |ω Ip | 2 |In | 2 exp − tr [ω(D − AX )(D − AX ) ] 2 pn 1 lnf (D|A, X, ω) ∝ ln ω − ω tr [DD ′ − 2AX ′ D ′ + AX ′ XA′ ] 2 2 f(ω):
(B.1) (B.2)
f (ω) = G(ϑ0 , ρ0 ) ∝ ω ϑ0 −1 e−ωρ0 ln f (ω) ∝ (ϑ0 − 1) ln ω − ωρ0
(B.3) (B.4)
f (A|Υ) = tN(ZA , Ip ⊗ Υ−1 ) r 1 pX ln vi − tr((A − ZA )Υ(A − ZA )′ ) ln f (A|Υ) ∝ 2 i=1 2
(B.5)
f(A):
f(c):
1 f (c) = tN(0r,1 , ς0 Ir ) ∝ exp(− tr(c′ ς0−1 Ir c)) 2 1 −1 ln f (c) ∝ − ς0 tr(c′ c) 2
(B.6)
(B.7) (B.8)
f(υ): f (υ|α0, β0 ) =
r Y i=1
ln f (υ) ∝ f(g):
Gυi (αi,0 , βi,0 ) ∝
r X i=1
r Y
i=1 r X
(αi,0 − 1) ln υi −
α
−1
(B.9)
υi βi,0
(B.10)
υi i,0 e−υi βi,0
i=1
f (g|ψ) = tN(0n , ψ −1 In ) n 1 ln f (g) ∝ ln ψ − ψ tr(g ′ g) 2 2
68
(B.11) (B.12)
f(ψ): f (ψ) = Gψ (ζ0 , η0 ) ∝ ψ ζ0 −1 e−ψη0 ln f (ψ) ∝ (ζ0 − 1) ln ψ − ψη0
(B.13) (B.14)
f(W): f (W ) = tN(MW , In ⊗ Ξ−1 W ) ∝ n 1 ′ −1 ′ ∝ |Ξ| 2 exp − tr (In ) (W − MW )ΞW (W − MW ) 2 r nX 1 ln f (W ) ∝ ln ξi − tr (W − MW )ΞW (W − MW )′ 2 i=1 2
(B.15) (B.16)
f(ξ): f (ξ|κ0, ν0 ) =
r Y i=1
ln f (ξ) ∝
Gξi (κi,0 , νi,0 ) ∝
r X i=1
r Y
i=1 r X
(κi,0 − 1) ln ξi −
κ
−1
(B.17)
ξi νi,0
(B.18)
ξi i,0 e−ξi νi,0
i=1
f(h): 1 f (h) = tN(0r,1 , τ0 Ir ) ∝ exp(− tr(h′ τ0−1 Ir h)) 2 1 −1 ′ ln f (h) ∝ − τ0 tr(h h) 2
(B.19) (B.20)
f(si ): f (si ) = U(0, n) = ln f (si ) ∝ ln
1 n
1 n−0
(B.21) (B.22)
f(li ): f (li ) = U(0, n − si ) = ln f (li ) ∝ ln
1 n − si
69
1 n − si − 0
(B.23) (B.24)
C Obrazove´ sekvence
10
10
10
10
10
20
20
20
20
20
30
30
30
30
30
40
40
40
40
40
50
50
50
50
50
60
60
60
60
20
40
60
20
40
60
20
40
60
60 20
40
60
10
10
10
10
10
20
20
20
20
20
30
30
30
30
30
40
40
40
40
40
50
50
50
50
50
60
60
60
60
20
40
60
20
40
60
20
40
60
40
60
10
10
10
10
20
20
20
20
20
30
30
30
30
30
40
40
40
40
40
50
50
50
50
50
60
60
60
60
40
60
20
40
60
20
40
60
40
60
10
10
10
10
20
20
20
20
20
30
30
30
30
30
40
40
40
40
40
50
50
50
50
50
60
60
60
60
40
60
20
40
60
20
40
60
20
40
60
20
40
60
20
40
60
60 20
40
Obr´azek C.1: studie 1: vyhodnocovan´a sekvence
70
60
60 20
10
20
40
60 20
10
20
20
60
10
10
10
10
20
20
20
20
30
30
30
30
40
40
40
40
50
50
50
50
60
60
60
20
40
60
20
40
60
60 20
40
60
10
10
10
10
20
20
20
20
30
30
30
30
40
40
40
40
50
50
50
50
60
60
60
20
40
60
20
40
60
40
60
10
10
10
20
20
20
20
30
30
30
30
40
40
40
40
50
50
50
50
60
60
60
40
60
20
40
60
40
60
10
10
10
20
20
20
20
30
30
30
30
40
40
40
40
50
50
50
50
60
60
60
40
60
20
40
60
20
40
60
20
40
60
20
40
60
60 20
40
60
Obr´azek C.2: studie 2: vyhodnocovan´a sekvence
71
60
60 20
10
20
40
60 20
10
20
20
Literatura [1] M. Caglar, G. Gedik, and E. Karabulut, “Differential renal function estimation by dynamic renal scintigraphy: influence of background definition and radiopharmaceutical,” Nuclear Medicine Communications, vol. 29, no. 11, p. 1002, 2008. ˇ ıdl and A. Quinn, The Variational Bayes Method in Signal Processing. Springer, [2] V. Sm´ 2006. [3] O. Tich´y, “Konvoluˇcn´ı parametrizace v anal´yze scintigrafick´ych obrazov´ych sekvenc´ı,” ˇ 2009. V´yzkumn´y u´ kol, FJFI CVUT. [4] A. Drastich, “Zobrazovac´ı syst´emy v l´ekaˇrstv´ı,” Brno: Ediˇcn´ı stˇredisko VUT Brno, 1990. [5] V. Ullman, “Nuklearn´ı medic´ına: Radioisotopov´a scintigrafie.” v´yukov´a str´anka, http://www.sweb.cz/AstroNuklFyzika/Scintigrafie.htm. ´ [6] K. Kupka, “Zobrazovac´ı metody nukle´arn´ı medic´ıny.” Uvodn´ ı text o zobrazovac´ıch metod´ach, http://unm.lf1.cuni.cz/. ˇ amal, “Uvod ´ [7] M. S´ do studia nukle´arn´ı medic´ıny.” v´yukov´e slidy na pˇredmˇet Nuklearn´ı medic´ına, http://unm.lf1.cuni.cz/. [8] H. Francois, Dynamic renal imaging in obstructive renal pathology. European Association of Nuclear Medicine, 2009. [9] J. Bassingthwaighte, G. Holloway, et al., “Estimation of blood flow with radioactive tracers,” in Seminars in Nuclear Medicine, vol. 6, pp. 141–161, Elsevier, 1976. [10] E. Durand, M. Blaufox, K. Britton, O. Carlsen, P. Cosgriff, E. Fine, J. Fleming, C. Nimmon, A. Piepsz, A. Prigent, et al., “International Scientific Committee of Radionuclides in Nephrourology (ISCORN) consensus on renal transit time measurements,” in Seminars in nuclear medicine, vol. 38, pp. 82–102, Elsevier, 2008. [11] A. Kuruc, J. Caldicott, and S. Treves, “Improved Deconvolution Technique for the Calculation of Renal Retention Functions.,” COMP. AND BIOMED. RES., vol. 15, no. 1, pp. 46–56, 1982. [12] H. Attias and C. Schreiner, “Blind source separation and deconvolution: the dynamic component analysis algorithm,” Neural Computation, vol. 10, no. 6, pp. 1373–1424, 1998.
72
[13] J. Fleming and P. Kemp, “A comparison of deconvolution and the Patlak-Rutland plot in renography analysis,” Journal of Nuclear Medicine, vol. 40, no. 9, p. 1503, 1999. [14] C. Beckmann and S. Smith, “Probabilistic independent component analysis for functional magnetic resonance imaging,” IEEE transactions on medical imaging, vol. 23, no. 2, pp. 137–152, 2004. [15] R. Lawson, “Application of mathematical methods in dynamic nuclear medicine studies,” Physics in medicine and biology, vol. 44, pp. R57–R98, 1999. ˇ ıdl and M. S´ ˇ amal, “Robust detection of linear part of patlak-rutland plots.” UTIA ´ [16] V. Sm´ ˇ v.v.i, Research Report 2243, 2008. AV CR, [17] P. Veltri, A. Vecchio, and V. Carbone, “Proper orthogonal decomposition analysis of spatio-temporal behavior of renal scintigraphies,” Physica Medica, 2009. [18] A. Peters, “The kinetic basis of glomerular filtration rate measurement and new concepts of indexation to body size,” European Journal of Nuclear Medicine and Molecular Imaging, vol. 31, no. 1, pp. 137–149, 2004. [19] M. Huˇskov´a, “Bayesovsk´e met´ody,” Praha, SPN, 1985. [20] M. Beal, Variational Algorithms for Approximate Bayesian Inference. PhD thesis, University College London, 2003. [21] A. Dempster, N. Laird, D. Rubin, et al., “Maximum likelihood from incomplete data via the EM algorithm,” Journal of the Royal Statistical Society, vol. 39, no. 1, pp. 1–38, 1977. ˇ ıdl, The Variational Bayes Approach in Signal Processing. PhD thesis, University [22] V. Sm´ of Dublin, Trinity College, 2004. [23] R. Neal and G. Hinton, “A view of the EM algorithm that justifies incremental, sparse, and other variants,” Learning in Graphical Models, vol. 89, pp. 355–368, 1998. [24] R. Kass and A. Raftery, “Bayes Factors.,” Journal of the American Statistical Association, vol. 90, no. 430, 1995. [25] T. Minka, A family of algorithms for approximate Bayesian inference. PhD thesis, Massachusetts Institute of Technology, 2001. [26] T. Minka, “Old and new matrix algebra useful for statistics. notes,” 2000. [27] C. Andrieu, N. De Freitas, A. Doucet, and M. Jordan, “An Introduction to MCMC for Machine Learning,” Machine Learning, vol. 50, pp. 5–43, 2003. [28] S. Kullback and R. Leibler, “On information and sufficiency,” Annals of Mathematical Statistics, vol. 22, no. 1, pp. 79–86, 1951.
73
[29] Z. Ghahramani and M. Beal, “Variational Inference for Bayesian Mixtures of Factor Analysers,” Advances in Neural Information Processing Systems 12: Proceedings of the 1999 Conference, 2000. [30] M. Sato, “Online Model Selection Based on the Variational Bayes,” 2001. ˇ amal, M. K´arn´y, M. Surov´a, E. Maˇr´ıkov´a, and Z. Dienstbier, “Rotation to simple [31] M. S´ structure in factor analysis of dynamic radionuclide studies.,” Physics in Medicine & Biology, vol. 32, no. 3, pp. 371–382, 1987. [32] O. Tich´y, “Anal´yza scintigrafick´ych obrazov´ych sekvenc´ı v l´ekaˇrsk´e diagnostice,” 2008. ˇ Bakal´aˇrsk´a pr´ace, FJFI CVUT.
74
