USG/Klasyczna rekonstrukcja

Klasyczna rekonstrukcja obrazu

Schemat nadawczy jest następujący. W każdym nadaniu przetworniki subapertury nadawczej generują wiązkę o stałym ognisku. Subaperturę nadawczą w [math]n[/math]-tym nadaniu tworzą przetworniki od [math]n[/math]-tego do [math](n+N_{tr})[/math]-tego (gdzie [math]N_{tr}[/math] - to liczba przetworników subapertury nadawczej). Zakładamy, że ognisko położone jest w linii prostej od środka subapertury nadawczej. Kształt wiązki otrzymujemy przez zastosowanie odpowiednich opóźnień nadawczych pomiędzy przetwornikami. Od momentu nadania wszystkie przetworniki (pełna apertura) rejestrują sygnał. Rejestracja sygnału trwa do zadanego momentu (określonego przez głębokość pomiarową jaką chcemy zbadać). W naszym przypadku subapertura nadawcza przesuwana jest w kolejnych chwilach o 1 przetwornik.

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

NT=192 # Liczba przetworników w pełnej aperturze
Ntr=64 # Subapertura nadawcza

Rf1= 40*(1e-3) # Położenie ogniska wiązki nadawczej od środka subapertury nadawczej

Ponadto, do rekonstrukcji wykorzystywać będziemy prędkość dźwięku w ośrodku. Zakładamy, że jest ona stała. Dla pomiarów w fantomie cystowym proszę założyć

c=1540

zaś dla pomiarów w fantomie nitkowym

c=1490 # prędkość dźwięku w wodzie destylowanej

Będziemy korzystać z bibliotek:

import numpy as np
import pylab as py
import scipy.signal as sig
from PIL import ImageTk, Image

Dane RF

Do dyspozycji mamy dwa zestawy danych RF z pomiarów na fantomie nitkowym (plik usg1_nitki.npy) i na fantomie cystowym (plik usg2_cysty.npy). Zaczniemy od wczytania oraz podejrzenia surowych danych RF (przed rekonstrukcją obrazu). Dane możemy wczytać za pomocą poleceń

RF=np.load('usg1_nitki.npy')

Dostaliśmy tablicę o wymiarach [math]NT\times N \times NT-N_{tr} [/math] gdzie [math]N[/math] odpowiada czasowi rejestracji danych z pojedynczego nadania. Proszę podejrzeć dane zebrane dla kilku strzałów (np. pierwszego i ostatniego). Dla przejrzystości proszę podejrzeć jedynie pierwsze 300 próbek:

py.subplot(2,1,1)
py.imshow(RF[:,:300,0]
py.subplot(2,1,2)
py.imshow(RF[:,:300,-1]
py.show()

Jak zmieniają się sygnały między kolejnymi strzałami? Czy na podstawie samych surowych danych widoczna jest oczekiwana struktura (obrazowane obiekty)? Jak zmienia się energia sygnałów dla różnych przetworników odbiorczych zależnie od odległości od środka subapertury?

Opóźnienia nadawczo-odbiorcze

Z danych zebranych dla pojedynczego nadania rekonstruować będziemy pojedynczą linię obrazu. Wykorzystamy do tego sygnału z tych samych przetworników, które nadawały (niepelna apertura odbiorcza). Wartość natężenia pojedynczego piksela otrzymamy sumując sygnały z tych przetworników przesunięte zgodnie z określonymi opóźnieniami odbiorczymi. Sygnał [math]S_{(i,j)}[/math] dla piksela w [math]i[/math]-tej linii i głębokości [math]j[/math] otrzymujemy jako

[math]S_{(i,j)}=\sum^{N_{tr}-1}_{k=0} s_{i,k}(j+t_{i}(k)))[/math]

gdzie [math]s_{i,k}(t)[/math] - sygnał z [math]i[/math]-tego nadania i [math]k[/math]-tego przetwornika w chwili [math]t[/math];
[math]t_{i}(k)[/math] - opóźnienie nadawczo-odbiorcze [math]k[/math]-tego przetwornika w [math]i[/math]-tym nadaniu.

W klasycznej metodzie do rekonstrukcji wykorzystujemy takie same opóźnienia jak te użyte przy nadawaniu. Funkcja powinna dla zadanej liczby przetworników i odległości ogniska położonego w punkcie [math](x,y)[/math] generować tablicę opóźnień, gdzie opóźnienie dla przetwornika położonego w punkcie [math](a,b)[/math] określone jest jako [math]\sqrt{(x-a)^2+(y-b)^2}/c [/math] Proszę pamiętać, że opóźnienia w powyższym wzorze są podane w sekundach - podczas gdy do dalszych obliczeń wygodniejsze może okazać się wykorzystanie wartości liczonych w okresach próbkowania.

Rekonstrukcja obrazu

Mając już do dyspozycji opóźnienia nadawcze dla interesującego nas punktu ogniskowania możemy przystąpić do rekonstrukcji obrazu. Proszę napisać skrypt dokonujący w pętli rekonstrukcji kolejnych linii obrazu.

Zrekonstruowany obraz przedstawiamy w skali decybelowej. Można wykonać to np. przy użyciu przykładowego kodu:

obraz
indices = obraz>0
indices2= obraz<0
indices= indices+indices2 #wybieramy indeksy którym odpowiadają niezerowe wartości - w przeciwnym wypadku operacja logarytmowania byłaby niejednoznaczna
obraz=np.abs(obraz)/np.max(np.abs(obraz[indices]))
obraz[indices] = 10*np.log(obraz[indices])

Zakres dynamiki

Istotny wpływ na wygląd ostatecznego obrazu ma dobór zakresu dynamiki wyświetlanych wartości. Na początek proszę wyświetlić obraz w pełnej dynamice.
Do wyświetlenia obrazu wykorzystamy bibliotekę PIL, która pozwala nam na łatwe interpolowanie obrazu (przeskalowanie osi) - bezpośrednio po rekonstrukcji mamy obraz z zaburzonymi proporcjami.

Frame = Image.fromarray(np.uint8(cm.Greys(obraz.transpose())*255))
width=NT*pitch
depth=N*c/fs
Frame = Frame.resize((int(N*width/depth)*2, N), Image.BICUBIC)
Frame.show()

Następnie proszę rozważyć ograniczenie dynamiki obrazu np. od -60 dB do 0. Ograniczenie dynamiki jest równoważne z zastosowaniem funkcji progowej

indices = obraz <-low
obraz[indices]=-low
indices = obraz >= 0
obraz[indices]=0

Czy na obrazach widoczna jest oczekiwana struktura? Jak zakres dynamiki wpływa na obraz?

Filtrowanie obrazu

Aby otrzymać czytelny obraz często niezbędne jest przefiltrowanie obrazu. Proszę wykonać filtrowanie dolno- i górnoprzepustowe obrazu wzdłuż głębokości, np. według poleceń:

b, a = sig.butter(10, 0.2, 'highpass')    
obrazek = sig.filtfilt(b, a, obrazek,axis=1)
b, a = sig.butter(10, 0.9, 'lowpass')    
obrazek = sig.filtfilt(b, a, obrazek,axis=1)

Jak zmienił się obraz po zastosowaniu filtrów? Proszę zbadać jak parametry filtra wpływają na końcowy obraz.

Obwiednia obrazu

W większości praktycznie stosowanych urządzeń przez obraz BMode rozumie się obwiednię zrekonstruowanego sygnału. Obwiednie można wykonać na różne sposoby; jedna z popularnych metod oparta jest o transformację Hilberta. Przykładowy skrypt zamieniający każdą linię obrazu na jej obwiednię

for p in range(np.shape(RF)[2]):
    obraz[p,:]=np.abs(sig.hilbert(obraz[p,:]))

Proszę porównać obraz zaraz po rekonstrukcji oraz obwiednię obrazu. Czym się one różnią? Czy informacja utracona w czasie przejścia do obwiedni mogła być istotna z punktu widzenia obrazowania medycznego?

Położenie ogniska

W czasie rekonstrukcji założyliśmy, że ognisko wiązki nadawczej położone było zawsze w określonej odległości od środka subapertury nadawczej. Proszę zbadać:

1. Jak zmienia się ostrość obrazu w zależności od głębokości? Na jakiej głębokości rozdzielczość obrazu wydaje się być najlepsza?
2. Proszę zrekonstruować obraz ponownie zakładając inną odległość ogniska nadawczego (np. dwa razy mniejszą). Jak zmienił się obraz?