Tablica NumPy nie ma położenia fizycznego

Segmentacja znajduje się kilka centymetrów od anatomii, dla której ją narysowano. Przyczyna rzadko leży tam, gdzie szuka się jej najpierw: wartości maski są poprawne, rysunek wykonano prawidłowo, a kod zakończył działanie bez błędu. Awaria nastąpiła wcześniej — maskę potraktowano tak, jakby miała położenie w ciele pacjenta, choć zawierała jedynie położenie w tablicy. Ta sekcja dotyczy jednego faktu, który zapobiega wielu błędom tej klasy: tablica NumPy jest pojemnikiem na liczby bez położenia w ciele pacjenta.

Ten sam fakt rozstrzyga mniej widowiskowe kwestie. Wczytujesz obraz CT do NumPy i otrzymujesz tablicę o kształcie (100, 512, 512). Inna osoba wczytuje inne badanie i również otrzymuje (100, 512, 512). Czy obie tablice są porównywalne? Czy można odjąć jedną od drugiej? Nie da się tego stwierdzić na podstawie samych tablic.

Ten sam kształt, inny świat

Dwie tablice o identycznym kształcie mogą przedstawiać:

• tego samego pacjenta i tę samą serię, ale różne rekonstrukcje (inne odstępy);
• tego samego pacjenta zeskanowanego dwukrotnie (inne położenie w skanerze);
• zupełnie różnych pacjentów;
• obraz CT ciała obejmujący 30 cm i preparat o długości 5 cm.

Tablica nie zawiera żadnej z tych informacji. Zawiera liczbę wymiarów, kształt, typ danych i wartości — ale nie milimetry, orientację, tożsamość pacjenta ani informację o tym, co jest górą. Aby rozumować przestrzennie, trzeba przypisać jej geometrię.

Ten sam indeks wskazuje inny punkt fizyczny, gdy zmienia się początek układu, odstępy lub kierunek; sam kształt nie zapewnia wyrównania w przestrzeni pacjenta.

Trzy układy współrzędnych

Obraz medyczny funkcjonuje niejawnie w trzech układach współrzędnych, a ich mylenie jest źródłem większości błędów geometrii.

Array / voxel indices
(col i, row j, slice k)
integers, shape W,H,D
        |
        | affine matrix
        v
Physical / patient coords
(x, y, z) in mm, LPS
continuous, real numbers
        |
        | window/level, zoom
        v
Display coords
screen pixels
view-dependent

Pierwszym układem jest sama tablica: pixel_array[k, j, i], gdzie kolejność — w przyjętej w tym rozdziale konwencji odwzorowania indeksów na pacjenta — odczytuje się jako (przekrój, wiersz, kolumna). Tablicę indeksują liczby całkowite, a jej kształt ma postać (slices, rows, columns). Indeksy nie mówią nic o rozmiarze fizycznym. Przycięcie tablicy lub zmiana jej próbkowania zmienia indeksy, a wraz z nimi odwzorowanie indeksów na przestrzeń fizyczną, chyba że początek układu, odstępy, kierunek i rozmiar zostaną poprawnie przeniesione lub obliczone ponownie, tak aby każde zamierzone położenie fizyczne pozostało skojarzone z reprezentującym je wokselem. Drugim układem jest fizyczna przestrzeń pacjenta, mierzona w ciągłych milimetrach. DICOM definiuje układ współrzędnych pacjenta jako LPS: +x → Left, +y → Posterior, +z → Head — jest to skrót odnoszący się do zwykłych danych człowieka w pozycji BIPED. Wartości ImagePositionPatient są współrzędnymi w milimetrach LPS, natomiast wartości ImageOrientationPatient są bezwymiarowymi cosinusami kierunkowymi, które wskazują, dokąd w tej przestrzeni zwrócone są osie wierszy i kolumn obrazu; nie są podane w milimetrach. Dopiero tutaj można porównywać geometrię różnych obiektów. Trzeci układ — współrzędne wyświetlania, piksele ekranu, zależne od widoku — jest interfejsem interakcji: pomiary fizyczne definiuje się w fizycznej przestrzeni pacjenta lub obrazu, a współrzędne ekranu zwykle służą użytkownikowi do wyboru punktów w przestrzeni pacjenta, które chce zmierzyć. Dalej nie będziemy zajmować się tym interfejsem.

Z tego rozróżnienia wynikają dwie pułapki graniczne. 3D Slicer działa wewnętrznie w układzie RAS (Right-Anterior-Superior) i przelicza LPS↔RAS na styku z danymi zewnętrznymi; pominięcie tego przeliczenia albo wykonanie go dwukrotnie odbija lub przemieszcza współrzędne wzdłuż osi x i y. NumPy w ogóle nie ma układu współrzędnych — ma tylko osie tablicy. Właśnie dlatego sama ndarray nie ma znaczenia geometrycznego, dopóki nie przypisze się jej początku układu, odstępów i kierunku. SimpleITK pozostaje w LPS, zgodnie z DICOM, dlatego służy jako pomost między tablicami a przestrzenią fizyczną.

Pierwszy i drugi układ łączy macierz afiniczna, którą jawnie zbudujemy w następnej sekcji.

Geometria składa się z czterech elementów

Pełną regularną siatkę obrazu opisują cztery wielkości. Początek układu, odstępy i kierunek definiuje się w przestrzeni pacjenta, natomiast rozmiar określa zasięg przestrzeni indeksów, z której prowadzi odwzorowanie:

• Początek układu — współrzędne fizyczne woksela o indeksie (0,0,0).
• Odstępy — krok między środkami sąsiednich próbek: (sx, sy, sz) w mm wzdłuż trzech osi. Jest to odległość między sąsiednimi punktami próbkowania, a nie w każdym przypadku fizyczny obszar — objętość tkanki — zajmowany przez woksel.
• Kierunek — macierz 3×3, której kolumny są cosinusami kierunkowymi trzech osi wolumenu w przestrzeni pacjenta. Macierz jednostkowa oznacza, że osie wolumenu są zgodne z wybranymi osiami fizycznymi, niekoniecznie z osiami samego skanera; macierz niejednostkowa może opisywać skośność (pochylenie gantry), permutację osi lub odbicie.
• Rozmiar — liczba próbek wzdłuż każdej osi: (slices, rows, columns), która wyznacza zakres istniejących indeksów.

Te cztery wielkości tworzą cały słownik geometrii obrazów medycznych. Początek układu odpowiada na pytanie gdzie, odstępy — jak daleko od siebie, kierunek — w którą stronę, a rozmiar — ile. Każda operacja w tym rozdziale — porządkowanie przekrojów, porównywanie CT z SEG, zmiana próbkowania — polega na operowaniu tymi czterema wielkościami.

„Takie same wymiary” niczego nie dowodzą

Warto zapamiętać ten fakt; wykorzystamy go w rozdziale dwukrotnie:

Dwa obrazy zajmują tę samą siatkę tylko wtedy, gdy ich rozmiar oraz początek układu + odstępy + kierunek są zgodne w granicach uprzednio określonych tolerancji. Jest to równość siatek. Porównywalność przestrzenna jest słabszą właściwością: dwie różne siatki w tym samym układzie pacjenta można nadal porównywać za pośrednictwem ich współrzędnych lub po zmianie próbkowania. Identyczny kształt tablicy niczego nie dowodzi w odniesieniu do żadnej z tych właściwości.

Najpierw zastosujemy tę zasadę, aby rozstrzygnąć, czy segmentacja już współdzieli siatkę CT, a zatem nie wymaga zmiany próbkowania (sekcja o wspólnej przestrzeni), a następnie, aby ustalić, czy dwa wolumeny są współrejestrowane. Należy wyrobić sobie następujący odruch: przed każdą operacją woksel po wokselu na dwóch obiektach najpierw porównaj ich geometrię, a nie kształt.

Stawka staje się konkretna przy odejmowaniu. Załóżmy, że obliczasz diff = array_A - array_B dla dwóch obrazów CT o kształcie (100, 512, 512) i otrzymujesz wiarygodnie wyglądające liczby. Zanim uznasz diff za rzeczywiste porównanie woksel po wokselu, trzeba odpowiedzieć na dwa pytania. Czy obrazy przedstawiają tego samego pacjenta? Obraz różnicowy dwóch różnych pacjentów jest pozbawiony znaczenia. Czy ich siatki są zgodne — mają ten sam rozmiar, początek układu, odstępy i kierunek w granicach uprzednio określonych tolerancji? Jeśli siatki się różnią, woksel [k, j, i] w A znajduje się w innym położeniu fizycznym niż woksel [k, j, i] w B, więc odejmowanie według tych samych indeksów łączy woksele zajmujące różne punkty fizyczne: jest to porównanie bez rejestracji, z niedopasowanymi indeksami, a nie samo w sobie zniekształcenie geometryczne, ponieważ żadnej wartości w tablicach nie przemieszczono ani nie poddano zmianie próbkowania. Dokładnie tak dzieje się po ponownym eksporcie tego samego pacjenta i tej samej serii z innymi odstępami — oba obrazy pozostają porównywalne poprzez współrzędne lub zmianę próbkowania, ale nie przez zestawienie identycznych indeksów. Sam kształt tablicy nie dostarczył żadnej informacji potrzebnej do rozstrzygnięcia tej kwestii.

Tablica nadal nie ma więc położenia. Po odczytaniu metadanych DICOM ma natomiast przepis, według którego można je wyznaczyć. Następna sekcja buduje ten przepis — przekształcenie afiniczne wraz z rozwiązanym przykładem liczbowym.

Źródła

• DICOM PS3.3, Image Plane Module (C.7.6.2) oraz Image Position and Image Orientation (C.7.6.2.1.1) — ImagePositionPatient, ImageOrientationPatient i odwzorowanie indeksu piksela na współrzędne pacjenta. (bieżące wydanie)
• SimpleITK, Fundamental Concepts — początek układu, odstępy, rozmiar, kierunek i odwzorowanie indeksu na punkt fizyczny. (dokumentacja)
• 3D Slicer, Coordinate systems — konwencje LPS i RAS oraz ich wzajemna relacja. (dokumentacja)