Pułapka wierszy i kolumn w DICOM

Zapisujesz arr[x, y], aby zaindeksować kolumnę x i wiersz y obrazu CT — a anatomia okazuje się transponowana. Albo odczytujesz PixelSpacing, zakładasz, że pierwsza liczba jest odstępem wzdłuż x, a zmiana próbkowania rozciąga obraz w niewłaściwym kierunku. Oba przypadki są tym samym błędem w dwóch przebraniach i oba wynikają z pomyłki, od której zaczął się ten rozdział: indeksy są położeniami w tablicy, ale nic nie wymusza zgodności nazw osi indeksów pomiędzy narzędziami. Nazewnictwo wierszy i kolumn w DICOM krzyżuje się samo ze sobą, a następnie krzyżuje się ponownie w zestawieniu z NumPy i SimpleITK. Ta sekcja ustala konwencje w jednym miejscu, aby transpozycja wydarzyła się na papierze, a nie w potoku przetwarzania.

Ten sam woksel w trzech systemach nazewnictwa

Trzy różne narzędzia opisują jeden fizyczny woksel za pomocą trzech różnych kolejności indeksów. Wszystkie wskazują ten sam woksel, ale inaczej go indeksują:

DICOM pixel       (i, j)     =  (column i, row j)
SimpleITK index   (x, y, z)  =  (column i, row j, slice k)   # NOT NumPy order
NumPy array       [k, j, i]  =  (slice k, row j, column i)
# ^ for the array returned by sitk.GetArrayFromImage on a scalar 3D image

NumPy i SimpleITK indeksują ten sam woksel w przeciwnych kolejnościach. W tym tkwi pułapka, z jednym ważnym ograniczeniem zakresu: samo NumPy nie przypisuje osiom tablicy znaczenia przekroju, wiersza ani kolumny. arr[k, j, i] jest kolejnością w skalarnej tablicy 3D utworzonej przez sitk.GetArrayFromImage(img) — konwencją konwersji, a nie właściwością każdej ndarray. Podobnie i i j w DICOM są liczonymi od zera indeksami kolumny i wiersza przesyłanej macierzy obrazu; słowa „wiersz” i „kolumna” oznaczają tutaj położenia w macierzy, a nie kierunki anatomiczne.

Powtarzający się symbol V oznacza jeden fizyczny woksel; między DICOM, NumPy i SimpleITK zmienia się wyłącznie jego zapis.

Krzyżowe powiązanie wyjaśnione wprost

Nieporozumienie ma dwie warstwy. Warstwa 1 — krzyżują się nazwy. W geometrii opisanej w poprzedniej sekcji:

• pierwsza trójka ImageOrientationPatient = kierunek wiersza = kierunek wzrostu indeksu kolumny i;
• pierwsza wartość PixelSpacing = odstęp między wierszami = krok wykonywany, gdy indeks wiersza j zwiększa się o 1.

Cosinus kierunkowy wiersza łączy się więc z odstępem między kolumnami (jest to krok wzdłuż i), a cosinus kierunkowy kolumny z odstępem między wierszami (krokiem wzdłuż j). Słowo „wiersz” występuje w obu określeniach, ale odnosi się do różnych osi. Nie jest to drobiazg terminologiczny — odwrócenie tych powiązań powoduje zmianę próbkowania lub rejestrację wzdłuż niewłaściwej osi. Wartość PixelSpacing równa [0.8, 0.7] pokazuje to wyraźnie: krok przy wzroście indeksu kolumny i o 1 wynosi 0.7 — jest to druga wartość, odstęp między kolumnami — ponieważ i rośnie wzdłuż wiersza. Użycie pierwszej wartości rozciąga lub ściska obraz przy zmianie próbkowania wzdłuż niewłaściwej osi. Ten cichy błąd wygląda niewinnie, dopóki nie nałoży się maski.

Warstwa 2 — kolejności indeksów krzyżują się między bibliotekami. Piksel DICOM o indeksach (i = column, j = row) po konwersji skalarnego obrazu 3D przez sitk.GetArrayFromImage(img) znajduje się w tablicy NumPy, w której przekroje są pierwszym wymiarem, pod adresem numpy[k, j, i], czyli numpy[slice, row, column]. arr[0] jest zatem przekrojem 0, a w obrębie przekroju zapis arr[k, j, i] wskazuje wiersz j i kolumnę i: jest to transpozycja naiwnego zapisu arr[x, y].

DICOM: pixel (i=col, j=row)
        | same voxel
        +--> NumPy: arr[k, j, i] = slice, row, column
        |
        +--> SimpleITK: GetPixel(x=i, y=j, z=k)
             = column, row, slice

NumPy and SimpleITK: opposite index order

W razie wątpliwości pamiętaj:

W DICOM i oznacza kolumnę, a j wiersz; PixelSpacing ma kolejność [row, col]; NumPy używa [k, j, i] — w skalarnej tablicy 3D zwracanej przez GetArrayFromImage — natomiast SimpleITK (x=i, y=j, z=k).

Wersja operacyjna brzmi: zawsze wyświetlaj GetSpacing()/GetOrigin()/GetDirection(), zamiast przyjmować założenia, i jawnie wykonuj konwersję za pomocą sitk.GetArrayFromImage / GetImageFromArray na granicy NumPy↔SimpleITK. GetImageFromArray tworzy nowy obraz z domyślnym początkiem układu, odstępami i kierunkiem: zanim uznasz, że przestrzennie odpowiada źródłu tablicy, skopiuj albo jawnie ustaw zamierzony początek układu, odstępy i kierunek.

Konkretna kontrola, którą możesz wykonać

import SimpleITK as sitk
img = sitk.ReadImage("series_dir")          # or build from DICOM
arr = sitk.GetArrayFromImage(img)           # -> NumPy [k, j, i]

print(img.GetSize())                        # (x=cols, y=rows, z=slices)  -- SimpleITK order
print(arr.shape)                            # (slices, rows, cols)        -- NumPy order (reversed)
print(img.GetSpacing())                     # (sx, sy, sz) in mm, LPS

Zwróć uwagę na odwrócenie kolejności: img.GetSize() ma postać (x, y, z), natomiast arr.shape ma postać (z, y, x). Woksel o indeksie SimpleITK (x=10, y=20, z=5) znajduje się pod adresem arr[5, 20, 10]. Jeśli zapomnisz o odwróceniu kolejności, po cichu transponujesz wolumen. Jest to konwencja konwersji GetArrayFromImage dla skalarnego obrazu 3D — nie wewnętrzna reguła porządku pamięci ani indeksowania NumPy dla dowolnych tablic.

Jak rozpoznawać orientację na podstawie cosinusów kierunkowych

Dwie trójki IOP wskazują, w którą stronę względem pacjenta są zwrócone osie wierszy i kolumn. W standardowej orientacji osiowej wzrost indeksu kolumny i wskazuje lewą stronę pacjenta (+x w LPS — od prawej do lewej strony pacjenta), a wzrost indeksu wiersza j wskazuje kierunek posterior (+y). Daje to osiową wartość IOP wykorzystaną w przykładzie liczbowym z poprzedniej sekcji:

ImageOrientationPatient = [1, 0, 0,  0, 1, 0]   # row=+x(L), col=+y(P), axial plane

Seria czołowa lub strzałkowa ma inne cosinusy, natomiast seria skośna (pozyskana przy nachyleniu) ma cosinusy niewyrównane z osiami. Nigdy nie zakładaj, że seria jest osiowa — odczytaj cosinusy, ponieważ macierz kierunku stanowi część geometrii decydującej o tym, czy dwa obiekty współdzielą przestrzeń (sekcja o wspólnej przestrzeni).

Konwencje z tej sekcji zapewniają poprawność geometrii w płaszczyźnie. Zanim można zaufać wolumenowi 3D, trzeba rozstrzygnąć jeszcze jedną kwestię geometryczną: w jakiej kolejności układają się przekroje? Odpowiedzią ponownie jest geometria, a następna sekcja wyprowadza tę odpowiedź.

Ź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, Equation C.7.6.2.1-1) — semantyka indeksów pikseli (kolumna i, wiersz j), ImageOrientationPatient i kolejność PixelSpacing. (bieżące wydanie)
• SimpleITK, Fundamental Concepts — dualność indeksów i przestrzeni fizycznej, GetPixel, GetSize, GetSpacing oraz granica konwersji tablicy/obrazu NumPy↔SimpleITK z odwróconą kolejnością wymiarów. (dokumentacja)