Nem statisztikai iteratív megoldások

Tankönyv Fizikusoknak » Az orvosi képalkotás matematikai alapjai » Képrekonstrukció algebrai összefüggésekkel » Nem statisztikai iteratív megoldások

Az előző szakaszban bevezetett jelölésekkel, a képrekonstrukció alapegyenlete felírtható a következőképpen:
$$\mathbf{y}=\mathbf{A}\mathbf{x}$
Ne felejtsük el, hogy az A mátrix nem négyzetes, ritka mátrix, és elemeinek száma akár 10¹².

A Moore-Penrose pszeudoinverz

Mivel az A a mátrix nem négyzetes, a legtöbb szokásos eljárás (pl. Gauss-Siedel iteráció) közvetlenül nem alkalmazható. Akár túl-, akár alulhatározott az egyenletünk, a keresett megoldás kritériumát felírthatjuk a következőképpen:
$$ \underset{\mathbf{x}}{min}\left \| \mathbf{y}-\mathbf{A}\boldsymbol{\mathbf{x}} \right \|$

ahol ||.|| választásunknak megfelelő norma. Tegyük fel, hogy $$ \widetilde{\mathbf{x}}$ minimalizálja a fenti kifejezést, továbbá válasszunk skalárszorzatot normának. Ekkor egy tetszőleges
$$ \mathbf{v}$ vektorra és egy t valós számra igaznak kell lennie, hogy
$$ \partial _{t}\left \| \left (\mathbf{y}-\mathbf{A}\left (\widetilde{\mathbf{x}} -t\mathbf{v} \right ) \right )\right \|\mid _{t=0} =0$
A deriváltat kiszámítva
$$ \mathbf{A}\mathbf{v} \left (\mathbf{y}-\mathbf{A}\left (\widetilde{\mathbf{x}} -t\mathbf{v} \right ) \right )\mid _{t=0} =0$
azaz
$$ \mathbf{v} \mathbf{A}^{T}\left (\mathbf{y}-\mathbf{A}\widetilde{\mathbf{x}} \right ) =0$
Ennek minden v-re teljesülnie kell, tehát:
$$ \mathbf{A}^{T}\mathbf{y}= \mathbf{A}^{T}\mathbf{A}\widetilde{\mathbf{x}}$

Az A^TA mátrix már négyzetes, és ha létezik inverze, akkor:
$$ \widetilde{\mathbf{x}}= \left (\mathbf{A}^{T}\mathbf{A} \right )^{-1}\mathbf{A}^{T}\mathbf{y}$
Az A⁺ Moore-Penrose pszeudoinverz:
$$ \mathbf{A}^{+}= \left (\mathbf{A}^{T}\mathbf{A} \right )^{-1}\mathbf{A}^{T}$

Az A⁺ kiszámításához használhatunk UV dekompozíciót, vagy más hagyományos eljárást, illetve kiszámítható a következő egyszerű (de nem annyira hatékony) iteratív képlettel is (Ben-Israel- and Cohen-iteráció):
$$\mathbf{A}_{(k)}^{+}=2\mathbf{A}_{(k-1)}^{+} +\mathbf{A}_{(k-1)}^{+}\mathbf{A}_{}\mathbf{A}_{(k-1)}^{+}$

A Kaczmarz-iteráció

Az PET, SPECT és CT rendszermátrixinak mérete a manapság könnyen elérhető numerikus számítási kapacitásokhoz képest nagy, illetve a hagyományos algoritmusok nem túl hatékonyak a ritkamátrixokon végzett műveletekkel. Egy egyszerű és gyors, továbbá stabil algoritmus a Kaczmarz-iteráció. Lényege az, hogy egy nxm dimenziós A mátrix tulajdonképpen n db m dimenziós hipersík egyenlete, melyek metszéspontja az x megoldás. Ha a rendszer túlhatározott és zajjal terhelt, nem feltétlenül lesz pontosan egy metszéspont, ha alulhatározott, elképzelhető, hogy valamilyen tértartományig szűkíthetjük le az ismeretlent. A Kaczmarz-iteráció lényege, hogy úgy próbálunk közelíteni a megoldás tartománya felé, hogy minden iterációs lépésben merőlegesen levetítjük az x pontot a mátrix egy soraként előálló hipersíkra.
Vegyünk példaként egy 2x2-es A mátrixot a_ii elemekkel. Ekkor a kövezkező két elemű egyenes-egyenletrendszert kapjuk:
$$ \left.\begin{matrix} y_{1}=a_{11}x_{1}+a_{12}x_{1} \\ y_{2}=a_{21}x_{2}+a_{22}x_{2} \end{matrix}\right\}$
melyek metszéspontja a megoldás. Az ábrán szemléltettük az iteráció első négy lépését.

A Kaczmarz-iteráció két dimenzióban négy iterációs lépésig

Ebben a kétdimenziós esetben jól látszik, hogy

ha az egyenesek egymásra merőlegesek volnának, azaz a mátrix sorai lineárisan függetlenek, akkor 2 lépésben megkapnánk az egzakt megoldást
ha az egyenesek párhuzamosak, azaz az egyenletrendszer ellentmondást tartalmaz, az iteráció vég nélkül oszcillál a két megoldás között
ha a két egyenes kis szöget zár be egymással, azaz a sorok kevéssé függetlenek, akkor az iteráció lassan konvergál, minél függetlenebbek, annál gyorsabban

A Kaczmarz-iteráció előnye, hogy a mátrixnak egyszerre mindössze egyetlen sorát használjuk fel, nincsen szükség a teljes mátrix kezelésére egyszerre.

Nézzük a vetítés lépéseit most már általános esetre, az i index jelölje, hogy a mátrix melyik sorát használjuk éppen, azaz hányadik sorvektorát. x^k legyen a k. lépésben kapott megoldás. A k+1. iterációs lépésben rá kell vetítenünk ezt a pontot az i. egyenlet által meghatározott hipersíkra. A hipersík egyenlete:
$$ \mathbf{A}^{i}\mathbf{x}=y_{i}$
Korábban láttuk, hogy egy ilyen egyenlet az $$ \mathbf{A}^{i}/\sqrt{\mathbf{A}^{i}\mathbf{A}^{i}}$ irányú egységvektorra merőleges síkot írja le, melynek távolsága az origótól $$ y_{i}/\sqrt{\mathbf{A}^{i}\mathbf{A}^{i}}$ . Ha tehát a k. lépésben kapott eredményünket rá akarjuk vetíteni az i. hipersíkra, $$ \mathbf{A}^{i}/\sqrt{\mathbf{A}^{i}\mathbf{A}^{i}}$ irányban kell haladnunk, konvenció szerint $$ -\alpha$ távolságot. Az egyszerűbb képletek érdekében sűrítsük bele ebbe a távolságba az egységvektor normálizációját is, azaz felírhatjuk a következő pont helyét:
$$ \mathbf{x}^{k+1}=\mathbf{x}^{k}-\alpha \mathbf{A}^{i}$

Az új pontnak rajta kell lennie a az i. hipersíkon, tehát:
$$ \mathbf{A}^{i}\left (\mathbf{x}^{k}-\alpha \mathbf{A}^{i} \right )=y_{i}$
Ebből kifejezhető $$ \alpha$
$$ \alpha=\frac{\mathbf{A}^{i}\mathbf{x}^{k}-y_{i}}{\mathbf{A}^{i}\mathbf{A}^{i}}$
A k+1. iterációs lépésben kapott megoldás tehát:
$$ \mathbf{x}^{k+1}=\mathbf{x}^{k}-\frac{\mathbf{A}^{i}\mathbf{x}^{k}-y_{i}}{\mathbf{A}^{i}\mathbf{A}^{i}} \mathbf{A}^{i}$

A Kaczmarz-iteráció az orvosdiagnosztikai irodalomban ART (Algebraic Reconstruction Technique) néven is megtalálható. Az eljárás egyszerű, a korszerű rekonstrukciók között már csak továbbfejlesztett változataival találkozhatunk. Érdekesség, hogy tekinthető a Kaczmarz-iteráció a később tárgyalt OSEM iterációs séma egy extrém kis részhalmazú esetének.

Következő fejezetünk a mátrixegyenlet olyan bővítésével foglalkozik, melyben az egyenlet elemeit valamilyen valószínűségi változó várható értékének tekintjük.

Export

Medical Imaging

Site Language: English

Nem statisztikai iteratív megoldások

A Moore-Penrose pszeudoinverz

A Kaczmarz-iteráció

Sidebar

Medical Imaging

Site Language: English