USG/Doppler: Różnice pomiędzy wersjami
Linia 29: | Linia 29: | ||
==Estymator autokorelacyjny== | ==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.<br> | ||
+ | Średnia częstość kołowa <math>\bar{\omega}</math> widma dopplerowskiego <math>P(\omega)</math> może być zdefiniowana jako <br><br> | ||
+ | <math>\bar{\omega}=\frac{\int^{\infty}_{-\infty}\omega P(\omega)\mbox{d}\omega}{P(\omega)\mbox{d}\omega}</math> | ||
+ | <br><br> | ||
+ | Jednocześnie, przy założeniu słabej stacjonarności sygnału dopplerowskiego, na mocy [[Twierdzenie_Wienera-Chinczyna|twierdzenia Wienera-Chinczyna]], funkcję autokorelacji <math>R(\tau)</math> tego sygnału możemy wyrazić jako <br><br> | ||
+ | <math>R(\tau)=\int^{\infty}_{-\infty} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math> | ||
+ | <br><br> | ||
+ | Różniczkując powyższe dostajemy <br><br> | ||
+ | <math>\dot{R}(\tau)=j\int^{\infty}_{-\infty} \omega P(\omega)e^{j\omega\tau}\mbox{d}\omega</math><br><br> | ||
+ | oraz <br><br> | ||
+ | <math>\ddot{R}(\tau)=-\int^{\infty}_{-\infty} \omega^{2} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math><br><br> | ||
+ | Stąd, częstotliwość średnia możemy być również przedstawiona jako <br><br> | ||
+ | <math>\bar{\omega}=-j\frac{\dot{R}(0)}{R(0)}</math><br><br> | ||
+ | W celu uproszczenia obliczeń przyjmuje się często następujące uproszczenie | ||
+ | <br><br><math>R(\tau)= |R(\tau)|e^{j\phi(\tau)}=A(\tau)e^{j\phi(\tau)} </math><br><br> | ||
+ | gdzie <math>A(\tau)</math> jest funkcją parzystą i <math>e^{j\phi(\tau)}</math> jest funkcją nieparzystą. Wtedy | ||
+ | <br><br><math>\dot{R}(\tau)=(\dot{A}(\tau)+jA(\tau)\dot{\phi}(\tau)e^{j\phi(\tau)} </math><br><br> | ||
+ | i stąd | ||
+ | <br><br><math>\dot{R}(0)=jA(0)\dot{\phi}(0) </math><br><br> | ||
+ | oraz | ||
+ | <br><br><math>R(0)=A(0) </math><br><br> | ||
+ | Korzystając ze wcześniejszych równań dostajemy | ||
+ | <br><br><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><br><br> | ||
+ | gdzie <math>T_{PRF}</math> - czas między kolejnymi strzałami. <br> | ||
+ | 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 <br><br> | 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 <br><br> | ||
<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><br><br> | <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><br><br> |
Wersja z 10:01, 30 cze 2016
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
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)