USG/Doppler
Spis treści
Metoda dopplerowska
Do dyspozycji mamy zestaw zrekonstruowanych danych RF z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:
f0 = 4e6 # Częstotliwość nadawcza [Hz]
z_step = 1.8970189701897084e-05 # Odległość między kolejnymi punktami w głębokości na siatce rekonstruowanego obrazu [m]
T_PRF = 1./1000 #czas między kolejnymi nadaniami [s]
Obrazowany jest przekrój fantomu przepływowego złożonego z rurek umieszczonych w materiale tkankopodobnym.
Pompa wymusza jednostajny przepływ płynu krwiopodobnego znajdującego się w rurkach. Będziemy starali się wykorzystać metodę dopplerowską do stworzenia mapy obrazującej zwrot i prędkość przepływu. Średnie odchylenie dopplerowskie szacować będziemy przy użyciu estymatora autokorelacyjnego.
Demodulacja sygnału RF
Estymator autokorelacyjny stosowany jest zdemodulowanym sygnale (I/Q). Sygnał RF rejestrowany bezpośrednio z przetworników odbiorczych jest sygnałem rzeczywistym (zerowa część urojona) o pasmowej charakterystyce widmowej. Informacja istotna z naszego punktu widzenia (czyli np. nie szum elektroniczny) położona jest w paśmie, którego środek znajduje się częstotliwości nadawczej; szerokość tego pasma zależna jest od fizycznych właściwości głowicy ultradźwiękowej i nazywana jest pasmem przenoszenia głowicy. Demodulacja jest operacją pozwalającą na przekształcenie takiego sygnału w sygnał zespolony o paśmie położonym wokół częstotliwości zerowej. Taki zabieg pozwala na zmniejszenie częstotliwości próbkowania - potencjalnie poniżej częstotliwości Nyquista określonej dla składowych sygnału niezdemodulowanego. Demodulacja jest operacją odwracalną, dlatego można myśleć o niej jako o bezstratnej kompresji sygnału (przy założeniu, że istotnie informacja poza pasmem jest pomijalna jako szum).
Demodulację możemy przeprowadzić zarówno na danych surowych RF, jak i na danych po rekonstrukcji. W naszym wypadku zastosujemy to drugie rozwiązanie.
Na początku stworzymy funkcję dokonującą demodulacji każdego z obrazów, tj. dla każdego próbki o współrzędnych [math](x,y)[/math]
[math]IQ(x,y)=RF(x,y)\cdot e^{-i\cdot2\pi f_{0}t} [/math]
gdzie [math]IQ[/math] - sygnał po demodulacji; [math]RF[/math] - sygnał przed demodulacją; [math]t[/math] - czas odpowiadający momentowi akwizycji próbki z danej głębokości; możemy przyjąć uproszczone założenie, że:
[math]t=2y/c [/math]; jak widać, pierwszy wymiar (szerokość) jest w naszej procedurze nieistotny.
Po demodulacji dane należy przefiltrować. W naszym wypadku wystarczający powinien być filtr dolnoprzepustowy o częstotliwości odcięcia równej 5.55MHz
Do wykonania:
- Zbadać widmo amplitudowe sygnału przed i po demodulacji (po filtrowaniu). Jak zmieniło się widmo? Gdzie jest położona średnia widma? Czy rozkłady dla ujemnych i dodatnich częstotliwości są swoimi odbiciami?
- Porównać widmo amplitudowe sygnału zdemodulowanego przed i po filtracji. Jakich składowych w częstości się pozbyliśmy (nie licząc szumu)?
- Zmienić w demodulacji wartość częstotliwości nośnej (np. zmniejszyć o połowę) i ponownie porównać widmo przed i po demodulacji.
Estymator autokorelacyjny
Rys teoretyczny
Dotychczas traktowaliśmy nasze obrazy jako dwuwymiarowy sygnał. Będziemy jednak dążyć do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie obszaru pomiarowego. Wymagać to będzie od nas uwzględnienia wymiaru czasowego. W poniższych rozważaniach pomijamy wymiary przestrzenne - całe rozumowania przeprowadzone są dla sygnału [math]s(i)[/math] na ustalonej głębokości i szerokości. Należy mieć na uwadze, że w praktyce ciężko mówić o estymacji prędkości w punkcie, ponieważ w naszej procedurze estymacji uwzględniamy efektywnie informacje z pewnego obszaru pomiarowego (energia danego piksela zawiera informację o obiektach z pewnego niepunktowego obszaru).
Średnią prędkość będziemy mogli obliczyć korzystając ze średniego przesunięcia dopplerowskiego [math]\Delta f[/math] w oparciu o szkolny wzór na częstotliwość Dopplera:
[math]\Delta f=\frac{2f_{0}v\mbox{cos}\theta}{c}[/math]
gdzie [math]v[/math] średnia prędkość w punkcie pomiarowym; [math]\theta[/math] kąt między kierunkiem przepływu a wiązką nadawczą.
Przesunięcie dopplerowskie jest różnicą między średnią częstotliwością [math]f_0[/math] sygnału nadanego a średnią częstotliwością [math]f_s[/math] sygnału odebranego (pochodzącego od rozpraszaczy przemieszczających się w kierunku do lub od głowicy). Przypomnijmy, że zdemodulowany sygnał jest sygnałem o widmie pasmowym położonym w pobliżu częstości zerowej. Pozwala nam to przyjąć, że średnie przesunięcie częstotliwości odpowiada średniej częstotliwości tego sygnału
[math]\Delta f = f_s-f_0 = f_s[/math]
Średnią częstotliwość możemy estymować na wiele sposób - np. licząc średnią ważoną z widma otrzymanego przez dyskretną transformację Fouriera. My zastosujemy do tego estymator autokorelacyjny[1], oparty na czasowej reprezentacji sygnału.
Średnia częstość kołowa [math]\bar{\omega}[/math] widma dopplerowskiego [math]P(\omega)[/math] może być zdefiniowana jako
[math]\bar{\omega}=\frac{\int^{\infty}_{-\infty}\omega P(\omega)\mbox{d}\omega}{P(\omega)\mbox{d}\omega}[/math]
Jednocześnie, przy założeniu słabej stacjonarności sygnału dopplerowskiego, na mocy twierdzenia Wienera-Chinczyna, funkcję autokorelacji [math]R(\tau)[/math] tego sygnału możemy wyrazić jako
[math]R(\tau)=\int^{\infty}_{-\infty} P(\omega)e^{j\omega\tau}\mbox{d}\omega[/math]
Różniczkując powyższe dostajemy
[math]\dot{R}(\tau)=j\int^{\infty}_{-\infty} \omega P(\omega)e^{j\omega\tau}\mbox{d}\omega[/math]
oraz
[math]\ddot{R}(\tau)=-\int^{\infty}_{-\infty} \omega^{2} P(\omega)e^{j\omega\tau}\mbox{d}\omega[/math]
Stąd, częstotliwość średnia możemy być również przedstawiona jako
[math]\bar{\omega}=-j\frac{\dot{R}(0)}{R(0)}[/math]
W celu uproszczenia obliczeń przyjmuje się często następujące uproszczenie
[math]R(\tau)= |R(\tau)|e^{j\phi(\tau)}=A(\tau)e^{j\phi(\tau)} [/math]
gdzie [math]A(\tau)[/math] jest funkcją parzystą i [math]e^{j\phi(\tau)}[/math] jest funkcją nieparzystą. Wtedy
[math]\dot{R}(\tau)=(\dot{A}(\tau)+jA(\tau)\dot{\phi}(\tau)e^{j\phi(\tau)} [/math]
i stąd
[math]\dot{R}(0)=jA(0)\dot{\phi}(0) [/math]
oraz
[math]R(0)=A(0) [/math]
Korzystając ze wcześniejszych równań dostajemy
[math]\bar{\omega}=\dot{\phi}(0)\approx \frac{\phi(T_{PRF})-\phi(0)}{T_{PRF}}=\frac{\phi(T_{PRF})}{T_{PRF}}=\frac{1}{T_{PRF}}\frac{\mbox{Im}(R(T_{PRF}))}{\mbox{Re}(R(T_{PRF}))} [/math]
gdzie [math]T_{PRF}[/math] - czas między kolejnymi strzałami.
Estymator Millera-Rochwargera
Wartość [math]R(T_{PRF})[/math] estymować możemy np. w oparciu o estymator Millera-Rochwargera [2]. Jest to estymator maksymalnej wiarygodności (maximum likehood estimator). Estymator ten jest skonstruowany dla par obserwacji zespolonego procesu [math]s[/math] próbkowanych w równych odstępach czasu (w naszym przypadku [math]T_{PRF}[/math]). W naszym przypadku jako pary obserwacji traktuje się pary próbek z kolejnych par następujących po sobie pomiarów (oddzielonych wartością [math]T_{PRF}[/math]). Otrzymujemy wtedy oszacowanie wartości autokorelacji w oknie czasowym długości [math]N-1[/math]
[math]\bar{\omega}=\frac{1}{T_{PRF}}\mbox{arctan}\frac{Im(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)}}{Re(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})} [/math]
Pamiętając zależność między średnią częstotliwością [math]f_s[/math] a średnią częstością kołową [math]fs=(1/2\pi)\bar{\omega}[/math] dostajemy
[math]f_s=\frac{1}{2\pi T_{PRF}}\mbox{arctan}\frac{Im(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)}}{Re(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})} [/math]
Zadanie
Proszę przygotować funkcję estymującą częstotliwość średnią w oparciu o powyższy estymator dla jednowymiarowego zespolonego sygnału.
Prezentacja Kolor
Po estymacji średnich prędkości w obszarach odpowiadających punktom na siatce użytej do rekonstrukcji obrazu, możemy nałożyć taką mapę na obraz B-mode w celu uzyskania obrazu "Kolor". Mając tablicę z danymi B-mode oraz tablicę prędkości, możemy otrzymać połączony obraz za pomocą poniższego skryptu:
#BMode - tablica zrekonstruowanych danych po przefiltrowaniu, obwiedni itp.
#flow - tablica przepływów
Frame = Image.fromarray(np.uint8(cm.bone(BMode)*255))
flowMask = np.copy(flow)
flowMask = np.abs(flowMask)
flowMask = flowMask/(np.max(flowMask))*255
flow = flow+np.abs(np.min(flow))
flow = flow/np.max(flow)
flow = Image.fromarray(np.uint8(cm.jet(flow)*255))
flowMask = Image.fromarray(np.uint8(flowMask), 'L')
flow.putalpha(flowMask)
Frame.paste(flow, (0,0), flow)
Zadanie
Zaimplementować funkcję otrzymującą na wejściu trójwymiarową tablicę zawierającą obrazy RF z [math]N[/math] kolejnych chwil pomiarowych i zwracającą obrazy z nałożoną na nią mapą prędkości obliczonych przy użyciu estymatora autokorelacyjnego.
Filtracja obrazu
Mapa częstości estymowana metodą autokorelacyjną jest dość wrażliwa na błędy estymacji powodowane m.in. obecność w sygnale informacji z dużego obszaru pomiarowego. Stosować można kilka metod mających na celu poprawę wynikowego obrazu - progowanie (zignorowanie względnie małych prędkości), rozpoznawanie ruchu i filtrowanie.
Rozpoznawanie ruchu
Jedną z pierwszy obserwacji jakiej można dokonać, to zauważenie, że nasza mapa jest niezerowa w punktach w których nie spodziewamy się żadnego ruchu. Pewnym rozwiązaniem tego problemu może być zastosowanie prostego kryterium rozpoznawania ruchu. Możemy przyjąć, że sygnał pochodzący od struktur pozostających w spoczynku powinien być stały w czasie. W praktyce sygnały takie różnić będą się głównie o składową szumu elektronicznego. Jednocześnie, sygnał pochodzący od ruchomych struktur będzie charakteryzować się względnie dużą zmiennością w czasie.
Zadanie
Zaimplementować procedurę zerującą estymatę prędkości w punkcie, jeśli średnia różnica między wartościami sygnału w tym punkcie w kolejnych chwilach czasu jest względnie mała (próg proszę dobrać eksperymentalnie - zaczynając np. od progu 10% średniej różnicy w obrazie).
Filtr medianowy
Dobrym rozwiązaniem w tego typu przypadkach jest zastosowanie filtru medianowego. Filtr medianowy przekształca wartość w punkcie na medianę wartości sygnału w ustalonej liczbie sąsiednich próbek (w obu wymiarach):
[math]s_{med}(x,y) = \sum^{K-1}_{i=0}\sum^{K-1}_{j=0}s(x+i,y+j) [/math]
gdzie [math]K[/math] - rozmiar okna filtru.
Filtr taki usuwa skrajne wartości w dużo większym stopniu niż np. filtr oparty o średnią arytmetyczną. W bibliotece scipy istnieje gotowa funkcja filtrująca tablicę filtrem medianowym:
scipy.signal.medfilt2d(array, K)
- ↑ Kasai, Chihiro, et al. "Real-time two-dimensional blood flow imaging using an autocorrelation technique." IEEE Trans. Sonics Ultrason 32.3 (1985): 458-464.
- ↑ Miller, Kenneth, and M. Rochwarger. "A covariance approach to spectral moment estimation." IEEE Transactions on Information Theory 18.5 (1972): 588-596.
