USG/PWI

Obrazowanie falą płaską

W obrazowaniu fala płaską (PWI - Plane Wave Imaging) wykres opóźnień nadawczych względem przetworników ma postać prostej - przetworniki odpalane są w równych odstępach czasu, aby uformować wiązkę o (w przybliżeniu) płaskim froncie falowym. W przeciwieństwie do rozważanej wcześniej metody klasycznej, wiązka taka jest nieogniskowana - stąd oczekujemy podobnej ostrości / rozdzielczości na całej głębokości pomiarowej. O froncie falowym w PWI możemy myśleć jak o froncie fali pochodzącej ze źródła punktowego (ognisko wirtualne) oddalonego istotnie daleko za aperturą.
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 sytuacja w której wszystkie przetworniki nadają jednocześnie. Obrót wykresu opóźnień nadawczych odpowiada przesunięciu ogniska wirtualnego.

Parametry potrzebne do dalszej pracy:

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

na = 3 # Liczba nadań
theta=[-10,0,10] # kąty dla kolejnych 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. Tak jak na poprzednich ćwiczeniach, proszę zacząć od podejrzenia tablicy RF dla kilku kolejnych nadań. Czy z samych surowych danych RF da się określić kształt frontu falowego?

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ść natężenia. 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.

maxWidth=0.02 # Połowa szerokości obrazowania [m] 
l2mbd = c/(fs*2)
ar=int(2*maxWidth/l2mbd)+1 #Określamy gęstość poziomą siatki zgodnie z twierdzeniem ???
maxdepth=0.0610304 #maksymalna głębokość [m]
mindepth=0.00298 # minimalna głębokość [m]
br=(maxdepth-mindepth)/((maxWidth*2)/ar)

X0 = np.linspace(-maxWidth,maxWidth,ar)
R0 = np.linspace(mindepth,maxdepth,br)
Grid_div_x,Grid_div_z = np.meshgrid(X0,R0,indexing='xy') #tablice ze współrzędnymi

Podobnie jak wcześniej wartość natężenia w punkcie (x,y) będziemy liczyć jako sumę natężeń 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 spowrotem do danego przetwornika odbiorczego. W przypadku PWI możemy założyć, że
[math]t_{tr}(x,y) = (x\cdot sin(\theta)+y\cdot cos(\theta))/c [/math]

Proszę stworzyć funkcję, która dla zadanej siatki generować będzie wartości natężeń dla pojedynczego strzału z danego kąta. Podobnie jak wcześniej, interesuje nas obraz w skali decybelowej z ograniczoną dynamiką.

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.