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

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)