USG/PWI

Obrazowanie falą płaską

W obrazowaniu fala płaską (ang. PWI - Plane Wave Imaging) wektor opóźnień nadawczych ma postać prostej, co pozwala uformować wiązkę o płaskim froncie falowym. W odróżnieniu do rozważanej wcześniej metody klasycznego beamformingu, wiązka taka jest nieogniskowana. O froncie falowym w PWI możemy myśleć jak o froncie fali pochodzącej ze źródła punktowego (ognisko wirtualne) położonego względnie daleko za aperturą (w tym sensie, że różnice w czasach przejścia od ogniska wirtualnego do różnych przetworników nadawczych są pomijalne).

Schemat ilustrujący nadanie falą płaską pod różnymi kątami

O ile nie ograniczają nas możliwości techniczne sprzętu (np. mała moc obliczeniowa), w PWI staramy się nadawać i odbierać pełną aperturą. Punktem wyjścia jest fala płaska nadawana pod kątem zerowym względem apertury. Takiemu schematowi nadawczemu odpowiada wektor opóźnień określony funkcją stałą (wszystkie przetworniki nadają w tym samym momencie). Obrót wykresu opóźnień nadawczych odpowiada zmianie kąta nadawczego. Do stworzenia obrazu ultradźwiękowego wystarczą nam dane zebrane dla pojedynczego nadania (z jednego kąta). W praktyce jednak wykonuje się kilka pomiarów (nadań, strzałów itp.) dla kilku rożnych kątów - pozwala to uzyskać poprawę zarówno rozdzielczości punktowej obrazu, jak i poprawę kontrastu i zwiększenie stosunku sygnału do szumu (SNR - Signal-to-Noise Ratio).

Parametry potrzebne do dalszej pracy:

f0 = 5.5e6 # Częstotliwość nadawcza [Hz]
fs = 50e6 # Częstotliwość próbkowania [Hz]
pitch = 0.21e-3 # Deklarowana odległość między środkami przetworników nadawczo-odbiorczych (tzw. pitch) [m]

theta = [ -1.74532925e-01 , -1.39626340e-01,  -1.04719755e-01,  -6.98131701e-02,
  -3.49065850e-02,  -2.77555756e-17,   3.49065850e-02 ,  6.98131701e-02,
   1.04719755e-01  , 1.39626340e-01  , 1.74532925e-01] # kąty dla kolejnych nadań
na = len(theta) # Liczba nadań

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

Surowe dane RF

Do dyspozycji mamy ponownie dwa pliki z pomiarów na fantomie nitkowym (usg2_nitki.npy) i cystowym (usg2_cysty.npy). W każdym pliku znajduje się trójwymiarowa tablica [math]NT\times N \times na [/math]
gdzie [math]N[/math] odpowiada czasowi rejestracji (maksymalnej głębokości obrazowanej) danych z pojedynczego nadania.
Kluczowe dla dalszych obliczeń będzie wyznaczenie chwili zerowej, czyli momentu rozpoczęcia akwizycji danych. Na różnych urządzeniach przyjmowane są różne konwencje - czasami akwizycja danych rozpoczyna się w momencie, gdy rozpoczyna się nadawanie; czasami akwizycja rozpoczyna się jeszcze przed nadawaniem.
Zadanie: proszę wyświetlić surowe dane RF (funkcja plot) dla kilku kolejnych nadań i spróbować określić (w przybliżeniu) w jakim momencie rozpoczyna się akwizycja danych: a) przed rozpoczęciem nadawania b) w momencie rozpoczęcia nadawania c) w chwili gdy nadała połowa przetworników

Rekonstrukcja obrazu z pojedynczego strzału

W przypadku klasycznej rekonstrukcji dokonywaliśmy najpierw przesunięcia w czasie próbek RF, a następnie przesunięty (syntezowany) obraz interpolowaliśmy do siatki odpowiadającej rzeczywistym proporcjom obrazowanej struktury. W przypadku PWI wygodniej będzie nam najpierw określić siatkę, a następnie dla kolejnych punktów na tej siatce liczyć wartość pikseli. Poniższy kod wygeneruje dwie tablice, w których przechowywane są współrzędne poziome i pionowe w [m] kolejnych punktów siatki. Siatka jest generowana przy założeniu, że punkt (0,0) znajduje się w środku apertury. Cztery próbki na długość fali powinny wystarczyć do otrzymania wystarczająco dobrej jakości obrazu.

halfWidth = 0.02 # Połowa szerokości obrazowania [m] 
lmbda = c/fs
step=lmbda/4.
maxdepth = N*c/(fs*2) # maksymalna głębokość [m]
mindepth = 0.005 # minimalna głębokość [m]

X0 = np.arange(-halfWidth,halfWidth,step)
R0 = np.arange(mindepth,maxdepth,step)
Grid_div_x,Grid_div_z = np.meshgrid(X0,R0,indexing='xy') # tablice ze współrzędnymi

Podobnie jak wcześniej sygnał w punkcie (x,y) będziemy liczyć jako sumę sygnałów z każdego przetwornika odbiorczego odpowiednio opóźnionych:

[math]S(x,y)=\sum^{NT-1}_{k=0} s_k(t_{tr}(x,y)+t_k(x,y)) [/math],

gdzie [math]s_k(i)[/math] to sygnał z danego przetwornika odbiorczego; [math]t_{tr}(x,y)[/math] czas przejścia fali do punktu obrazowania; [math]t_k(x,y)[/math] czas przejścia fali odbitej z punktu do przetwornika [math]k[/math].
Dla ustalonego przetwornika [math]k[/math] i punktu (x,y) opóźnienie liczone jest jako czas przejścia fali od początku nadania do punktu i z powrotem do danego przetwornika odbiorczego. W przypadku PWI możemy założyć, że dla kąta nadawczego [math]\theta[/math]

[math]t_{tr}(x,y) = (x\cdot sin(\theta)+y\cdot cos(\theta))/c + t_{start} + d [/math],

gdzie [math]d[/math] jest pewnym przesunięciem stałym dla danego nadania - wyznaczonym przez położenie chwili zerowej (moment rozpoczęcia akwizycji).

Do wykonania

1. Stworzyć funkcję, która dla zadanej siatki obliczać będzie wartości energii dla pojedynczego strzału z danego kąta. Podobnie jak wcześniej, interesuje nas obraz w skali decybelowej z ograniczoną dynamiką.
2. Porównać obraz syntezowany dla zadanej wcześniej siatki z obrazem syntezowanym dla siatki o dwukrotnie większych odległościach między punktami

Apodyzacja po stronie odbiorczej

Standardową metodą poprawy jakości obrazu jest zastosowanie apodyzacji po stronie odbiorczej. W czasie rekonstrukcji natężenie w danym punkcie liczyliśmy jako sumę sygnałów z pojedynczych przetworników odbiorczych. Sygnały z pojedynczych przetworników zawierają jednak informację pochodzącą nie tylko z rozproszenia w danym punkcie. Ponadto, dla dwóch różnych przetworników odbiorczych energia sygnału pochodząca z rozproszenia z danego punktu niekoniecznie musi stanowić taką samą cześć całkowitej energii sygnału. W celu zbadania tego zjawiska:

1. Stworzyć funkcję, która dla zadanego punktu liczy sygnał jako sumę sygnałów tylko z wybranej części przetworników odbiorczych (np. wybranych 12 przetworników). Proszę porównać otrzymane obrazy dla kilku wybranych subapertur odbiorczych (np. pierwsze 12 przetworników, środkowe 12 przetworników i ostatnie 12 przetworników - licząc od lewego końca apertury). Jak bardzo różnią się te obrazy między sobą i jak różnią się od obrazu otrzymanego z pełnej apertury odbiorczej?

Możemy zauważyć, że więcej informacji o danym punkcie zawierają przetworniki położone bliżej tego punktu niż przetworniki położone dalej. Apodyzacja po stronie odbiorczej polega na zastosowaniu przy rekonstrukcji wag, które będą odpowiednio proporcjonalne do odległości przetwornika od punktu pomiarowego.

1. Proszę zaimplementować zmodyfikowany algorytm rekonstrukcji wykorzystujący jakąś wybraną formę apodyzacji odbiorczej, gdzie sumowane sygnały z pojedynczych przetworników odbiorczych będą ważone wartościami niemalejącej funkcji od odległości przetwornika od punktu pomiarowego.

Złożenie (compounding) obrazów z nadań pod różnymi kątami

Jednym z prostych sposobów poprawienia obrazu jest złożenie ze sobą (np. przez średnią arytmetyczną) kilku obrazów z różnych kątów nadawczych.

1. Zaimplementować procedurę uśredniającą obrazy otrzymane z kilku kątów nadawczych.
2. Zbadać obraz zrekonstruowany dla wybranego niezerowego kąta nadawczego i spróbować odpowiedzieć na pytanie: czy dla wszystkich punktów na siatce rekonstrukcja jest sensowna? Czy są punkty, które należałoby pominąć w procesie rekonstrukcji/sumowania przy składaniu obrazów z wielu kątów?

Literatura

1. Tanter, Mickael, and Mathias Fink. "Ultrafast imaging in biomedical ultrasound." Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 61.1 (2014): 102-119.
2. Cikes, Maja, et al. "Ultrafast cardiac ultrasound imaging: Technical principles, applications, and clinical benefits." JACC: Cardiovascular Imaging 7.8 (2014): 812-823.
3. Montaldo, Gabriel, et al. "Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography." Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 56.3 (2009): 489-506.
4. Smith, Warren J. Modern optical engineering. Tata McGraw-Hill Education, 1966.