Dlaczego obrazy i maski wymagają odmiennej zmiany próbkowania

Zmieniasz próbkowanie CT i jego maski guza, aby umieścić je na wspólnej siatce 1 mm. Do obu obiektów stosujesz ten sam interpolator. Maska zawiera teraz wartości wokseli 0.4 i 0.7 — „40% guza, 60% wątroby” — a cechy kształtu i granicy mogą się po cichu zmienić. Ta sekcja wyjaśnia, dlaczego zmiana próbkowania obrazu i maski jest celowo asymetryczna, co naprawdę robią interpolatory, dlaczego w potoku TRACE-CT ograniczanie zakresu HU następuje po zmianie próbkowania oraz dlaczego ta kolejność stanowi część definicji pomiaru. Asymetria jest kluczowym rozróżnieniem; kolejność operacji prowadzi do kolejnych rozdziałów o pomiarach (rozdział 5).

Wartości ciągłe a etykiety dyskretne

Zmiana próbkowania jest nieunikniona w radiomice, gdzie potoki normalizują dane do wspólnych odstępów, oraz podczas przenoszenia masek. Tę samą operację trzeba jednak wykonać inaczej dla obrazu i maski ze względu na to, co oznaczają zapisane w nich liczby.

Obraz zawiera wielkość ciągłą — HU, czyli pole skalarne o przybliżonym znaczeniu fizycznym — dlatego interpolacja między wartościami sąsiednich wokseli jest uzasadniona: rzeczywista wartość w punkcie między dwoma wokselami może wiarygodnie leżeć między ich wartościami.

Maska zawiera etykiety dyskretne. Interpolowanie „połowy guza i połowy wątroby” między dwoma wokselami z etykietami daje wartość, która nie jest żadną z tych etykiet — ma sens tylko wtedy, gdy właśnie zamierzasz uzyskać reprezentację ułamkową. Niektóre formaty celowo przechowują zajętość lub prawdopodobieństwo dla każdego woksela i jest to poprawna reprezentacja; błędem jest dopuszczenie, aby binarna maska trafiła do ekstrakcji radiomicznej jako przypadkowe etykiety ułamkowe. Dlatego próbkowanie obu obiektów zwykle zmienia się za pomocą różnych interpolatorów:

• Obrazy: interpolacja liniowa (lub B-spline) — gładka i fizycznie uzasadniona.
• Maski: metoda najbliższego sąsiada — każdy woksel wyjściowy otrzymuje etykietę najbliższego woksela wejściowego, dzięki czemu istniejący dyskretny zbiór etykiet zostaje zachowany bezpośrednio. Jest to również ścieżka masek w PyRadiomics.

Reguła uwzględniająca etykiety: metoda najbliższego sąsiada zachowuje dyskretny zbiór etykiet bez dodatkowego etapu. IBSI dopuszcza ponadto liniową (trójliniową) interpolację maski, po której następuje jawna binaryzacja z określonym progiem — jest to poprawna alternatywa, jeżeli próg zostanie zadeklarowany i przypięty. Nie wolno pozostawić binarnej mapy etykiet w postaci ułamkowej po interpolacji: narusza to integralność etykiet i może zmienić cechy zależne od granicy.

L oznacza interpolację liniową, a NN — przeniesienie etykiety metodą najbliższego sąsiada. Znacznik wartości ułamkowej pokazuje, czego nie wolno przekazać dalej jako etykiety binarnej; liniowa lub trójliniowa interpolacja maski również jest poprawna, jeżeli następuje po niej jawna ponowna binaryzacja z zadeklarowanym i przypiętym progiem.

Porównanie interpolatorów

Interpolator	Zastosowanie	Działanie	Ryzyko
Najbliższy sąsiad	maska / etykiety	każdy woksel wyjściowy otrzymuje etykietę najbliższego woksela wejściowego; etykiety całkowitoliczbowe zostają zachowane	schodkowa granica (zgodnie z oczekiwaniem)
Liniowy (trójliniowy)	obraz (bezpieczny wybór domyślny)	wypukła kombinacja ważona 8 otaczających wokseli — pozostaje w ich zakresie; metoda szybka, bez przekroczenia zakresu	wygładzenie (jego stopień zależy od kontekstu)
B-spline	obraz (zamrożony wybór TRACE-CT)	gładka, kawałkami wielomianowa rekonstrukcja	możliwe przekroczenie zakresu wartości wejściowych przy ostrych granicach

Precyzyjna definicja siatki

Siatkę wokseli definiują trzy wielkości z otwarcia rozdziału: początek układu (fizyczne położenie woksela [0,0,0]), odstępy (krok fizyczny wzdłuż każdej osi wokseli — przypadek TRACE-CT RIDER zestawia serie zrekonstruowane z grubością przekroju 1.25 mm i 5.0 mm) oraz kierunek (orientacja osi w przestrzeni pacjenta). Zmiana próbkowania oznacza wyznaczenie wartości obrazu na nowej siatce — o innych odstępach, początku układu lub kierunku. Aby obliczyć wartość każdego nowego woksela, odczytuje się woksele wejściowe otaczające odpowiadające mu położenie fizyczne i łączy ich wartości. Kluczowe pytanie brzmi, jak je połączyć, a odpowiedź zależy od tego, co oznaczają liczby: pole ciągłe czy etykiety dyskretne.

Gęstsza siatka nie oznacza większej rozdzielczości

Przypadek TRACE-CT RIDER zestawia rekonstrukcję o grubości 5 mm z rekonstrukcją o grubości 1.25 mm u tego samego pacjenta — są to dwie różne rekonstrukcje, a nie dwa próbkowania jednego identycznego obrazu. Umieszczenie serii 5 mm na siatce 1 mm daje w przybliżeniu pięć razy więcej miejsc próbkowania wzdłuż z — dokładna liczba zależy od obrazowanego zasięgu i konwencji punktów końcowych — ale nie dostarcza nowej anatomii. Więcej próbek ≠ więcej informacji. Siatka 1 mm wyznacza wartości interpolacji próbek 5 mm; nie może odtworzyć szczegółów o wysokiej częstotliwości wzdłuż z, które rekonstrukcja 5 mm utraciła już podczas samej rekonstrukcji (przypadek rekonstrukcji z rozdziału 1). Żaden kierunek zmiany próbkowania nie tworzy rozdzielczości pozyskanej podczas akwizycji i oba mogą zmieniać wartości: zwiększenie częstości próbkowania wypełnia gęstszą siatkę wartościami wywnioskowanymi, które mogą być wygładzone, przekraczać zakres lub w inny sposób różnić się od pomiaru natywnego; zmniejszenie częstości próbkowania usuwa szczegóły i może rozmywać lub powodować aliasing. Zmiana próbkowania serii 1.25 mm do wspólnej izotropowej siatki 1 mm w TRACE-CT jest decyzją projektową dotyczącą spójności potoku — nie zmianą, która pozostawia każdą wartość nietkniętą.

Co naprawdę robi B-spline i dlaczego może przekraczać zakres

• Najbliższy sąsiad kopiuje najbliższy woksel wejściowy. Wynik jest nieciągły, ale bezpieczny dla etykiet.
• Interpolacja liniowa lub trójliniowa oblicza wypukłą kombinację ważoną ośmiu otaczających wokseli — w 3D są to 8 wierzchołków otaczającego sześcianu, a wagi wynikają z interpolacji trójliniowej. Wynik nigdy nie może wyjść poza ich zakres: wypukła kombinacja liczb jest ograniczona przez ich minimum i maksimum. Stopień wygładzenia zależy od kontekstu — próbkowania, struktur i liczby połączonych operacji.
• B-spline rekonstruuje obraz jako gładki model kawałkami wielomianowy (funkcję sklejaną o zwartym nośniku — każdy fragment wielomianowy „widzi” tylko kilka wokseli), a następnie wyznacza wartości tego modelu w nowych położeniach. Kuszące jest wyobrażenie sobie surowych próbek jako lokalnych współczynników sterujących funkcji sklejanej, lecz należy traktować to wyłącznie jako intuicję: implementacje wyprowadzają współczynniki wielomianowe z próbek, a szczegóły różnią się między bibliotekami. Bezpieczne stwierdzenie brzmi: woksel wyjściowy nie jest średnią danych wejściowych, lecz wartością zbudowanego z nich gładkiego, kawałkami wielomianowego modelu.

Na tym opiera się cała intuicja: metoda liniowa odpowiada na pytanie „jaka jest tutaj średnia ważona?”, a B-spline — „jaką wartość wskazuje tutaj gładki, kawałkami wielomianowy model zbudowany z mojego sąsiedztwa?”. Interpolacja gładka wyższego rzędu nie jest ograniczona do lokalnego zakresu próbek, dlatego przy ostrych przejściach może go przekraczać — na granicy między powietrzem o wartości −1000 HU a tkanką o wartości +100 HU może dać wartości poniżej −1000 lub powyżej +100. Przypomina to elastyczną linijkę wyginaną wokół narożnika: dopasowana krzywa może wyjść poza zakres wyznaczony przez dwie próbki. To właśnie przekroczenie zakresu — właściwość gładkiej interpolacji wyższego rzędu na ostrych krawędziach, a nie błąd interpolatora. Wybór B-spline jest decyzją projektową, nie prawem: w TRACE-CT jest to zamrożony interpolator obrazu, wybrany ze względu na wierność w gładkich obszarach, ale w innych zastosowaniach równie zasadnie używa się interpolacji liniowej lub innych metod. Wymogiem naukowym jest uzasadnienie, zamrożenie i zaraportowanie interpolatora, ponieważ determinuje on cechy, a jeśli wybrano B-spline — również poradzenie sobie z przekroczeniem zakresu. Temu służy opisane niżej ograniczanie wartości.

Ograniczanie wartości: przykład liczbowy

TRACE-CT ogranicza wartości HU do [-1000, +400] po zmianie próbkowania:

input HU      clipped HU
−1200    →    −1000     (below range: pinned to floor)
−850     →    −850      (inside range: unchanged)
 120     →     120      (inside range: unchanged)
 900     →     400      (above range: pinned to ceiling)

Ograniczanie zakresu modyfikuje intensywności obrazu i nic więcej — maska, liczba wokseli i geometria pozostają nienaruszone. Każda wartość woksela poniżej −1000 staje się równa −1000, a każda wartość powyżej +400 staje się równa +400. Jedna operacja daje dwa skutki: (1) wymusza uprzednio określone granice analizy, co może zmienić obserwowany rozkład intensywności i jego dyskretyzację; (2) ustala zakres intensywności dostępny podczas dyskretyzacji, co ma znaczenie, ponieważ cechy tekstury dyskretyzowane z użyciem binWidth zależą od tego zakresu (rozdział 6).

Ograniczanie zakresu nie jest ponowną segmentacją. PyRadiomics udostępnia oddzielne ustawienie resegmentRange, które usuwa z maski woksele o intensywnościach spoza określonego zakresu — zmienia zatem przynależność wokseli do ROI. Ograniczanie zakresu zmienia wartości wnoszone przez te woksele. TRACE-CT ogranicza intensywności obrazu poza PyRadiomics i pozostawia resegmentRange ustawione na None, co oznacza brak zmiany maski. Pomylenie tych operacji po cichu zmienia zarówno rozkład intensywności, jak i ROI.

Kolejność operacji stanowi część definicji pomiaru

Czy resample → clip jest tym samym co clip → resample? Nie — te operacje nie są algebraicznie równoważne, co pokazuje krótki eksperyment myślowy. Rozważ jednowymiarową granicę z próbkami wejściowymi …, −1000, −1000, +350, +500, …:

• Najpierw ograniczanie: +500 → +400; próbki mają teraz wartości …, −1000, −1000, +350, +400, …. Następnie interpolujemy w punkcie środkowym między +350 i +400: metoda liniowa daje +375, a B-spline — pobliską wartość wynikającą z gładkiego dopasowania. Ograniczenie wartości przed interpolacją zmienia próbki dostępne dla interpolatora — struktura wartości sprzed ograniczenia poza tym zakresem (czy było to +500, czy +900?) zostaje usunięta, zanim interpolator w ogóle ją zobaczy. Nie ogranicza to wyniku do [-1000, +400]: interpolacja wyższego rzędu, taka jak B-spline, nadal może przekroczyć zakres próbek, które już ograniczono.
• Najpierw zmiana próbkowania: interpolator otrzymuje wartości bez ograniczenia (+350, +500) i oblicza gładkie oszacowanie między nimi oraz w ich otoczeniu (+425 w punkcie środkowym, z możliwością przekroczenia +500). Dopiero potem ograniczenie sprowadza wartości końcowe do [-1000, +400] (+425 → +400).

Te same dane wejściowe i ta sama siatka docelowa dają inne zmierzone intensywności na granicy, dlatego obliczona z nich cecha tekstury również może się różnić. Żadna kolejność nie jest „błędna”; są to różne definicje pomiaru. W TRACE-CT zamrożono następującą kolejność:

DICOM pixels → HU (RescaleSlope/Intercept)
  → external 1 mm isotropic resampling (B-spline image, NN mask)
  → clip to [-1000, +400] HU
  → PyRadiomics

i zamrożono ją z podanych powodów: jest to zatwierdzona implementacja i konfiguracja potoku TRACE-CT, sprawdzana testami regresji w przypiętym commicie (testy obejmują kolejność zmiana próbkowania, potem ograniczanie, oraz konfigurację). Ta walidacja ustanawia, czym jest potok — nie stanowi twierdzenia o poprawności klinicznej, optymalności ani o wystąpieniu przekroczenia zakresu B-spline w każdej serii. Kolejność gwarantuje natomiast, że ograniczanie obejmuje końcowe wartości próbek bezpośrednio przed ekstrakcją. Jest to zamrożona definicja pomiaru TRACE-CT, a nie twierdzenie, że ograniczanie po zmianie próbkowania jest zawsze jedyną poprawną kolejnością. Uniwersalna jest dyscyplina: określ kolejność, przypnij ją i nie zmieniaj uprzednio określonego głównego potoku po zobaczeniu wyników bez zadeklarowania nowej analizy — prospektywna zmiana metody jest uzasadniona, ale tworzy nową definicję pomiaru, a nie tę samą, ponieważ zmiana kolejności zmienia mierzoną wielkość. Jest to zabezpieczenie przed dostrajaniem, później ujęte jako studium przypadku „najpierw mierz, potem interweniuj” w rozdziale 6.

Bez podwójnej zmiany próbkowania

Z tej samej logiki wynika jeszcze jedna zamrożona zasada. TRACE-CT zmienia próbkowanie na 1 mm zewnętrznie, przed wywołaniem PyRadiomics. Dlatego plik parametrów PyRadiomics nie może ustawiać resampledPixelSpacing. W przeciwnym razie PyRadiomics ponownie zmieniłby próbkowanie obrazu, który już przeszedł tę operację, i zastosowałby do danych drugi operator interpolacji. Dwa przebiegi interpolacji nie stanowią tego samego pomiaru co jeden: dodatkowy etap zmienia potok i może zmienić próbki, a działanie zależy od interpolatora — metoda najbliższego sąsiada w ogóle nie wygładza, a metoda liniowa nie przekracza zakresu tak jak może to robić B-spline. Jest to zamrożona zasada projektowa TRACE-CT: zewnętrzna zmiana próbkowania albo zmiana próbkowania w PyRadiomics, nigdy obie. Stanowi decyzję projektową zapobiegającą niezamierzonemu powtórzeniu etapu, a nie uniwersalne prawo zakazujące każdemu potokowi dwóch zmian próbkowania. TRACE-CT sprawdza tę zasadę, wczytując plik parametrów przez faktycznie używany ekstraktor w przypiętej wersji i kontrolując, czy efektywna wartość resampledPixelSpacing wynosi None (sekcja o definicji operacyjnej z rozdziału 5).

Asymetryczna zmiana próbkowania w kodzie

import SimpleITK as sitk
# image: B-spline; mask: nearest-neighbour — to the SAME reference geometry
resampled_img  = sitk.Resample(img,  ref_img, sitk.Transform(), sitk.sitkBSpline,
                               0.0, sitk.sitkInt16)
resampled_mask = sitk.Resample(mask, ref_img, sitk.Transform(), sitk.sitkNearestNeighbor,
                               0, sitk.sitkUInt8)

Zwróć uwagę na celową asymetrię: ten sam obraz odniesienia, to samo przekształcenie, ale inny interpolator i inna wartość tła dla obrazu oraz maski. W ten sposób pomost SimpleITK z otwarcia rozdziału spełnia swoje zadanie — przenosi geometrię razem z obrazem, dzięki czemu zmiana próbkowania jest operacją geometryczną, a nie przestawianiem elementów tablicy.

Dwa znaczenia efektu objętości częściowej: fizyka akwizycji a cyfrowa zmiana próbkowania

Określenie efekt objętości częściowej jest używane w odniesieniu do dwóch różnych zjawisk. Ich rozdzielenie pozwala zachować poprawny opis przyczynowy.

Fizyczny efekt objętości częściowej na etapie akwizycji i rekonstrukcji. Zrekonstruowany woksel lub przekrój reprezentuje skończony zasięg przestrzenny — rekonstrukcja o grubości 5 mm ma skończoną czułość w kierunku prostopadłym do płaszczyzny, rzędu grubości przekroju, a nie czułość nieskończenie cienkiej płaszczyzny. Dokładny profil czułości zależy od akwizycji i rekonstrukcji; nie musi być jednorodną średnią z obszaru o grubości 5 mm. Tkanki znajdujące się w tym profilu czułości wpływają na zrekonstruowaną wartość, dlatego grube przekroje mieszają struktury i tracą szczegóły przestrzenne. Co najważniejsze, następuje to wewnątrz skanera i procesu rekonstrukcji: zanim powstanie obraz cyfrowy, informacja jest już utracona. Żadna późniejsza zmiana próbkowania, niezależnie od kierunku, nie może jej odzyskać.

Cyfrowa zmiana próbkowania, w tym zmiana na rzadszą siatkę. Operacja ta działa na już zrekonstruowanym obrazie cyfrowym: interpolacja estymuje wartości na innej siatce. Zmiana próbkowania na rzadszą siatkę może usuwać szczegóły, wygładzać, powodować aliasing albo rozmywać granice — dokładny efekt zależy od interpolatora (metoda najbliższego sąsiada niczego nie uśrednia, lecz wybiera woksel źródłowy; B-spline i metoda liniowa obliczają gładkie oszacowania). Jest to operacja matematyczna na próbkach, a nie ponowny pomiar mieszaniny tkanek.

anatomy
  → scanner acquisition + reconstruction (physical partial volume: finite
    voxel/slice support; thick reconstructions mix structures — loss happens HERE)
  → reconstructed voxel values (the digital image)
  → later digital interpolation / resampling (grid change; interpolator-dependent)

Kolejność ma znaczenie, ponieważ przypisuje nieodwracalną utratę do właściwego etapu: zmiana próbkowania serii 5 mm na siatkę 1 mm nie może odtworzyć szczegółów w kierunku prostopadłym do płaszczyzny, które rekonstrukcja 5 mm już utraciła — jedynie ponownie wyznacza wartości zachowanych próbek na gęstszej siatce (zob. gęstsza siatka nie oznacza większej rozdzielczości powyżej). Cyfrowa zmiana próbkowania na rzadszą siatkę wprowadza natomiast własną, zależną od interpolatora zmianę ostrości granic i tekstury — na granicy guza rzadsza siatka sprawia, że reprezentacja granicy staje się bardziej rozmyta albo schodkowa, co zmienia objętość, ostrość krawędzi i cechy tekstury. Ten skutek na dalszych etapach jest oczekiwanym działaniem potoku, a nie błędem. Właśnie dlatego każda analiza radiomiczna obejmująca zmianę próbkowania musi, odpowiednio do przypadku, raportować siatkę docelową (rozmiar, początek układu, odstępy, kierunek), przekształcenie, interpolator, mające zastosowanie wartości domyślne i kolejność operacji — bez tego cecha jest nieodtwarzalna.

Główny wniosek rozdziału

Wszystko w tym rozdziale prowadzi do jednego wniosku:

Wartość radiomiczna jest funkcją próbek obrazu, reprezentacji ROI, pełnej geometrii, przetwarzania wstępnego i implementacji cechy. Interpolacja i dyskretyzacja są wyborami dotyczącymi przetwarzania wstępnego, a nie częściami geometrii. Zmiana któregokolwiek z tych elementów zmienia definicję pomiaru — i wartość cechy może się zmienić, nawet jeśli anatomia pozostaje niezmieniona.

Krótki wariant wyostrza ten wniosek. Zmień próbkowanie obrazu i maski z 1 mm na izotropową siatkę 2 mm, stosując metodę najbliższego sąsiada do maski i B-spline do obrazu, a następnie ponownie oblicz cechę tekstury: jej wartość się zmienia. Nowa wartość nie jest „błędna” — jest cechą w innej geometrii, a zatem zasadnie inną wielkością. Błędem byłoby podanie nowej wartości bez zarejestrowania nowych odstępów i interpolatora albo porównanie jej z wartością obliczoną przy 1 mm tak, jakby oba wyniki stanowiły ten sam pomiar. Zmiana próbkowania może zmienić wartości cech i zawsze zmienia definicję pomiaru, nawet jeśli konkretna wartość pozostanie taka sama. Wymagana dyscyplina polega na przypięciu i zaraportowaniu geometrii docelowej, aby każde porównanie dotyczyło wielkości wyznaczonych w takich samych warunkach. Stwierdzenie, że cecha „zmieniła się wskutek zmiany próbkowania”, opisuje problem tylko wtedy, gdy nie zauważysz zmiany potoku.

Dlatego IBSI ustala przepływ przetwarzania, PyRadiomics pozwala przypiąć parametry zmiany próbkowania, a TRACE-CT zatwierdza plik parametrów (configs/pyradiomics_params.yaml) z ustalonymi wartościami binWidth i odstępów, zamiast pozostawiać je ustawieniom domyślnym. Cecha nie jest właściwością guza; jest właściwością guza w ramach określonego potoku. Kolejne rozdziały precyzują to stwierdzenie.

Ostatnim krokiem w tym rozdziale jest przerwać czytanie i spojrzeć: studium przypadku zmienia wszystkie powyższe pojęcia w elementy, które można sprawdzić wzrokiem.

Źródła

• Zwanenburg A, Vallières M, Abdalah MA, et al. The Image Biomarker Standardisation Initiative: standardized quantitative radiomics for high-throughput image-based phenotyping. Radiology. 2020;295(2):328-338. PMID 32154773. doi:10.1148/radiol.2020191145. (PubMed; opublikowany pełny tekst: PMC7193906) — źródło IBSI dotyczące algorytmów interpolacji i efektów objętości częściowej w masce ROI.
• Dokumentacja PyRadiomics — Customizing the extract: resampledPixelSpacing, interpolator, resegmentRange i pliki parametrów. (pyradiomics.readthedocs.io/en/latest/customization.html)
• Dokumentacja SimpleITK — Fundamental Concepts: Resampling oraz dokumentacja API ResampleImageFilter. (simpleitk.readthedocs.io/en/master/fundamentalConcepts.html#resampling)
• DICOM PS3.3, CT Image Module — Rescale Slope (0028,1053) i Rescale Intercept (0028,1052) wraz z równaniem jednostek wyjściowych. (dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.8.2.html)
• Zhao B, Schwartz LH, Kris MG, Riely GJ. Coffee-break lung CT collection with scan images reconstructed at multiple imaging parameters (Version 3) [Dataset]. The Cancer Imaging Archive. doi:10.7937/K9/TCIA.2015.U1X8A5NR. (doi.org/10.7937/K9/TCIA.2015.U1X8A5NR) — tożsamość i licencja publicznej kolekcji RIDER; etykiety grubości przekrojów 1.25 mm i 5.0 mm użyte w tej sekcji są celami projektu TRACE-CT, a nie ogólnymi ustawieniami zbioru danych.
• Potok radiomiczny TRACE-CT (ograniczony dowód z przypadku) — zamrożona implementacja przetwarzania wstępnego, plik parametrów (configs/pyradiomics_params.yaml) i testy dryfu w commicie 7aaf308c76ed2718fd0b364d1cde856969711b72. (commit 7aaf308)