USG/Doppler: Różnice pomiędzy wersjami
m |
m |
||
| Linia 38: | Linia 38: | ||
<math>\Delta f=\frac{2f_{0}v\mbox{cos}\theta}{c}</math> | <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ą. | gdzie <math>v</math> średnia prędkość w punkcie pomiarowym; <math>\theta</math> kąt między kierunkiem przepływu a wiązką nadawczą. | ||
| − | Przesunięcie dopplerowskie jest różnicą między średnią częstotliwością <math>f_0</math> sygnału nadanego a średnią częstotliwością <math>f_s</math> sygnału odebranego (pochodzącego od rozpraszaczy przemieszczających się w kierunku do lub od głowicy). Przypomnijmy, że zdemodulowany sygnał jest sygnałem o widmie pasmowym położonym w pobliżu częstości zerowej. Pozwala nam to przyjąć, że średnie przesunięcie częstotliwości odpowiada średniej częstotliwości tego sygnału | + | Przesunięcie dopplerowskie jest różnicą między średnią częstotliwością <math>f_0</math> sygnału nadanego a średnią częstotliwością <math>f_s</math> sygnału odebranego (pochodzącego od rozpraszaczy przemieszczających się w kierunku do lub od głowicy). Przypomnijmy, że zdemodulowany sygnał jest sygnałem o widmie pasmowym położonym w pobliżu częstości zerowej. Pozwala nam to przyjąć, że średnie przesunięcie częstotliwości odpowiada średniej częstotliwości tego sygnału <math>\Delta f = f_s-f_0 = f_s</math>. |
| − | <math>\Delta f = f_s-f_0 = f_s</math>. | ||
Ś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<ref>Kasai, Chihiro, et al. "Real-time two-dimensional blood flow imaging using an autocorrelation technique." IEEE Trans. Sonics Ultrason 32.3 (1985): 458-464.</ref>, oparty na czasowej reprezentacji sygnału. | Ś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<ref>Kasai, Chihiro, et al. "Real-time two-dimensional blood flow imaging using an autocorrelation technique." IEEE Trans. Sonics Ultrason 32.3 (1985): 458-464.</ref>, 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 | Ś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> | <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 [[Twierdzenie_Wienera-Chinczyna|twierdzenia Wienera-Chinczyna]], funkcję autokorelacji <math>R(\tau)</math> tego sygnału możemy wyrazić jako: | 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: | ||
<math> R(\tau)=\int^{\infty}_{-\infty} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math> | <math> R(\tau)=\int^{\infty}_{-\infty} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math> | ||
Wersja z 11:22, 3 lip 2016
Spis treści
Metoda dopplerowska
Do dyspozycji mamy zestaw zrekonstruowanych danych RF (link) 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]
PRF = 1000 # Pulse Repetition Frequency [Hz]
T_PRF = 1./PRF # 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 na zdemodulowanym sygnale kwadraturowym I/Q (ang. In-Phase/Quadrature). 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 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.
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 sygnał należy przefiltrować dolnoprzepustowo w zakresie pasma podstawowego (tzw. baseband). 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
Dotychczas traktowaliśmy rekonstruowane obrazy jako sygnały dwuwymiarowe. Teraz, do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie, konieczna będzie analiza wielu następujących po sobie obrazów, które zostały zebrane ze stałym interwałem czasowym (tzw. PRF - Pulse Repetition Frequency). W poniższych rozważaniach pomijamy wymiary przestrzenne - całe rozumowania przeprowadzone są dla sygnału [math]s(i)[/math] na ustalonej głębokości i szerokości. 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 (energia danego piksela zawiera informację o obiektach z pewnego obszaru). Średnią prędkość będziemy mogli obliczyć korzystając ze średniego przesunięcia dopplerowskiego [math]\Delta f[/math] w oparciu o szkolny wzór 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ą. Przesunięcie dopplerowskie jest różnicą między średnią częstotliwością [math]f_0[/math] sygnału nadanego a średnią częstotliwością [math]f_s[/math] sygnału odebranego (pochodzącego od rozpraszaczy przemieszczających się w kierunku do lub od głowicy). Przypomnijmy, że zdemodulowany sygnał jest sygnałem o widmie pasmowym położonym w pobliżu częstości zerowej. Pozwala nam to przyjąć, że średnie przesunięcie częstotliwości odpowiada średniej częstotliwości tego sygnału [math]\Delta f = f_s-f_0 = f_s[/math]. Ś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[1], 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 (odwrotność PRF).
Estymator Millera-Rochwargera
Wartość [math]R(T_{PRF})[/math] estymować możemy np. w oparciu o estymator Millera-Rochwargera [2]. Jest to estymator maksymalnej wiarygodności (maximum likehood estimator). Estymator ten jest skonstruowany dla par obserwacji zespolonego procesu [math]s[/math] próbkowanych w równych odstępach czasu (w naszym przypadku [math]T_{PRF}[/math]). W naszym przypadku jako pary obserwacji traktuje się pary próbek z kolejnych par następujących po sobie pomiarów (oddzielonych wartością [math]T_{PRF}[/math]). Otrzymujemy wtedy oszacowanie wartości autokorelacji w oknie czasowym długości [math]N-1[/math]
[math]\bar{\omega}=\frac{1}{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]
Pamiętając zależność między średnią częstotliwością [math]f_s[/math] a średnią częstością kołową [math]fs=(1/2\pi)\bar{\omega}[/math] dostajemy
[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]
Zadanie
Proszę przygotować funkcję estymującą częstotliwość średnią w oparciu o powyższy estymator dla jednowymiarowego zespolonego sygnału.
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)
Zadanie
Zaimplementować funkcję otrzymującą na wejściu trójwymiarową tablicę zawierającą obrazy RF z [math]N[/math] kolejnych chwil pomiarowych i zwracającą obrazy z nałożoną na nią mapą prędkości obliczonych przy użyciu estymatora autokorelacyjnego.
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.
Zadanie
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)
Zadanie
Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji)
- ↑ Kasai, Chihiro, et al. "Real-time two-dimensional blood flow imaging using an autocorrelation technique." IEEE Trans. Sonics Ultrason 32.3 (1985): 458-464.
- ↑ Miller, Kenneth, and M. Rochwarger. "A covariance approach to spectral moment estimation." IEEE Transactions on Information Theory 18.5 (1972): 588-596.
