USG/Doppler

Metoda dopplerowska

Do dyspozycji mamy zestaw danych RF z 31 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach: UWAGA: te parametry niekoniecznie są zgodne z Prawdą.

f0 = 5.5e6 # Częstotliwość nadawcza [Hz]
fs = 9e6 # Częstotliwość próbkowania [Hz]
pitch = 0.21e-3 # odległość pomiędzy środkami kolejnych przetworników [m]

NT = 192 # Liczba przetworników w pełnej aperturze
Ntr = 192 # subapertura nadawcza

T_PRF = ????

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). 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 musimy zrekonstruować 31 obrazów (z pełną dynamiką jasności, przed liczeniem obwiedni) - proszę wykorzystać w tym celu funkcje przygotowane na poprzednich zajęciach.

Następnie 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

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ść 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)