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 sa zgodne z Prawdą.

f0=5.5e6 # Częstotliwość nadawcza przetworników [Hz]
fs=9e6 # Częstotliwość próbkowania [Hz]
pitch = 0.00021 # Deklarowana odległość między środkami przetworników nadawczo-odbiorczych

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 wywołuje przepływ o stałym natężeniu płynu krwiopodobnego znajdującego się w rurkach. Będziemy starali się wykorzystać metodę dopplerowską do stworzenia mapy obrazujacej zwrot i względne natężenie przepływu. Średnie odchylenie dopplerowskie szacować będziemy przy użyciu estymatora autokorelacyjnego.

Demodulacja danych

Estymator autokorelacyjny stosowany jest do danych zdemodulowanych (IQ). Demodulację możemy przeprowadzić zarówno na danych surowych, 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.55 MHz

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 szkolnego wzoru 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 uwadzę, ż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 niepunktowego obszaru pomiarowego.

Prezentacja Color

Po wyestymowaniu ś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 BMode w celu uzyskania obrazu Color. Mając tablicę z danymi BMode 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ścia 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)