Tomografia komputerowa a układy równań. Metoda Kaczmarza.

Równania liniowe oraz układy równań liniowych są jedną z tych rzeczy, których musimy się uczyć w szkole. Nieraz się pewnie zastanawiając po co to komu potrzebne? To, że umiemy i cały czas udoskonalamy metody rozwiązywania układów równań liniowych legło u podstaw wynalezienia tomografu komputerowego. Matoda Kaczmarza, bo o niej mowa, jest ponadto świetnym przykładem tego, że warto zajmować się (matematycznymi) rzeczami z pozoru nie mającymi zastosowania w prawdziwym życiu.

Równania liniowe są najprostszym rodzajem równań, które uczymy się rozwiązywać na lekcjach matematyki. Nieco trudniejsze do rozwiązania są układy równań liniowych. Te spotykane w szkole zazwyczaj składają się z dwu równań z dwiema niewiadomymi. Uczymy się, że jeżeli Zosia za dwie (takie same) bułki oraz batona zapłaciła 5 zł, a Jola za taką samą bułkę i takiego samego batona zapłaciła 3 zł, to taki problem możemy rozwiązać przy pomocy układu równań
\[\left\{\begin{array}{lllll}
2x & + & y & = & 5\\ x & + & y & = & 3.
\end{array}\right.\]

Znacznie rzadziej spotykamy się w szkole z układami równań, mającymi więcej niż dwie niewiadome. Z takimi studenci matematyki oraz np. kierunków technicznych zaprzyjaźniają się na pierwszym roku studiów. Jednakże bardzo rzadko wspomina się o zastosowaniach, skupiając się głównie na umiejętności ich rozwiązywania.

Wraz z rozwojem nauki zmuszeni jesteśmy do rozwiązywania nowych i często coraz bardziej złożonych obliczeniowo problemów. Niektóre z nich, jak np. modele pogody prowadzą właśnie do układów równań liniowych. Dlatego warto rozwijać znane i wymyślać nowe metody np. rozwiązywania układów równań liniowych.

W tym wpisie opiszemy jeden z przykładów, gdzie właśnie one znalazły zastosowanie (przez krótki czas, później używano lepszych metod) oraz użytą wówczas metodę ich (w tym przypadku przybliżonego) rozwiązywania. Dokładniej, zajmiemy się bardzo ważnym wynalazkiem jakim jest tomograf komputerowy, za konstrukcję którego Sir Godfrey Hounsfield oraz Allan Cormack otrzymali w 1979 Nagrodę Nobla w dziedzinie medycyny.

Tomografia, to sposób zobaczenia co jest wewnątrz obiektu bez otwierania go.
A. Markoe ,,Analytic Tomography”

Mimo iż tomograf został skonstruowany na przełomie lat sześćdziesiątych i siedemdziesiątych ubiegłego wieku, to jego historia zaczęła się tak naprawdę znacznie wcześniej.

W 1895 roku Wilhelm Röntgen odkrył promieniowanie znane dzisiaj (w Polsce i Niemczech) jako promieniowanie rentgenowskie lub (na świecie) jako promieniowanie X. Za swoje odkrycie otrzymał w 1901 Nagrodę Nobla. Poniżej jedno z pierwszych zdjęć rentgenowskich przedstawiające dłoń Alberta von Köllikera. Wykonał je sam Röntgen w 1896 roku w Instytucie Fizyki Uniwersytetu w Würzburgu.

Zwykłe prześwietlenie pozwala zobaczyć bardzo dobrze struktury (np. kości), które absorbują sporo promieniowania. Gorzej np. z narządami wewnętrznymi, które absorbują go znaczniej mniej i jeżeli w linii między źródłem promieniowania a detektorem znajduje się kość, to przysłoni inne narządy na rentgenogramie.

A co jeżeli zrobimy więcej prześwietleń? Taka, w skrócie, jest idea działania tomografu komputerowego. Urządzenie wysyła promieniowanie rentgenowskie pod różnymi kątami, które przeszedłszy przez badany obiekt (np. pacjenta) trafia do detektora. Jesteśmy w stanie zmierzyć natężenie promieniowania wejściowego oraz wyjściowego i na tej podstawie wywnioskować jak wygląda wewnętrzna struktura badanego obiektu.

To, że coś takiego jest możliwe wiemy z pracy z 1917 roku austriackiego matematyka Johana Radona, który udowodnił, że gdy na płaszczyźnie jest określona funkcja spełniająca pewne warunki, to jesteśmy w stanie ją odtworzyć mając informacje o wartościach całek z tej funkcji wzdłuż prostych na płaszczyźnie. Wprowadził w tym celu coś co nazywamy dzisiaj transformatą Radona. Wyniki Radona zostały ponownie odkryte przez Cormacka w jego pracach z 1963 oraz 1964. Przez długi czas ,,matematyczna” część tomografów komputerowych opierała się właśnie na wykorzystaniu tego co odkryli Radon oraz Cormack. Z pewnością i o tym wpis się na tym blogu pojawi.

Jednakże pierwsze tomografy działały nieco inaczej. Rozwiązywały układ równań liniowych wykorzystując metodę przybliżonego ich rozwiązywania wymyśloną (i opublikowaną w 1937) przez polskiego matematyka Stefana Kaczmarza i ponownie odkrytą w pracy Gordona, Bendera oraz Hermana. Zanim jednak omówimy metodę Kaczmarza, to najpierw opiszemy w jaki sposób układy równań liniowych mogą pojawiać się w kontekście tomografii komputerowej. Idea jest prosta.

Gdy promieniowanie przechodzi przez pewien jednorodny ośrodek, pokonując w nim drogę \(\alpha\), to zgodnie z prawem Lamberta-Beera, jego natężenie zmienia się wg wzoru
\[I=I_0e^{-k\alpha},\]
gdzie \(I_0\) – natężenie wejściowe, \(I\) – natężenie wyjściowe, \(k\) – współczynnik absorpcji (zależny od ośrodka).
W przypadku, gdy promieniowanie przechodzi przez kilka ośrodków, otrzymujemy
\[I=I_0e^{-k_1\alpha_1}e^{-k_2\alpha_2}\cdot\ldots\cdot e^{-k_n\alpha_n}=I_0e^{-(k_1\alpha_1+\ldots+k_n\alpha_n)}.\]
Dzieląc obustronnie przez \(I_0\), a następnie logarytmując otrzymujemy równanie liniowe
\[-\ln{\frac{I}{I_0}}=k_1\alpha_1+\ldots+k_n\alpha_n.\]
Znając \(I\) oraz \(I_0\) jesteśmy w stanie obliczyć wartość \(\ln{\frac{I}{I_0}}\).

Innymi słowy, każdy promień przechodzący przez np. obiekt taki jak poniżej prowadzi do pewnego równania liniowego. W tym przypadku z trzema niewiadomymi.

Znając punkt oraz kąt wyemitowania jesteśmy w stanie stwierdzić przez które kwadraciki promień przeszedł oraz jaką pokonał w nich drogę. Jedynymi niewiadomymi są współczynniki absorpcji.

Teraz pokażemy w jaki sposób można to wykorzystać w tomografii. Dla prostoty posłużymy się modelem dwuwymiarowym.

Pierwszym krokiem jest podzielenie obiektu naszych badań na mniejsze części, tj. piksele. Są one małymi kwadratami. W przypadku trójwymiarowym mielibyśmy sześciany zwane wokselami.

Traktujemy każdy piksel tak jakby składał się z jednorodnego materiału. Jest to przybliżenie, ale dzięki niemu jesteśmy w stanie otrzymać przybliżony obraz wnętrza obiektu.

Wiemy już, że z każdym promieniem przechodzącym przez obiekt jesteśmy w stanie stowarzyszyć równanie liniowe. Naświetlanie go promieniami pod różnymi kątami prowadzi do układu równań.

Niewiadomymi są współczynniki absorpcji wszystkich pikseli. Drogę pokonaną przez dany promień w każdym pikselu jesteśmy w stanie wyliczyć znając kąt oraz punkt wyemitowania. Im lepszą chcemy otrzymać dokładność, tym na więcej pikseli powinniśmy podzielić nasz obiekt.

Otrzymanie zadowalającej rozdzielczości wiąże się z rozwiązaniem sporego układu równań. Pierwsze tomografy pozwalały obrazować przekroje mózgu przy rozdzielczości 80\(\times\)80 pikseli. Tak otrzymany układ równań ma 1600 niewiadomych. Pocieszające może być, że znaczna część promieni przechodzi jedynie przez znikomą część pikseli. Dzięki temu tak otrzymany układ równań ma wiele współczynników równych 0.

Oczywiście to tylko matematyczna idea, od której często wszystko się zaczyna i w której wiele obiektów jest wyidealizowanych. W przypadku zastosowań często należy się zmierzyć jeszcze z napotykanymi problemami praktycznymi.

Tylko jak taki układ rozwiązać? Z matematycznego punktu widzenia o układach równań liniowych wiemy wszystko. Mając jakikolwiek układ równań liniowych znamy metody pozwalające stwierdzić czy taki układ ma w ogóle jakiekolwiek rozwiązanie. A jeżeli ma, to w pełni opisać zbiór jego rozwiązań. Podstawowe metody, których uczymy się w szkole oraz na studiach, mają jedną wadę. Kompletnie nie nadają się do rozwiązywania dużych układów równań nawet przy użyciu komputerów. Duże układy pojawiają się w wielu zagadnieniach praktycznych, oprócz tomografii, np. w modelach prognozy pogody. Dlatego tak ważne jest abyśmy, jako cała ludzkość, potrafili ulepszać znane i wymyślać nowe, efektywniejsze metody ich rozwiązywania. Zanim opiszemy na czym polega metoda Kaczmarza, musimy sobie przypomnieć kilka informacji o układach równań liniowych.

Ze szkoły wiemy, że układ dwu równań z dwiema niewiadomymi opisuje dwie proste w układzie współrzędnych. Lub jedną, jeśli oba równania opisują tę samą prostą. Gdy takie proste się przecinają, to współrzędne punktu przecięcia są rozwiązaniem. Może się również zdarzyć, że proste nie mają punktów wspólnych. Wówczas układ jest sprzeczny. Trzecią możliwością jest ta gdy równania opisują tę samą prostą. Taki układ ma nieskończenie wiele rozwiązań.

Ciekawiej to wygląda, gdy mamy układ równań z trzema niewiadomymi. Każde takie równanie opisuje pewną płaszczyznę zanurzoną w przestrzeni trójwymiarowej. Jeżeli układ ma jedno rozwiązanie, to płaszczyzny opisane jego równaniami przecinają się w dokładnie jednym punkcie.

Płaszczyzny mogą się również przecinać wzdłuż prostej.

Podobnie mogą w ogóle nie mieć punktów wspólnych jak i wszystkie równania mogą opisywać tę samą płaszczyznę. Zatem, z geometrycznego punktu widzenia, zbiór rozwiązań układu równań liniowych z trzema niewiadomymi może być zbiorem pustym, punktem, prostą bądź płaszczyzną.

Gdy równanie liniowe ma \(n\) niewiadomych, to opisuje \((n-1)\)-wymiarową przestrzeń zanurzoną w \(n\)-wymiarowej (matematycy by powiedzieli \((n-1)\)-wymiarową hiperpłaszczyznę).

Metoda Kaczmarza jest algorytmem, który z układem równań liniowych stowarzysza ciąg jego \((K_n)\) przybliżonych rozwiązań. Jeżeli układ jest niesprzeczny, to ciąg ten jest zbieżny do jednego z rozwiązań.

Idea jest bardzo prosta. Na początek wybieramy dowolny punkt \(K_1\) przestrzeni, w której są zanurzone hiperpłaszczyzny opisane równaniami układu. Następnie rzutujemy ten punkt prostopadle na hiperpłaszczyznę \(\Pi_1\) opisaną pierwszym równaniem. Otrzymany punkt rzutujemy prostopadle na hiperpłaszczyznę \(\Pi_2\) opisaną drugim równaniem itd. Gdy dojdziemy do hiperpłaszczyzny opisanej ostatnim równaniem, to rzutujemy znowu na tę opisaną pierwszym.

W przypadku układu dwu równań z dwiema niewiadomymi ową zbieżność ciągu \((K_n)\) łatwo zobrazować. Aczkolwiek rozwiązywanie takiego układu algorytmem Kaczmarza jest niezbyt optymalne.

Metoda Kaczmarza swą praktyczną użyteczność zawdzięcza temu, że mając \(K_n\) jesteśmy w stanie łatwo wyznaczyć \(K_{n+1}\) korzystając ze wzoru \[K_{n+1}=K_n-\dfrac{\textbf{a}_i\circ K_n-b_i}{||\textbf{a}_i||^2}\cdot\textbf{a}_i,\]
gdzie \(i=k\mod m\) (oczywiście \(m\) jest liczbą równań). Wyprowadzenie tego wzoru oraz dowód poprawności algorytmu Kaczmarza, na końcu.

Jak już wspomnieliśmy, algorytm Kaczmarza został szybko zastąpiony przez wydajniejsze metody otrzymywania tomogramów oparte na wynikach Radona oraz Cormacka (tu nie chodziło o układy równań). Jednakże dzięki swej prostocie i ulepszeniom nadal znajduje zastosowanie.

Jest to również przykład tego, że warto się zajmować rzeczami, które z pozoru wydają się niezbyt potrzebne. Zarówno Kaczmarz jak i Radon w swych pracach nie podawali żadnych zastosowań. Swoją drogą praca Kaczmarza przez długi czas nie wzbudzała zbyt wielkiego zainteresowania matematyków. Pierwsze jej cytowanie pojawiło się najprawdopodobniej dopiero w 1948 roku. Choć spory wpływ na to mogła mieć II wojna światowa, której zresztą Kaczmarz nie przeżył. Okoliczności jego śmierci nie są znane, ale wiele wskazuje na to, że po napaści Niemiec na Polskę zginął już w pierwszych dniach wojny podczas bombardowania Niska. Z czasem praca Kaczmarza została doceniona i doczekała się już ponad tysiąca cytowań.

Wartą odnotowania ciekawostką o Kaczmarzu jest to, że był również nauczycielem w szkole średniej. Mimo posiadania stopnia dr. habilitowanego musiał zdać egzamin aby móc nauczać w szkole. W komisji egzaminacyjnej zasiadł m.in. sam Stefan Banach. Więcej szczegółów (m.in. ciekawe pytania egzaminacyjne) można znaleźć tutaj.

Zanim przejdziemy do dowodu metody Kaczmarza spójrzmy jak wyglądały jedne z pierwszych tomogramów.

Rezultat może nie jest oszałamiający, lecz mówimy o początku lat siedemdziesiątych ubiegłego wieku. Wówczas był to ogromny postęp.

Metoda Kaczmarza – dowód

Udowodnimy teraz, że ciąg punktów \((K_n)\) otrzymany przy użyciu metody Kaczmarza rzeczywiście dąży do rozwiązania, o ile wyjściowy układ równań jest niesprzeczny.

Niech układ równań liniowych (mający \(N\) niewiadomych)\[\textbf{a}_i\circ\textbf{x}=b_i,\]
gdzie \(i\in\{1,\ldots,m\}\), będzie niesprzeczny, a \(G\in\mathbb{R}^N\) jego dowolnym rozwiązaniem. Przez \(\Pi_i\) oznaczmy hiperpłaszczyznę w \(\mathbb{R}^N\) opisaną równaniem \[\textbf{a}_i\circ\textbf{x}=b_i\text{ dla }i\in\{1,\ldots,m\}.\]

Z kontrukcji ciągu \((K_n)\) oraz twierdzenia Pitagorasa otrzymujemy \[|GK_{n+1}|^2+|K_nK_{n+1}|^2=|GK_n|^2.\]

Stąd \[|GK_n|\geqslant |GK_{n+1}|,\]
a więc ciąg \((|GK_n|)\) jest nieujemny oraz nierosnący czyli ograniczony. Jako ciąg monotoniczny oraz ograniczony, jest zbieżny. Z tego wynika, że \((K_n)\) jest ograniczony. Na mocy twierdzenia Bolzano-Weierstrassa zawiera podciąg zbieżny \((K_{n_j})\). Oznaczmy przez \(G_0\) jego granicę.

Ponieważ wyjściowy układ równań opisuje skończenie wiele hiperpłaszczyzn w \(\mathbb{R}^N\), to istnieje takie \(i\in\{1,\ldots,m\}\), że nieskończenie wiele wyrazów ciągu \((K_{n_j})\) leży na hiperpłaszczyźnie \(\Pi_i\), która jest zbiorem domkniętym w \(\mathbb{R}^N\). Zatem \(G_0\in\Pi_i\).

Z tego, że \(|GK_{n+1}|^2+|K_nK_{n+1}|^2=|GK_n|^2\) oraz ze zbieżności ciągu \((|GK_n|)\) mamy \[\lim\limits_{n\to\infty}|K_{n+1}K_n|=0.\] Stąd oraz z tego, że \((K_{n_j})\) jest zbieżny wynika, że \((K_{n_j+1})\) również jest zbieżny.

Ale \(K_{n_j+1}\in\Pi_{i+1}\) gdy \(m>i\) oraz \(K_{n_j+1}\in\Pi_1\) dla \(i=m\).

Zatem \(G_0\in\Pi_{i+1}\) lub \(G_0\in\Pi_1\) w zależności od \(i\). Dostajemy więc, że \(G_0\in\Pi_i\) dla każdego \(i\in\{1,\ldots,m\}\), czyli jest to rozwiązanie układu.

Przypomnijmy, że ciąg \((|GK_n|)\) jest zbieżny dla dowolnego rozwiązania \(G\) układu. W szczególności dla \(G=G_0\) mamy
\[|GK_{n_j}|\to 0,\]
a stąd \(K_n\to G_0\).

Rzut punktu na hiperpłaszczyznę – wyprowadzenie wzoru

Jak wspomnieliśmy wcześniej, praktyczna użyteczność metody Kaczmarza wynika tego, że łatwo opisać wzorem rzut punktu \(
\textbf{p}\in\mathbb R^n\) na hiperpłaszczyznę \(\Pi\subset\mathbb R^n\) opisaną równaniem \[\textbf{a}\circ\textbf{x}=b\] dla pewnych \(\textbf{a}\in\mathbb R^n\) oraz \(b\in\mathbb R\).

Wektor \(\textbf{a}\) jest prostopadły do płaszczyzny \(\Pi\). Ponieważ rzutujemy \(\textbf{p}\) prostopadle na \(\Pi\), to wektor \(\textbf{p}-\textbf{p}_0\), gdzie \(\textbf{p}_0\) jest szukanym rzutem, jest równoległy do \(\textbf{a}\). Lub równy \(\textbf{0}\), gdy \(\textbf{p}\in\Pi\). Czyli\[\textbf{p}-\textbf{p}_0=\lambda\textbf{a}\]
dla pewnego \(\lambda\in\mathbb R\). Ponieważ \[\textbf{p}_0=\textbf{p}-(\textbf{p}-\textbf{p}_0),\] to \[\textbf{p}_0=\textbf{p}-\lambda\textbf{a}.\] Ale \(\textbf{p}_0\in\Pi\), więc \[b=\textbf{a}\circ\textbf{p}_0=\textbf{a}\circ (\textbf{p}-\lambda\textbf{a})=\textbf{a}\circ\textbf{p}-\lambda||\textbf{a}||^2.\] Stąd otrzymujemy

\[\lambda=\dfrac{\textbf{a}\circ\textbf{p}-b}{||\textbf{a}||^2}.\]

Zatem \[\textbf{p}_0=\textbf{p}-\dfrac{\textbf{a}\circ\textbf{p}-b}{||\textbf{a}||^2}\cdot\textbf{a}.\]

Na koniec literatura oraz pozycje, które okazały się pomocne przy tworzeniu tego tekstu. Za nią natomiast treść oryginalnego artykułu Kaczmarza.

Literatura

S. Kaczmarz, Przybliżone rozwiązywanie układów równań liniowych – Angenäherte Auflösung von Systemen linearer Gleichungen, Bull. Int. Acad. Polon. Sci., Cl. Sci. Math., Ser. A, Sci. Nat. (1937), s. 22 – 24
A. Cegielski Bibliografia związana z pracą Kaczmarza
A. Cormack, Representation of a function by its line integrals, with some radiological applications, J. Appl. Phys. 34 (1963), s. 2722 – 2727
G. Hounsfield, Computerized transverse axial scanning (tomography). Part I. Description of system, Br J Radiol 46 (1973), s. 1016 – 1022
J. Ambrose, Computerized transverse axial scanning (tomography). Part 2. Clinical application, Br J Radiol. (46) 1973, s. 1023-1047
T.G. Feeman, The Mathematics of Medical Imaging, Springer (2015)
R. Gabor, R. Bender, G.T. Herman, Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography, J. Theor. Biol. 29 (1970), s. 471–481
A. Markoe, Analytic Tomography, Cambridge Univ. Press (2006)
J. Radon, Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten, Ber. Verh. Säche. Akad. Wiss. 69 (1917), s. 262–277