Porządkuj przekroje bez ufania nazwom plików

Masz 320 plików DICOM z jednej serii. Ich nazwy nie niosą informacji przestrzennej, a w InstanceNumber występują luki. W jakiej kolejności należy je ułożyć, aby otrzymać prawidłowy wolumen 3D? Poprzednie sekcje jeszcze tego nie rozstrzygnęły: przekształcenie afiniczne umieściło każdy przekrój w ciele pacjenta, ale wolumen jest stosem, a ułożenie przekrojów w złej kolejności zrywa zgodność między indeksem tablicy a kierunkiem przemieszczania się w przestrzeni pacjenta — późniejsze obliczenia zależne od z mogą być wtedy po cichu błędne. Odpowiedź wyprowadzana w tej sekcji polega na obliczeniu kolejności z położenia fizycznego — z geometrii, a nie z etykiet.

Dlaczego InstanceNumber i nazwy plików kłamią

Normatywne znaczenie InstanceNumber (0020,0013) jest skromne: DICOM definiuje go jedynie jako „liczbę identyfikującą ten obraz” w obrębie serii. Nie jest to współrzędna przestrzenna i nic w standardzie nie nakazuje, aby kodowała kolejność w przestrzeni pacjenta potrzebną do utworzenia wolumenu — ciągłość, monotoniczność anatomiczna i stabilność po ponownym eksporcie leżą poza jej kontraktem, a wartości mają znaczenie wyłącznie w obrębie jednej serii. Nazwy plików niosą jeszcze mniej informacji przestrzennej: ponowny eksport lub pobranie może zmienić nazwy bądź kolejność plików, a zwykłe sortowanie leksykograficzne umieszcza slice_10 przed slice_2. Żadna z tych etykiet nie zawiera geometrii płaszczyzny obrazu, więc żadna nie może ustalić kolejności anatomicznej.

Jeżeli ułożysz przekroje według InstanceNumber lub nazw plików i dopisze Ci szczęście, wolumen będzie wyglądał prawidłowo. Jeśli nie, przekroje zostaną subtelnie pomieszane: indeks tablicy przestanie odpowiadać postępowi w przestrzeni pacjenta, a obliczenia zależne od geometrii — sąsiedztwa w z, powierzchnie, tekstury — mogą zostać po cichu zniekształcone. Nie każda wielkość musi się zmienić — na przykład liczba wokseli pomnożona przez stałą objętość woksela jest niezmiennicza względem permutacji — a tablica wynikowa w żadnym z przypadków nie sygnalizuje uszkodzenia. To otwierająca rozdział awaria maski w niewłaściwej przestrzeni, przesunięta o jeden wymiar. Rozwiązaniem jest odrzucenie zaufania do etykiet i wyprowadzenie kolejności z fizycznego położenia każdego przekroju.

Wektor normalny przekroju

Trzecią osią wolumenu CT jest wektor normalny przekroju n — kierunek, wzdłuż którego układają się przekroje. Jest on iloczynem wektorowym cosinusów kierunkowych wiersza i kolumny:

n = cross([rx, ry, rz], [cx, cy, cz])

Dla serii osiowej w standardowej orientacji r = [1, 0, 0], c = [0, 1, 0] otrzymujemy n równy [0, 0, 1] — w przybliżeniu kierunek Superior, czyli kierunek „w górę stosu”. Znak wynika z uporządkowanych cosinusów kierunkowych, a nie z chronologii akwizycji. Dowolna orientacja osiowa nie musi mieć normalnej zwróconej ku Superior — dlatego metoda rzutuje położenia na obliczony wektor n, zamiast zakładać +z. Wszystko w tej sekcji sprowadza się do rzutowania położeń przekrojów na n.

Porządkowanie według rzutu na wektor normalny przekroju

Każdy przekrój ma własną wartość ImagePositionPatient p_k — fizyczne położenie lewego górnego woksela, jak ustaliliśmy przy budowie przekształcenia afinicznego. Zrzutuj to położenie na wektor normalny przekroju, aby uzyskać pojedynczą skalarną „wysokość” wzdłuż kierunku stosu:

t_k = dot(p_k, n)

Następnie posortuj przekroje rosnąco według t_k. Dla równoległych przekrojów o wspólnej, poprawnej i spójnej orientacji otrzymasz stos uporządkowany konsekwentnie wzdłuż jednego kierunku fizycznego, niezależnie od nazw plików i InstanceNumber.

slice A  IPP = [..., 40.0]  t = 40.0
    |
    v
slice B  IPP = [..., 41.25]  t = 41.25
    |
    v
slice C  IPP = [..., 42.5]  t = 42.5

sort by t ascending = order along the stack axis

Rzutowanie i sortowanie są konieczne, lecz same nie wystarczają. Poprawny stos wymaga także sprawdzenia geometrii: cosinusy kierunkowe muszą być skończone, ortonormalne i spójne między przekrojami; rzuty położeń muszą być unikatowe — zduplikowane lub niemal zduplikowane położenia są niejednoznaczne i nie dają się wzajemnie uporządkować; wymiary w płaszczyźnie i odstępy między pikselami muszą być spójne; przekroje lokalizacyjne (scout) trzeba wykluczyć. Ponadto sortowanie unieszkodliwia ponowny eksport wyłącznie wtedy, gdy eksport zachował geometrię: porządkuje zakodowane położenia niezależnie od tego, czy są prawidłowe.

Dokładnie takie obliczenia wykonują funkcje compute_slice_normal() i project_position() w TRACE-CT, a load_ct_volume_explicit() wykorzystuje je do uporządkowania metadanych przekrojów; następnie loader sprawdza wymiary, orientację, odstępy między pikselami, unikatowość położeń i regularność odstępów, zanim przydzieli pamięć na wolumen pikseli i go wypełni.

Przykład rozwiązany krok po kroku: trzy przekroje

Jako syntetyczne wyprowadzenie — skonstruowany przykład, a nie zmierzoną serię — rozważmy trzy przekroje (osiowe, więc n ≈ [0,0,1]) o następujących położeniach:

Plik przekroju	ImagePositionPatient (z)	t = dot(p, n)
img_002.dcm	42.5	42.5
img_001.dcm	40.0	40.0
img_003.dcm	41.25	41.25

Nazwy plików wskazują kolejność 1, 2, 3, lecz położenia są przemieszane. Sortowanie według t daje fizyczną kolejność stosu: img_001 (40.0) → img_003 (41.25) → img_002 (42.5). Odstęp między kolejnymi przekrojami jest różnicą kolejnych wartości t: 1,25 mm i 1,25 mm — jest to regularny stos.

Ostrzegawczy wariant hipotetyczny pokazuje, dlaczego metoda geometryczna jest warta stosowania. Załóżmy, że inna osoba porządkuje poprawnie wyeksportowaną serię według InstanceNumber i otrzymuje prawidłowy wolumen — wręcz identyczny z wynikiem porządkowania geometrycznego. Taka zgodność dowodzi jedynie, że etykieta i geometria są zgodne w tej jednej serii; nie dowodzi ogólnego bezpieczeństwa etykiety. Porządkowanie geometryczne wymaga kilku wierszy kodu i eliminuje pytanie, czy tym razem etykieta przypadkiem odpowiada geometrii.

Nieregularne odstępy: kiedy stos nie jest wolumenem

Po posortowaniu wartości t_k oblicz różnice t_{k+1} − t_k. Jeśli wszystkie są w przybliżeniu równe — w granicach tolerancji określonej z góry odpowiednio do precyzji kodowania i dalszego zastosowania — stos może reprezentować jedną regularną trzecią oś. Jeśli się różnią, występują nieregularne odstępy: założenia regularnej siatki, na których opierają się operacje 3D woksel po wokselu, przestają obowiązywać, odległości wzdłuż z nie odpowiadają temu, co sugerują indeksy tablicy, a każdy pomiar lub każda zmiana próbkowania traktująca stos jako jednorodny będzie po cichu błędna. Nierówne odstępy mogą wynikać z brakujących przekrojów, akwizycji z przeplotem albo innych skutków akwizycji i selekcji; sam wzorzec odstępów nie rozstrzyga przyczyny, a w szczególności nie oznacza sam przez się zmiennej grubości przekroju. Równe odstępy także są warunkiem koniecznym, lecz niewystarczającym: regularna siatka wymaga ponadto spójnej orientacji, spójnych odstępów i wymiarów w płaszczyźnie, unikatowych położeń oraz skończonej geometrii. Możliwe są wtedy dwa rozwiązania:

• zmienić próbkowanie do regularnej siatki, jeśli uzasadnia to fizyka (sekcja o zmianie próbkowania);
• albo potraktować przekroje jako niejednorodny stos i odrzucić operacje 3D woksel po wokselu, które zakładają regularność.

Każde z tych rozwiązań jest uzasadnione dopiero po ocenie pochodzenia stosu i jego fizycznej adekwatności.

W obu przypadkach zmierz odstępy; nie zakładaj ich wartości — i wykonaj kontrolę krzyżową. Oblicz rzuty odstępów między środkami sąsiednich przekrojów i porównaj je z SpacingBetweenSlices (0018,0088), jeżeli seria deklaruje ten atrybut: jego normatywne znaczenie to odległość między środkami sąsiednich przekrojów, natomiast różnice rzutów są odstępami wynikającymi z zakodowanych położeń przy spójnej orientacji. Żadne z tych pól nie jest nieomylną prawdą podstawową; są to dwa zapisy, które mogą sobie przeczyć, a istotna niezgodność nakazuje zbadać lub odrzucić serię, nie zaś wybrać ulubioną wartość. Taką dyscyplinę stosuje load_ct_volume_explicit() w TRACE-CT: odrzuca zduplikowane położenia i nieregularne odstępy wynikające z rzutów, zanim zwróci użyteczny wolumen.

Po ustaleniu kolejności przekrojów i regularności siatki wolumen jest geometrycznie poprawnym samodzielnym obiektem. Następne pytanie dotyczy relacji: kiedy segmentacja naprawdę współdzieli przestrzeń CT — i czego właściwie dowodzi FrameOfReferenceUID?

Źródła

• DICOM PS3.3, General Image Module (C.7.6.1, Table C.7-9) — InstanceNumber (0020,0013) zdefiniowany jako liczba identyfikująca ten obraz. (bieżące wydanie)
• DICOM PS3.3, Image Plane Module (C.7.6.2) — semantyka ImagePositionPatient, ImageOrientationPatient i SpacingBetweenSlices (0018,0088), czyli odległości między środkami sąsiednich przekrojów. (bieżące wydanie)
• SimpleITK, Fundamental Concepts — początek układu, odstępy, kierunek i rozmiar jako metadane geometrii regularnego obrazu. (dokumentacja)
• TRACE-CT — implementacja porządkowania przekrojów i walidacji wolumenu CT (compute_slice_normal(), project_position(), load_ct_volume_explicit()), ograniczona do wskazanego commitu. (commit 7aaf308)