USG/Doppler: Różnice pomiędzy wersjami

Z Brain-wiki
Linia 30: Linia 30:
 
==Estymator autokorelacyjny==
 
==Estymator autokorelacyjny==
 
===Rys teoretyczny===
 
===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>
+
Dotychczas traktowaliśmy nasze obrazy jako dwuwymiarowy sygnał. Będziemy jednak dążyć do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie obszaru pomiarowego. Wymagać to będzie od nas uwzględnienia wymiaru czasowego. 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 niepunktowego obszaru).<br><br>
 +
Ś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:<br><br>
 +
<math>\Delta f=\frac{2f_{0}v\mbox{cos}\theta}{c}</math><br><br>
 +
gdzie <math>v</math> średnia prędkość w punkcie pomiarowym; <math>\theta</math> kąt między kierunkiem przepływu a wiązką nadawczą.<br><br>
 +
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
 +
<br><br><math>\Delta f = f_s-f_0 = f_s</math><br><br><br> Ś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.<br>
 
Średnia częstość kołowa <math>\bar{\omega}</math> widma dopplerowskiego <math>P(\omega)</math> może być zdefiniowana jako <br><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>
 
<math>\bar{\omega}=\frac{\int^{\infty}_{-\infty}\omega P(\omega)\mbox{d}\omega}{P(\omega)\mbox{d}\omega}</math>
Linia 52: Linia 57:
 
<br><br><math>R(0)=A(0) </math><br><br>
 
<br><br><math>R(0)=A(0) </math><br><br>
 
Korzystając ze wcześniejszych równań dostajemy
 
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>
+
<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>
 
gdzie <math>T_{PRF}</math> - czas między kolejnymi strzałami. <br>
Praca Miller
+
====Estymator Millera-Rochwargera====
 +
Wartość <math>R(T_{PRF})</math> estymować możemy np. w oparciu o estymator Millera-Rochwargera <ref>Miller, Kenneth, and M. Rochwarger. "A covariance approach to spectral moment estimation." IEEE Transactions on Information Theory 18.5 (1972): 588-596.</ref>. Jest to [[wikipedia:Maximum_likelihood_estimation|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>
 +
<br><br><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><br><br>
 +
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
 +
<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>
  
 
===Zadanie===
 
===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>
+
Proszę przygotować funkcję estymującą częstotliwość średnią w oparciu o powyższy estymator dla jednowymiarowego zespolonego sygnału.
<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>
 
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:<br><br>
 
<math>\Delta f=\frac{2f_{0}v\mbox{cos}\theta}{c}</math><br><br>
 
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==
 
==Prezentacja Kolor==

Wersja z 10:48, 30 cze 2016

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

Rys teoretyczny

Dotychczas traktowaliśmy nasze obrazy jako dwuwymiarowy sygnał. Będziemy jednak dążyć do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie obszaru pomiarowego. Wymagać to będzie od nas uwzględnienia wymiaru czasowego. 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 niepunktowego 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.

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)

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)
  1. Kasai, Chihiro, et al. "Real-time two-dimensional blood flow imaging using an autocorrelation technique." IEEE Trans. Sonics Ultrason 32.3 (1985): 458-464.
  2. Miller, Kenneth, and M. Rochwarger. "A covariance approach to spectral moment estimation." IEEE Transactions on Information Theory 18.5 (1972): 588-596.