1. Wstęp
Automatyczne wykrywanie krawędzi jest fundamentalną kwestią w analizie obrazów. W obrazowaniu serca możliwość automatycznego wykrycia granicy wsierdzia w obrazowaniu lewej komory pozwoliłaby na obiektywny pomiar objętości komór oraz odkształcenia mięśnia sercowego (strain). Zostało to osiągnięte w echokardiografii dzięki technice śledzenia plam (speckle tracking). Opracowanie niezawodnych metod automatycznego wykrywania granic jest trudnym zadaniem, dla którego nie znaleziono jeszcze powszechnie uznanego i niezawodnego rozwiązania w rezonansie magnetycznym serca (CMR). W rzeczywistości w praktyce klinicznej granice są albo rysowane ręcznie przez operatora, albo wykrywane przez oprogramowanie jako interfejs między mięśniem sercowym a jamą 1,2. W niniejszym artykule przedstawiamy inne podejście, w którym granice nie są „wykrywane”, lecz „śledzone”, tzn. obserwowane w czasie, zaczynając od jednego wiarygodnego, istniejącego obrysu chwilowego, który zazwyczaj – choć niekoniecznie – jest rysowany ręcznie przez doświadczonego operatora na pojedynczej klatce obrazu. Poszczególne punkty tworzące taki pierwszy wiarygodny obrys są śledzone w czasie poprzez wyszukiwanie tych samych cech, które znajdują się w sąsiedztwie danego punktu w kolejnych klatkach. Śledzone cechy mogą być granicą jama-tkanka lub elementami anatomicznymi, które różnią się w obrębie tkanki. Są one wyszukiwane za pomocą metod największej wiarygodności w dwóch obszarach zainteresowania pomiędzy dwiema klatkami.
Lokalne przemieszczenie z klatki na klatkę jest równoważne ocenie lokalnej prędkości (stosunku przemieszczenia do odstępu czasu). Automatyczna ocena prędkości w danym punkcie jest określana na podstawie porównania przemieszczenia danych obrazowych wokół tego punktu w dwóch kolejnych klatkach. Metody te, w kilku różnych sformułowaniach, były stosowane w wielu dziedzinach badawczych. W zaawansowanej analizie obrazu zalicza się je do ogólnej kategorii znanej jako przepływ optyczny (Optical Flow) 3,4. W obrazowaniu echograficznym, gdy takie prędkości są wykorzystywane do śledzenia ruchu fizjologicznego, określa się je zazwyczaj jako śledzenie plamek (Speckle Tracking) 5, 6, ale znajdują one zastosowanie również w każdej innej metodzie obrazowania, np. w CMR, gdzie metody te nazywa się śledzeniem cech (feature tracking) lub śledzeniem granic (border tracking).
2. Materiały i metody
Metoda śledzenia cech
Brzeg wsierdziowy lub epikardialny obrazu cine 2D CMR jest ręcznie obrysowywany na jednej dowolnej klatce (patrz ryc. 1). Można również obrysować cechy śródmięśniowe. Taki brzeg definiuje się następnie jako sekwencję N punktów, określonych przez ich pary współrzędnych (Proszę podać tekst źródłowy do tłumaczenia.i,yiz i=1...NŚledzenie krawędzi odbywa się poprzez śledzenie każdego pojedynczego punktu; proces ten opiera się na algorytmie hierarchicznym w wielu skalach oraz na połączeniu technik śledzenia 1D, które gwarantują wyższą dokładność, i śledzenia 2D, które jest niezbędne do prawidłowego wykrywania cech rozciągniętych przestrzennie w 2D.
Aby w pierwszej kolejności uchwycić duże przemieszczenie geometryczne brzegu, śledzenie przeprowadza się w kierunku prostopadłym do samego brzegu, gdzie granica między jamą a tkanką jest najlepiej rozpoznawalna. Śledzenie w tym kierunku wykonuje się przy użyciu metody cięć transmuralnych w następujący sposób (patrz ryc. 3). Rysuje się linię przecinającą ścianę, przechodzącą przez dany punkt i prostopadłą do niej. Piksele pobrane wzdłuż linii transmuralnej są rozmieszczone w kolumnach, gdzie każda kolumna odpowiada jednej klatce sekwencji obrazów. W ten sposób ewolucja wzdłuż cięcia transmuralnego może zostać przedstawiona dla wszystkich chwil jednocześnie w reprezentacji dwuwymiarowej, w której jedną osią jest odległość wzdłuż linii, a drugą oś czas (patrz ryc. 2). Reprezentacja ta jest podobna do tzw. trybu M w echokardiografii, natomiast w CMR odpowiada ona bardziej funkcji „scout”. Aby poprawić jakość analizy, w przypadku obrazów słabej jakości z niskim stosunkiem sygnału do szumu, reprezentacja czasoprzestrzenna jest budowana z wykorzystaniem linii cięcia transmuralnego o grubości 5 pikseli. Następnie śledzenie brzegu jest wykonywane na tak uzyskanej mapie czasoprzestrzennej.
W drugim kroku, aby uwzględnić dwuwymiarowe przemieszczenie krawędzi, dla każdego punktu niezależnie przeprowadza się standardowe śledzenie 2D (w oparciu o przepływ optyczny) w ruchomym oknie o rozmiarze MxM, które jest zawsze wyśrodkowane na wcześniej oszacowanym punkcie krawędzi. Śledzenie 2D wykonuje się w dwóch etapach, w których połowa pierwszej estymacji jest wykorzystywana do wyśrodkowania ruchomych okien w drugim przejściu śledzenia. Następnie rozmiar okna jest redukowany z 32 do 16 w dwóch dodatkowych przejściach.
Aby zwiększyć dokładność śledzenia ruchu wzdłuż krawędzi, który służy do szacowania rotacji i torsji, śledzenie 1D jest wykonywane na obrazach czasoprzestrzennych zbudowanych z grubych przekrojów „równoległych” do zakrzywionej krawędzi (patrz rycina 3). W każdym punkcie, niezależnie, piksele pobrane wzdłuż ruchomej krawędzi, wycentrowane w punktach tej krawędzi, są układane w kolumnach, gdzie każda kolumna odpowiada jednej klatce z sekwencji obrazów. Aby poprawić jakość analizy i jak najlepiej uchwycić cechy krawędzi, linię przedłużono o 5 pikseli w głąb tkanki (podwsierdzie). Śledzenie krawędzi jest następnie wykonywane na obrazie czasoprzestrzennym zgodnie z opisaną powyżej procedurą. Aby zapewnić spójność przestrzenną śledzonej krawędzi, do przesunięć obliczonych w sąsiednich punktach na każdym etapie zastosowano medianowy filtr 3-punktowy oraz gaussowski filtr 3-punktowy (o wagach 0,25, 0,5, 0,25).
Śledzenie wzdłuż dwuwymiarowego obrazu czasoprzestrzennego
W tej sekcji opisano procedurę śledzenia krawędzi wzdłuż jednego kierunku na obrazie dwuwymiarowym (zbliżonym do trybu M), rozpoczynając od znanej pozycji w danym momencie.
Oś X definiuje kierunek poziomy, a oś y kierunek pionowy. Kolumny są oznaczone jako xi, i=1...M, gdzie M to liczba kolumn na obrazie. Śledzenie polega na wyznaczeniu dyskretnej sekwencji liczb rzeczywistych yi=y(xi), zaczynając od znanego punktu yk odpowiadającego kolumnom xk.
Przesunięcie z punktu znanego yk do punktu yk+1 jest szacowane poprzez ocenę korelacji wzajemnej między całą kolumną w xk a całą kolumną w xk+1. Funkcja korelacji wzajemnej wykaże maksimum, którego pozycja wskazuje wartość przesunięcia pionowego niezbędnego do zmaksymalizowania podobieństwa między dwiema kolumnami; zatem yk+1 jest szacowane poprzez dodanie takiego przesunięcia do yk. Procedura ta jest powtarzana dla wszystkich par sąsiednich kolumn, a wynikiem jest oszacowanie całej krawędzi yi, i=1...M. Korelacja wzajemna jest tutaj obliczana przy użyciu algorytmu szybkiej transformaty Fouriera (Fast Fourier Transform) w celu skrócenia czasu obliczeń.
Pierwsza estymata yi jest następnie doprecyzowana iteracyjnie. Aby osiągnąć ten cel, wyodrębnia się podzbiór obrazu, pobierając kilka punktów powyżej i poniżej poprzedniej estymaty yi, a następnie generuje się nowy obraz, którego centrum odpowiada sekwencji yi, i wykorzystuje go do śledzenia korekcji. Proces doprecyzowania jest powtarzany do momentu, aż nie zostanie stwierdzona żadna korekta.
Ulepszony i bardziej naturalny wynik uzyskuje się następnie za pomocą końcowej procedury „snake” [5], aby w obrazie czasoprzestrzennym śledzić poziom jasności obrazu przechodzący przez punkt stały yk. Cały proces wykorzystuje okresowość czasową, aby zapewnić okresowy wynik i uniknąć efektu dryftu.
Ograniczenia techniczne śledzenia cech
Technika śledzenia krawędzi, podobnie jak każda metoda śledzenia plamek (speckle tracking), opiera się na kwantyfikacji zmian jasności pikseli między poszczególnymi klatkami. Narzuca to dolną granicę prędkości, wynikającą z konieczności zaobserwowania plamki, która w jednej klatce znajduje się w jednym pikselu, a w następnej przesuwa się do sąsiedniego piksela. Granica ta jest zatem określona wzorem:

Równanie 1
gdzie Δx to rozmiar piksela, a Δt to odstęp czasu między dwiema klatkami. Współczynnik k zależy od jakości algorytmu śledzenia oraz jego zdolności do oceny dynamicznych zmian w skali subpikselowej. Oznacza to, że prędkości znacznie przekraczające tę granicę są szacowane z dużą dokładnością, natomiast dokładność ta maleje, gdy wartości prędkości zbliżają się do tej granicy lub spadają poniżej niej.
Ograniczenie to oznacza również, że zwiększenie częstotliwości próbkowania klatek (zmniejszenie Δt) z jednej strony pozwala na łatwiejszą ocenę dużych prędkości i ich gwałtownych zmian (np. podczas faz izowolumetrycznych). Z drugiej strony, zwiększenie częstotliwości klatek (zmniejszenie Δt) podnosi tę granicę i wiąże się z obniżeniem dokładności oceny niższych prędkości, o ile nie towarzyszy temu podobny wzrost rozdzielczości przestrzennej (zmniejszenie Δx).
Przygotowanie obrazu fantoma
Przygotowano serię sztucznych, wygenerowanych komputerowo pętli, aby umożliwić testowanie procedury analizy obrazu w prostych i w pełni kontrolowanych warunkach. W tym celu przygotowano fantom w projekcji krótkiej osi idealnej lewej komory zgodnie z poniższym opisem.
Granice wsierdzia i osierdzia są reprezentowane przez dwa koncentryczne okręgi o promieniach odpowiednio R0(t) i R1(t). Obraz przygotowano poprzez nadanie pierścieniowi, który reprezentuje tkankę pomiędzy obiema granicami, jednolitego szarego koloru na czarnym tle. Następnie zastosowano liniowy filtr top-hat o rozmiarze 8x8, aby uniknąć niefizycznych nieciągłości.
Ruch epikardium wyrażono w [mm] jako R0(t)=10+5cos(2πt/T), gdzie T to okres bicia serca przyjęty jako T=1s. Teoretyczna kinematyka endokardium jest stała wzdłuż krawędzi i zależy wyłącznie od czasu; prędkość jest jedynie radialna i dana wzorem V0(t)= dR0/dt=-π sin(2πt/T), w [cm/s]. Odkształcenie procentowe, obliczone w stosunku do długości krawędzi w czasie zero, wynosi St0(t)=100x(R0(t)-R0(0))/R0(0)=100(cos(2πt/T) -1)/3, a szybkość odkształcenia wynika z równania (1) jako SR0(t)=10 V0/R0, w [s-1]. Przyjęto, że epikardium porusza się zgodnie ze stałą grubością, R1(t)= R0(t)+5mm, lub pozostaje w spoczynku R1(t)= R0(0)+5mm.
Każdy obraz jest kwadratem o rozmiarze 48mm, wycentrowanym na pierścieniu tkanki, i posiada rozdzielczość NxN. Przykładowe obrazy przedstawiono na rysunku 4, płyty a i b; profile czasowe odkształcenia i szybkości odkształcenia pokazano na rysunku 4, płyty c i d. Pętle są przygotowywane poprzez zmianę rozdzielczości N, częstotliwości klatek FR oraz rodzaju ruchu wsierdziowego.
Metoda śledzenia wsiercdzia jest stosowana do takich obrazów poprzez wyznaczenie na pierwszej klatce liczby Np punktów równomiernie rozmieszczonych wzdłuż okrągłego wsierdzia.
3. Reprezentatywne wyniki
Badanie na fantomie
W niniejszej pracy przeanalizowano zastosowanie metody analizy obrazu do obrazów fantomów generowanych komputerowo. Globalną miarę końcowego błędu obliczono za pomocą pierwiastka średniokwadratowego procentowej różnicy. Pierwiastek średniokwadratowy, błąd średni oraz maksymalny błąd odkształcenia wsierdzia zdefiniowano jako

Równanie 2
gdzie St0(t) jest wartością dokładną, St(t) jest wartością obliczoną poprzez analizę obrazu, a sumowanie obejmuje wszystkie klatki NF=FRxT. Tą samą definicję zastosowano dla promienia, prędkości i szybkości odkształcenia. Śledzenie jest w zasadzie niezależne od pozycji wzdłuż wsierdzia, a różnice między poszczególnymi punktami są znacznie mniejsze niż 1%.
Wyniki podsumowano w Tabeli I dla 15 fantomów o różnej rozdzielczości przestrzennej, częstotliwości klatek oraz typie ruchu granicy epikardium; przedstawiono również efekt zmiany liczby punktów wykorzystanych do śledzenia granicy endokardium.
Błędy są we wszystkich przypadkach bardzo małe dla wartości całkowych (promień i odkształcenie) i nieco większe dla wartości różniczkowych (prędkość i szybkość odkształcenia), które są powiązane z pochodną tych pierwszych. Było to oczekiwane, ponieważ operator pochodnej potęguje błędy. Jakość wyników pogarsza się wraz ze zmniejszeniem rozdzielczości; w rzeczywistości dokładność jest powiązana z rozmiarem piksela, który reprezentuje (w szerokim sensie) minimalne przemieszczenie możliwe do odczytania z jednej klatki do drugiej. Rozdzielczość czasowa nie wpływa znacząco na wyniki, dopóki liczba klatek na sekundę jest wystarczająca; przy bardzo wysokiej częstotliwości klatek wyniki nie ulegają poprawie, ponieważ przemieszczenia między poszczególnymi klatkami stają się mniejsze niż rozmiar piksela. Pokazuje to, że zwiększenie liczby klatek na sekundę ma niewielką lub żadną użyteczność, gdy nie towarzyszy mu zwiększenie rozdzielczości przestrzennej.
Jednak rozważony tutaj prosty ruch sinusoidalny nie wymaga ekstremalnej rozdzielczości czasowej. Podobnie, zastosowanie zaledwie 8 punktów jest wystarczające, aby odwzorować prosty, kolisty kształt wsierdzia. Wyniki dotyczące wsierdzia nie są w znacznym stopniu zależne od rodzaju ruchu, jakiego poddawany jest epikardium. Potwierdziliśmy również, że wyniki nie są istotnie wpływowane przez przyjęte filtrowanie obrazu.
Wizualna prezentacja wyników przedstawiona jest na Rysunku 4, gdzie obliczona granica endokardium w dwóch chwilach została nałożona na obrazy fantoma (tablice a i b). Odkształcenie oraz szybkość odkształcania są przedstawione na (tablicach c i d) dla przypadku nr 1 oraz przypadku nr 8 o niższej rozdzielczości. Odkształcenie i szybkość odkształcania w przypadku nr 1 (kwadraty) wykazują doskonałą zgodność z wartością teoretyczną, przy średnim błędzie wynoszącym odpowiednio 0,6% i 3%. Zgodność jest nieco mniejsza w przypadku nr 8, w którym rozdzielczość obrazu jest dwukrotnie mniejsza, a błędy dla odkształcenia i szybkości odkształcania wynoszą odpowiednio 0,9% i 4,5%.
Walidacja kliniczna 1.
Porównano obwodowy odkształcenie mięśnia sercowego (εcc) w całym przekroju środkowej części lewej komory (LV) przy użyciu technik Harmonic Phase Imaging (HARP) oraz FT u 191 pacjentów z dystrofią mięśniową Duchenne'a, podzielonych według wieku i stopnia niewydolności serca, oraz u 42 osób kontrolnych dopasowanych pod względem wieku. Przeprowadzono retrospektywną analizę off-line dopasowanych przekrojów z obrazowania tagowanego oraz SSFP. Dla całej populacji badanej (n=233) średnie wartości FT εcc (-13,3 ± 3,8 %) były silnie skorelowane z HARP εcc (-13,6±3,4 %), przy współczynniku korelacji Pearsona wynoszącym 0,899. Średnie wartości εcc u pacjentów z DMD wyznaczone za pomocą HARP (-12,52 ± 2,69 %) i FT (-12,16 ± 3,12 %) nie różniły się istotnie (p=NS). Podobnie, średnie wartości εcc osób kontrolnych wyznaczone za pomocą HARP (-18,85 ± 1,86) i FT (-18,81 ± 1,83) nie różniły się istotnie (p=NS). Wyciągnięto wniosek, że ocena εcc oparta na metodzie FT wykazuje wysoką korelację z εcc pochodzącym z obrazów tagowanych w dużej populacji pacjentów z DMD z szerokim zakresem niewydolności serca.
| Fantom nr | Częstotliwość klatek | Rozdzielczość. N | Ruchy osierdzia | Np | εR | εRmaksymalny | εV | εVmaksymalny | εSR | εSRmaksymalne | εSt | εStmaksymalne |
| 1 | 32 | 401 | z endo | 16 | 0.57 | 1.10 | 3.98 | 8.26 | 2.95 | 6.55 | 0.59 | 1.34 |
| 2 | 32 | 401 | brak ruchu | 16 | 0.27 | 0.60 | 3.62 | 8.53 | 2.90 | 5.99 | 0.75 | 1.29 |
| 3 | 64 | 401 | z endo | 16 | 0.35 | 0.62 | 1.75 | 4.92 | 3.08 | 10.47 | 0.34 | 0.86 |
| 4 | 64 | 401 | brak ruchu | 16 | 0.39 | 0.73 | 1.78 | 5.48 | 2.38 | 6.46 | 0.32 | 0.63 |
| 5 | 128 | 401 | z endo | 16 | 0.22 | 0.40 | 1.57 | 3.66 | 2.73 | 8.39 | 0.35 | 0.89 |
| 6 | 16 | 401 | z endo | 16 | 0.41 | 0.92 | 12.10 | 18.27 | 7.12 | 15.17 | 0.88 | 1.93 |
| 7 | 32 | 201 | z endo | 16 | 0.43 | 0.82 | 3.84 | 7.03 | 4.21 | 7.79 | 0.83 | 1.52 |
| 8 | 32 | 201 | brak ruchu | 16 | 0.49 | 0.90 | 4.00 | 8.54 | 4.54 | 10.84 | 0.93 | 1.64 |
| 9 | 32 | 101 | z endo | 16 | 2.58 | 3.87 | 5.70 | 17.44 | 9.22 | 18.27 | 4.44 | 6.77 |
| 10 | 32 | 401 | z endo | 64 | 0.32 | 0.73 | 3.98 | 8.09 | 3.11 | 8.86 | 0.48 | 1.10 |
| 11 | 32 | 401 | brak ruchu | 64 | 0.33 | 0.53 | 3.56 | 7.83 | 2.78 | 7.44 | 0.70 | 0.98 |
| 12 | 32 | 401 | z endo | 32 | 0.44 | 1.00 | 3.99 | 8.16 | 2.99 | 6.82 | 0.63 | 1.49 |
| 13 | 32 | 401 | z endo | 8 | 0.20 | 0.41 | 3.64 | 6.84 | 3.09 | 9.10 | 0.43 | 0.74 |
| 14 | 32 | 401 | brak ruchu | 8 | 0.22 | 0.41 | 3.36 | 7.42 | 2.76 | 5.93 | 0.58 | 1.01 |
| 15* | 32 | 401 | z endo | 16 | 1.35 | 2.42 | 5.54 | 9.68 | 7.02 | 16.58 | 2.62 | 4.46 |
Tabela 1. Analiza fantomowa śledzenia brzegów wsierdziowych: pierwiastka średniokwadratowego oraz maksymalne błędy procentowe [%] obliczono dla głównych wielkości w odniesieniu do różnych parametrów fantomów. Parametry zaznaczone pogrubieniem wskazują zmiany względem Fantomu nr 1. Rozważono zależność od częstotliwości klatek, rozdzielczości oraz liczby śledzonych punktów. Wpływ rodzaju ruchu osierdzia rozważono dla dwóch przypadków granicznych: gdy brzeg osierdzia się nie porusza (brak ruchu) lub porusza się wraz z wsierdziem (brak pogrubienia). Ostatni fantom (*) został zbudowany bez filtrowania podstawowych obrazów skokowych z gwałtowną zmianą jasności w jednym pikselu. Błędy przekraczające 10% zaznaczono pogrubieniem.

Rysunek 1. Obraz cMR lewej komory w przekroju podłużnym (lewy obraz) i w przekroju poprzecznym (prawy obraz), z naniesionym obrysem granicy wsierdzia.

Rycina 2. Reprezentacja czasoprzestrzenna sekwencji obrazów, w której przestrzeń przebiega wzdłuż przekroju transmuralnego. Przekrój transmuralny wykonano zgodnie z punktem początkowym na rycinie 3. Przedstawiono ewolucję czasową automatycznie śledzonego punktu początkowego.

Rysunek 3. Obraz lewej komory w widoku osi podłużnej, z cięciami transmuralnymi oraz cięciami równoległymi do zakrzywionej krawędzi.

Rysunek 4. Badanie z użyciem fantomów. Dwa obrazy (przypadek nr 2) przy maksymalnym rozszerzeniu (płyta a) i skurczu (płyta b), obliczone punkty granicy wsiekoliczne są nałożone na siebie. Odkształcenie (płyta c) oraz szybkość odkształcenia (płyta d) obliczone dla dwóch różnych fantomów (przypadki nr 1 i nr 8) przedstawiono w porównaniu z wartościami efektywnymi.

Rycina 5. Przykłady globalnego odkształcenia obwodowego (krzywa czarna) i odkształcenia segmentowego (krzywe kolorowe) u zdrowych pacjentów (a). Przykład globalnego odkształcenia obwodowego (krzywa czarna) i segmentowego odkształcenia obwodowego (krzywe kolorowe) u pacjentów z obniżoną funkcją lewej komory i blokiem lewej odnogi pęczka Hisa (b). Zwrócić uwagę na różny czas wystąpienia szczytowego odkształcenia obwodowego, co wskazuje na obecność dyssynchronii lewej komory