USG/Doppler

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.

Zdjęcie fantomu przepływowego.

Schemat obrazujący przekrój fantomu odpowiadający w przybliżeniu płaszczyźnie obrazowania.

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).

Widmo sygnału: a) przed demodulacją; b) po demodulacji; c)po demodulacji i filtrowaniu.

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:

1. 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?
2. 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)?
3. 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

Po demodulacji otrzymaliśmy sygnał o pasmowym widmie skoncentrowanym w pobliżu zera. Następnie chcemy estymować średnią częstotliwość (tj. częstotliwość średnią wążoną energiami dla poszczególnych częstotliwości) tego sygnału; przesunięcie częstotliwości średniej od zera świadczyć będzie o tym, że rozpraszacze poruszają się w kierunku od lub do głowicy. Ś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, 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.
Praca Miller

Zadanie

Następnie przygotowujemy skrypt realizujący estymator autokorelacyjny. Dotychczas traktowaliśmy nasze obrazy jako dwuwymiarowy sygnał. Na potrzeby estymacji mapy prędkości uwzględnić musimy jeszcze wymiar czasowy. Dla danego punktu (ustalona głębokość i szerokość) estymator szacuje średnią częstotliwość sygnału [math]s(t)[/math] w oknie czasowym długości [math]N[/math] jako

[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]

gdzie [math]T_{PRF}[/math] - czas między kolejnymi strzałami. Tak otrzymaną estymatę średniego przesunięcia [math]\Delta f=f_0-f_s[/math] dopplerowskiego możemy wykorzystać do obliczenia średniej prędkości przy pomocy szkolnego wzoru 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ą. 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 (tzw. objętość pomiarowa).

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)

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.
Proszę 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)