Brain-wiki - Wkład użytkownika [pl]

USG/Projekty

2016-12-14T19:10:41Z

Mlew: /* Wstęp */

=Wstęp=
Zaliczenie warsztatów wymaga realizacji jednego z poniżej zaproponowanych projektów zespołowych (2-3 os.). Celem projektów jest implementacja i przetestowanie wybranej metody obrazowania/przetwarzania sygnału USG.
Studenci mają samodzielnie zaimplementować algorytm korzystając z powołanych publikacji lub źródeł internetowych. W zależności od tematu, algorytmy będą testowane bezpośrednio na edukacyjnym systemie ultrasonografu lub na dostarczonych danych. Istnieje także możliwość samodzielnego zebrania danych z fantomu tkankowego.
Na zakończenie każdy zespół przygotuje i wygłosi krótką prezentację (10-15 slajdów), która przedstawi zrealizowany projekt i osiągnięte rezultaty.

==Implementacja OpenCL metod obrazowych==
Zadanie polega na zaimplementowaniu i optymalizacji w OpenCL jednej z metod rekonstrukcji obrazu B-mode: metody klasycznego beamformingu, metody PWI, metody STA. Rezultaty algorytmu równoległego w OpenCL muszą zostać zweryfikowane i porównane z implementacją w Python.
Integracja zaimplementowanego kodu ze skryptem służącym do odbierania danych surowych z edukacyjnego systemu USG, powinna pozwolić na uzyskanie obrazowania w czasie rzeczywistym.

==Short-lag Spatial Coherence==
Short-lag Spatial Coherence (SLSC) jest alternatywną metodą obrazowania tkanek, względem klasycznej prezentacji B-Mode. W odróżnieniu od standardowego podejścia, obraz nie przedstawia energii z rozproszenia, a stopień koherencji (spójności) przestrzennej pomiędzy sygnałami rejestrowanymi przez różne przetworniki odbiorcze. Projekt polegać będzie na zaimplementowaniu prostego toru obliczania obrazu SLSC z danych RF na podstawie pracy.<ref>Lediju, Muyinatu A., et al. "Short-lag spatial coherence of backscattered echoes: Imaging characteristics." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 58.7 (2011): 1377-1388.</ref>

==Obrazowanie tłumienia==
Na zajęciach badane było przesunięcie pasma częstotliwości sygnału odebranego względem sygnału nadanego powstałe w wyniku efektu Dopplera. Pewne przesunięcie widma obserwuje się jednak również dla obiektów statycznych i jest ono zależne od współczynnika tłumienia tkanki obrazowej. Efekt ten pozwala potencjalnie na tworzenie mapy współczynnika tłumienia i w konsekwencji - uzyskanie dodatkowej informacji na temat obrazowanej tkanki. Projekt polegać będzie na zaimplementowaniu jednej z metod tworzenia takiej mapy w oparciu o estymator częstotliwości średniej.
Przedstawienie metod obrazowania tłumienia można znaleźć w pracach dr Z.Klimondy.<ref>http://www.ippt.pan.pl/staff/zklim</ref>

==Obrazowanie metodą Synthetic Transmit Aperture==
Schemat Synthetic Transmit Aperture (STA) polega na schemacie nadawczym-odbiorczym: jeden przetwornik nadaje, wszystkie odbierają. Rekonstrukcja obrazu B-mode odbywa się w sposób analogiczny jak w metodzie PWI omówionej na zajęciach. Dokładny opis metody i wzory potrzebne do rekonstrukcji można znaleźć w pracy.<ref>Ihor Trots, Andrzej Nowicki, Marcin Lewandowski and Yuriy Tasinkevych (2011). Synthetic Aperture Method in Ultrasound Imaging, Ultrasound Imaging, Mr Masayuki Tanabe (Ed.), InTech, DOI: 10.5772/15986. Available from: http://www.intechopen.com/books/ultrasound-imaging/synthetic-aperture-method-in-ultrasound-imaging</ref>

==Coś z Dopplera?==

----

2016-12-14T18:39:52Z

Mlew: /* Wstęp */

=Wstęp=

'''Warsztaty z metod obrazowania ultradźwiękowego''' są autorskim programem<ref>z inicjatywy dr Marcina Lewandowskiego z Instytutu Podstawowych Problemów Techniki PAN; opracowanie materiałów przy współpracy z Tomaszem Steiferem</ref> obejmującym metody i algorytmy obrazowania USG, budowę i funkcje aparatury USG oraz wprowadzenie do programowania równoległego na procesorach GPU w środowisku OpenCL.
Dzięki dofinansowaniu otrzymanemu z ''Funduszu Innowacji Dydaktycznych'' UW opracowano materiały oraz zakupiono edukacyjny system ultrasonografu, który umożliwia zbierania surowych sygnałów i testowanie własnych algorytmów przetwarzania sygnałów.

===Opis===
Ultrasonografia stanowi najpowszechniej stosowaną modalność diagnostyczną we współczesnej medycynie. Kurs ma na celu zapoznanie studentów z szerokim spektrum metod i zastosowań ultradźwięków w diagnostyce medycznej. Metody te obejmują standardowe obrazowanie USG, obrazowanie parametryczne, elastografię, metody oceny przepływu krwi i inne. Zajęcia będą podzielone na część wykładową (ok. 1/3 czasu) i warsztatową (ok. 2/3 czasu).

Zagadnienia części wykładowej:
#Ultradźwięki – zagadnienia generacji, propagacji fal ultradźwiękowych, zjawiska falowe.
#Zagadnienia tworzenia obrazu USG; metody elektronicznego sterowania i ogniskowania wiązki; metoda beamformingu.
#Nowe metody syntetycznej apertury w obrazowaniu; obrazowanie 3D/4D.
#Wybrane inne techniki obrazowe (np. elastografia, obrazowanie tłumienia).
#Podstawy metod oceny i pomiaru przepływu krwi.
#Prezentacja wybranych zastosowań USG i metod dopplerowskich w praktyce medycznej.
#Programowanie równoległe, architektura procesorów CPU/GPU, środowiska CUDA/OpenCL, narzędzia.
#Budowa i funkcje aparatury USG, zagadnienia zapewnienia jakości, certyfikacji wyrobów.

W części warsztatowej studenci będą zespołowo implementować wybrane metody przetwarzania - m.in. klasyczną rekonstrukcję obrazu, rekonstrukcję w obrazowaniu falą płaską czy obrazowanie prędkości przepływu metodą dopplerowską. Studenci zapoznają się z algorytmami obróbki sygnałów ultradźwiękowych oraz zagadnieniami ich implementacji i optymalizacji. Praca w części warsztatowej odbywać będzie się w grupach (2-3 osoby).
Zwieńczeniem kursu będzie realizacja projektów zespołowych. Projekty te będą obejmować implementację własnych algorytmów przetwarzania surowych sygnałów ech ultradźwiękowych w oparciu o samodzielnie zebrane dane z uniwersalnej platformy ultrasonografu. Podstawą uzyskania pozytywnej oceny będzie wykonanie zadań z części warsztatowej oraz zaliczenie projektu (raport lub prezentacja).

===Oprogramowanie i konfiguracja do ćwiczeń===
Ćwiczenia przygotowane zostały pod język python. Korzystać będziemy z bibliotek:
#numpy
#scipy
#PIL http://www.pythonware.com/products/pil/
#PyOpenCL

Ponadto, konieczne jest zainstalowanie SDK OpenCL i sterowników dla posiadanych procesorów lub kart graficznych (np. Intel, AMD, Nvidia) oraz kompilatora (przy pracy na Windowsie prostym rozwiązaniem może być zainstalowanie darmowej wersji Visual Studio <ref>https://www.visualstudio.com/pl-pl/visual-studio-homepage-vs.aspx</ref>). Zaleca się zainstalowanie w pierwszej kolejności kompilatora oraz SDK OpenCL, a dopiero następnie instalowanie biblioteki PyOpenCL do Pythona.

==Wykłady==

#[[USG/Wyklad_Fizyka_USG|Podstawy fizyczne USG]]
#[[USG/Wyklad_Metody_USG|Metody USG]]
#[[USG/Wyklad_Aparatura_i_Aplikacje|Aparatura i aplikacje medyczne]]
#[[USG/Wyklad_OpenCL|Programowanie równoległe procesorów GPU w OpenCL]]

==Ćwiczenia==

#[[USG/Klasyczna_rekonstrukcja|Klasyczna rekonstrukcja obrazu]]
#[[USG/PWI|Obrazowanie falą płaską i rozbieżną]]
#[[USG/Doppler|Obrazowanie prędkości metodą dopplerowską]]
#[[USG/Parametryczne|Obrazowanie prędkości dźwięku]]
#[[USG/GPU|Wstęp do obliczeń równoległych na GPU]]

==Rozwiązania do ćwiczeń==
*[[USG/Rozwiązania|Rozwiązania]]

==Projekty zaliczeniowe==
*[[USG/Projekty|Propozycje tematów zaliczeniowych]]

----

2016-11-16T18:02:02Z

Mlew:

USG/Wyklad OpenCL

2016-11-16T17:49:57Z

Mlew: Utworzenie

USG

2016-11-16T17:36:24Z

Mlew:

=Pracownia USG=

===Wstępna konfiguracja===
Ćwiczenia przygotowane zostały pod język python. Korzystać będziemy z bibliotek:
#numpy
#scipy
#PIL http://www.pythonware.com/products/pil/
#PyOpenCL

Ponadto, konieczne jest zainstalowanie SDK OpenCL i sterowników dla posiadanych procesorów lub kart graficznych (np. Intel, AMD, Nvidia) oraz kompilatora (przy pracy na Windowsie prostym rozwiązaniem może być zainstalowanie darmowej wersji Visual Studio <ref>https://www.visualstudio.com/pl-pl/visual-studio-homepage-vs.aspx</ref>). Zaleca się zainstalowanie w pierwszej kolejności kompilatora oraz SDK OpenCL, a dopiero następnie instalowanie biblioteki PyOpenCL do Pythona.

==Wykłady==

#[[USG/Wyklad_Fizyka_USG|Podstawy fizyczne USG]]
#[[USG/Wyklad_Metody_USG|Metody USG]]
#[[USG/Wyklad_Aparatura_i_Aplikacje|Aparatura i aplikacje medyczne]]
#[[USG/Wyklad_OpenCL|Programowanie równoległe procesorów GPU w OpenCL]]

==Ćwiczenia==

#[[USG/Klasyczna_rekonstrukcja|Klasyczna rekonstrukcja obrazu]]
#[[USG/PWI|Obrazowanie falą płaską i rozbieżną]]
#[[USG/Doppler|Obrazowanie prędkości metodą dopplerowską]]
#[[USG/Parametryczne|Obrazowanie prędkości dźwięku]]
#[[USG/GPU|Wstęp do obliczeń równoległych na GPU]]

==Projekty zaliczeniowe==
#[[USG/Projekty|Propozycje tematów zaliczeniowych]]

==Rozwiązania do ćwiczeń==
#[[USG/Rozwiązania|Rozwiązania]]

USG

2016-11-16T17:28:32Z

Mlew:

USG/PWI

2016-07-27T17:41:27Z

Mlew: /* Obrazowanie falą płaską */

==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).

[[Plik:Rys_falaplaska.png|200px|thumb|right|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:

<source lang = python>
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]

na = 3 # Liczba nadań
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ń

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

===Surowe dane RF===
Do dyspozycji mamy ponownie dwa pliki z pomiarów na fantomie nitkowym ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg2_pwi_nitki.npy usg2_nitki.npy]) i cystowym ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg2_pwi_cysty.npy 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.

<source lang = python>
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,ar)
R0 = np.arange(mindepth,maxdepth,br)
Grid_div_x,Grid_div_z = np.meshgrid(X0,R0,indexing='xy') # tablice ze współrzędnymi
</source>

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
# 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ą.
# 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:
# 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.
# 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.

# Zaimplementować procedurę uśredniającą obrazy otrzymane z kilku kątów nadawczych.
# 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===
#Tanter, Mickael, and Mathias Fink. "Ultrafast imaging in biomedical ultrasound." Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 61.1 (2014): 102-119.
#Cikes, Maja, et al. "Ultrafast cardiac ultrasound imaging: Technical principles, applications, and clinical benefits." JACC: Cardiovascular Imaging 7.8 (2014): 812-823.
#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.
#Smith, Warren J. Modern optical engineering. Tata McGraw-Hill Education, 1966.

USG/PWI

2016-07-27T17:40:19Z

Mlew: /* Obrazowanie falą płaską */

==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). 
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.
[[Plik:Rys_falaplaska.png|200px|thumb|right|Schemat ilustrujący nadanie falą płaską pod różnymi kątami]]
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:

<source lang = python>
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]

na = 3 # Liczba nadań
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ń

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

===Surowe dane RF===
Do dyspozycji mamy ponownie dwa pliki z pomiarów na fantomie nitkowym ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg2_pwi_nitki.npy usg2_nitki.npy]) i cystowym ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg2_pwi_cysty.npy 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.

<source lang = python>
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,ar)
R0 = np.arange(mindepth,maxdepth,br)
Grid_div_x,Grid_div_z = np.meshgrid(X0,R0,indexing='xy') # tablice ze współrzędnymi
</source>

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
# 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ą.
# 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:
# 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.
# 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.

# Zaimplementować procedurę uśredniającą obrazy otrzymane z kilku kątów nadawczych.
# 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===
#Tanter, Mickael, and Mathias Fink. "Ultrafast imaging in biomedical ultrasound." Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 61.1 (2014): 102-119.
#Cikes, Maja, et al. "Ultrafast cardiac ultrasound imaging: Technical principles, applications, and clinical benefits." JACC: Cardiovascular Imaging 7.8 (2014): 812-823.
#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.
#Smith, Warren J. Modern optical engineering. Tata McGraw-Hill Education, 1966.

Plik:Rys falaplaska.png

2016-07-27T17:39:25Z

Mlew:

USG/Doppler

2016-07-05T07:41:26Z

Mlew:

=Metoda dopplerowska impulsowa (Pulsed Wave)=
Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:

<source lang = python>
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]
</source>

Obrazowany jest przekrój fantomu przepływowego złożonego z rurek umieszczonych w materiale tkankopodobnym.
[[Plik:Rys_fantom_flow.png|200px|thumb|right|Zdjęcie fantomu przepływowego.]]
[[Plik:Rys_fantom_flow_przekroj.png|200px|thumb|right|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 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.
[[Plik:Rys_demodulation.png|200px|thumb|right|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żdej 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). Dla naszych danych wystarczający powinien być filtr o częstotliwości odcięcia równej 5.55MHz.

===Zadanie===
# 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ę z pewnej objętości pomiarowej, co wynika z kształtu wiązki ultradźwiękowej).
Ś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 [[Fizyka_III/Fale_dźwiękowe#Efekt_Dopplera|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<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

<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:

<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 <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>

<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 sygnału zespolonego.

==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:

<source lang=python>
# 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)
</source>

===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:

<source lang=python>
scipy.signal.medfilt2d(array, K)
</source>

====Zadanie====
Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji).

----

USG/GPU

2016-07-03T12:31:44Z

Mlew:

=Wstęp do obliczeń równoległych na GPU=
W tej części ćwiczeń stworzymy pierwszy program wykorzystujący bibliotekę OpenCL do obliczeń na kartach graficznych.
Posłużymy się pythonowym wrapperem do OpenCL - biblioteką PyOpencl. Same kernele będziemy pisać w języku C, ale wszystkie dodatkowe procedury będziemy mogli pisać przy użyciu standardowych poleceń Pythona.
Zaczniemy od stworzenia prostego programu sumującego dwie tablice liczb zmiennoprzecinkowych i następnie mnożącego tę sumę przez pewien ustalony mnożnik.
Na sam początek importujemy bibliotekę:

<source lang=python>
import pyopencl as cl
</source>

Kod OpenCL przechowywać będziemy w postaci zmiennej tekstowej. Stworzymy prosty kernel, który dostaje wejściu dwie tablice o wymiarach odpowiednio <math>N\times M\times 2</math> i <math>N\times M</math> oraz zmienną z mnożnikiem. Kernel wykorzystywać będzie <math>N\cdot M</math> jednostek obliczeniowych. Pierwsza tablica zawiera nasze dane wejściowe (ostatni wymiar koduje indeks tablicy wejściowej); druga tablica będzie tablicą w której przechowywać będziemy wynik obliczeń. Przykładowy kernel zdefiniowany jest poniżej:

<source lang=python>
kernels="
__kernel void Sum(__global float *in, __global float *out, const float add)
{
//zaczynamy od zadeklarowania zmiennych
const int n = get_global_id(0); //indeks w pierwszym wymiarze danej jednostki roboczej
const int m = get_global_id(1); //indeks w drugim wymiarze danej jednostki roboczej
const int M = get_global_size(1);
const int nm = n*M + m; // indeks w tablicy wyjsciowej odpowiadajacy wartości o współrzędnych (n,m)
__private int index1; // indeks w zmiennej wejsciowej odpowiadajacy wartości o współrzędnych (n,m,0)
__private int index2; // indeks w zmiennej wejsciowej odpowiadajacy wartości o współrzędnych (n,m,1)
index1 = 0 + m * 2 + n * M * 2;
index2 = 1 + m * 2 + n * M * 2;
out[nm] = (in[index1]+in[index2]) * add;
}
"
</source>

W celu wywołania kernela będziemy musieli przygotować środowisko OpenCL, skompilować kod OpenCL (za pomocą odpowiednich funkcji w Pythonie) oraz zarezerwować miejsce w pamięci i przesłać dane do bufora pamięci na GPU.
Na początek sprawdzimy jakie procesory mamy do dyspozycji

<source lang=python>
platform = cl.get_platforms()
</source>

Proszę wypisać zmienną platform. Zależnie od sprzętu i oprogramowania do dyspozycji będziemy mieli jedną lub więcej platform (gdy np. mamy proces Intela i kartę graficzną AMD) wraz z określoną liczbą urządzeń na każdej platformie. Dostępne urządzenia na pierwszej platformie możemy podejrzeć wywołując:

<source lang=python>
print(platform[0].get_devices())
</source>

Wybierzemy jedno z urządzeń i skompilujemy nasz kernel do obliczeń na tym urządzeniu

<source lang=python>
my_device = [platform[0].get_devices()[0]] # ta czesc kodu musi byc dostosowana do sprzetu na ktorym sa prowadzone cwiczenia
ctx = cl.Context(my_gpu_device)
queue = cl.CommandQueue(ctx, device=my_gpu_device[0])
mod = cl.Program(ctx,kernels).build()
</source>

Musimy jeszcze przygotować dane. Na początek zdefiniujemy przykładowe dane - pamiętaj przy tym, że obliczenia wykonywać możemy na liczbach pojedynczej precyzji.

<source lang=python>
in = np.zeros((300,200,2))
in[:,:,0] = np.arange((300*200)).reshape((300,200))
in[:,:,1] = np.arange((300*200)).reshape((300,200)) - 100
in = in.astype(np.float32) # zmieniamy dane na pojedynczą precyzję
out = np.zeros((300,200)).astype(np.float32)
add = np.float32(add)
</source>

Następnie zarezerwujemy pamięć, przenosząc od razu dane wejściowe. Poniższe zmienne stanowią rodzaj wskaźnika do interesujących nas miejsc w pamięci urządzenia.

<source lang=python>
input = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=in)
output = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=out)
</source>

Jesteśmy już gotowi do wywołania naszego kernela:

<source lang=python>
event = mod.Sum(queue, (np.int32(300), np.int32(200), 1), (1,1,1), input, output, add)
# jak widać nasz kernel możemy wywołać jako metodę zmiennej mod. przyjmuje ona na wejściu potok (?) obliczeniowy, całkowite wymiary grup roboczych, wymiary ???? oraz wskaźniki do buforów)
event.wait()
</source>

Po wywołaniu kernela pozostaje nam tylko przenieść z powrotem dane wynikowe do pamięci hosta (tak aby były one dostępne do dalszych, standardowych obliczeń z poziomu Pythona).

<source lang=python>
event = cl.enqueue_copy(queue, out, output)
event.wait()
</source>

Proszę sprawdzić dane wynikowe i porównać je z analogicznymi danymi otrzymanymi z obliczeń liniowych:

<source lang=python>
out_liniowy = in[:,:,0] + in[:,:,1]
</source>

===Zadanie===
Proszę stworzyć kernel, który będzie dostawać dwie tablice liczb zmiennoprzecinkowych, a następnie liczyć sufit (ang. ceil) z wartości bezwzględnych elementów obu tablic i zwracać największy wspólny dzielnik elementów tych tablic. Porównać działanie (wyniki) programu z poniższą liniową implementacją:

<source lang=python>
...
</source>

==Profilowanie==
Główną przyczyną naszego zainteresowania obliczeniami na GPU jest potencjalny spadek czasu wykonywania wymaganych obliczeń. Z tego punktu widzenia istotna jest możliwość dokładnego określenia czasu obliczeń i porównania algorytmów liniowych z równoległymi. Praktyczny czas trwania obliczeń jest zależny od implementacji sprzętowej, parametrów komputera itp. Nie jest to pojęcie tożsame ze [[wikipedia:Złożoność obliczeniowa | złożonością obliczeniową]]. 
Do mierzenia czasu obliczeń liniowych wykorzystamy standardową bibliotekę time. Aby wykonać pojedynczy pomiar czasu działania funkcji y(x) należy wywołać kod poniższej postaci

<source lang=python>
from time import clock
start = clock() # wskazanie zegara w danej chwili
y(x)
end = clock()
t=end-start # zmierzony czas
</source>

Oczywiście, pomiar taki jest obarczony błędem pomiarowym. Z tego względu rozsądne szacunki należy wykonywać w oparciu o dużą liczbę pomiarów. 
Biblioteka time niestety nie współpracuje dobrze z obliczeniami równoległymi - zawyżając wartość efektywnego czasu obliczeń. Z tego względu, OpenCL dostarcza własnych narzędzi do profilowania. Aby możliwe było profilowanie danej kolejki należy deklarując użyć odpowiedniego parametru:

<source lang=python>
queue=cl.CommandQueue(ctx, device=my_gpu_device[0], properties=cl.command_queue_properties.PROFILING_ENABLE)
</source>

Aby zmierzyć czas wykonywania naszego kernela Sum z wcześniejszej części tekstu, wykonujemy

<source lang=python>
event = mod.Sum(queue, (np.int32(300), np.int32(200), 1), (1,1,1), input, output, add)
event.wait()
t = (event.profile.end-event.profile.start)*1e-9 # zmierzony czas
</source>

Zwróćmy uwagę na to, że wynik dotyczy wyłącznie samych obliczeń - nie uwzględnia on czasu potrzebnego na przeniesienie danych z pamięci hosta do pamięci urządzenia i z pamięci urządzenia do hosta. Czas komunikacji między urządzeniami jest względnie długi - sprawia to, że obliczenia równoległe są opłacalne dopiero przy odpowiednio dużych danych wejściowych. W oparciu o poniższy kod możemy zmierzyć, ile czasu trwa skopiowanie zawartości zmiennej in do bufora input w pamięci urządzenia:

<source lang=python>
start_event=cl.enqueue_marker(queue)
input = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=in)
t = (event.profile.end-start_event.profile.start)*1e-9 # zmierzony czas
</source>

Analogicznie zmierzyć możemy czas potrzebny na przeniesienie wyników z pamięci urządzenia do pamięci hosta.

===Zadanie===
#Dla kilku wybranych rozmiarów tablicy danych wejściowych zmierzyć (liczba powtórzeń >10000) czas potrzebny na:
#* przeniesienie danych z hosta do urządzenia;
#* wykonanie kernela Sum;
#* przeniesienie danych z urządzenia do hosta.
#:Czy czasy te rosną liniowo z rozmiarem tablicy?
# Zmierzyć (liczba powtórzeń >10000) czas trwania obliczeń funkcji np.add(x,y) dla kilku wybranych rozmiarów tablic danych wejściowych (rozmiary tożsame z rozmiarami z podpunktu pierwszego).
#:Czy wielkość ta zmienia się liniowo z rozmiarem tablicy?
# Na podstawie powyższego ocenić: (dla konkretnego komputera) od jakiego rozmiaru tablicy zrównoleglenie obliczeń zaczyna być opłacalne?

USG/GPU

2016-07-03T12:27:01Z

Mlew:

=Wstęp do obliczeń równoległych na GPU=
W tej części ćwiczeń stworzymy pierwszy program wykorzystujący bibliotekę OpenCL do obliczeń na kartach graficznych.
Posłużymy się pythonowym wrapperem do OpenCL - biblioteką PyOpencl. Same kernele będziemy pisać w języku C, ale wszystkie dodatkowe procedury będziemy mogli pisać przy użyciu standardowych poleceń Pythona.
Zaczniemy od stworzenia prostego programu sumującego dwie tablice liczb zmiennoprzecinkowych i następnie mnożącego tę sumę przez pewien ustalony mnożnik.
Na sam początek importujemy bibliotekę:

<source lang=python>
import pyopencl as cl
</source>

Kod OpenCL przechowywać będziemy w postaci zmiennej tekstowej. Stworzymy prosty kernel, który dostaje wejściu dwie tablice o wymiarach odpowiednio <math>N\times M\times 2</math> i <math>N\times M</math> oraz zmienną z mnożnikiem. Kernel wykorzystywać będzie <math>N\cdot M</math> jednostek obliczeniowych. Pierwsza tablica zawiera nasze dane wejściowe (ostatni wymiar koduje indeks tablicy wejściowej); druga tablica będzie tablicą w której przechowywać będziemy wynik obliczeń. Przykładowy kernel zdefiniowany jest poniżej:

<source lang=python>
kernels="
__kernel void Sum(__global float *in, __global float *out, const float add)
{
//zaczynamy od zadeklarowania zmiennych
const int n = get_global_id(0); //indeks w pierwszym wymiarze danej jednostki roboczej
const int m = get_global_id(1); //indeks w drugim wymiarze danej jednostki roboczej
const int M = get_global_size(1);
const int nm = n*M + m; // indeks w tablicy wyjsciowej odpowiadajacy wartości o współrzędnych (n,m)
__private int index1; // indeks w zmiennej wejsciowej odpowiadajacy wartości o współrzędnych (n,m,0)
__private int index2; // indeks w zmiennej wejsciowej odpowiadajacy wartości o współrzędnych (n,m,1)
index1 = 0 + m * 2 + n * M * 2;
index2 = 1 + m * 2 + n * M * 2;
out[nm] = (in[index1]+in[index2]) * add;
}
"
</source>

W celu wywołania kernela będziemy musieli przygotować środowisko OpenCL, skompilować kod OpenCL (za pomocą odpowiednich funkcji w Pythonie) oraz zarezerwować miejsce w pamięci i przesłać dane do bufora pamięci na GPU.
Na początek sprawdzimy jakie procesory mamy do dyspozycji

<source lang=python>
platform = cl.get_platforms()
</source>

Proszę wypisać zmienną platform. Zależnie od sprzętu i oprogramowania do dyspozycji będziemy mieli jedną lub więcej platform (gdy np. mamy proces Intela i kartę graficzną AMD) wraz z określoną liczbą urządzeń na każdej platformie. Dostępne urządzenia na pierwszej platformie możemy podejrzeć wywołując:

<source lang=python>
print platform[0].get_devices()
</source>

Wybierzemy jedno z urządzeń i skompilujemy nasz kernel do obliczeń na tym urządzeniu

<source lang=python>
my_device = [platform[0].get_devices()[0]] # ta czesc kodu musi byc dostosowana do sprzetu na ktorym sa prowadzone cwiczenia
ctx = cl.Context(my_gpu_device)
queue=cl.CommandQueue(ctx, device=my_gpu_device[0])
mod = cl.Program(ctx,kernels).build()
</source>

Musimy jeszcze przygotować dane. Na początek zdefiniujemy przykładowe dane - pamiętaj przy tym, że obliczenia wykonywać możemy na liczbach pojedynczej precyzji.

<source lang=python>
in=np.zeros((300,200,2))
in[:,:,0] = np.arange((300*200)).reshape((300,200))
in[:,:,1] = np.arange((300*200)).reshape((300,200)) - 100
in = in.astype(np.float32) # zmieniamy dane na pojedynczą precyzję
out = np.zeros((300,200)).astype(np.float32)
add = np.float32(add)
</source>

Następnie zarezerwujemy pamięć, przenosząc od razu dane wejściowe. Poniższe zmienne stanowią rodzaj wskaźnika do interesujących nas miejsc w pamięci urządzenia.

<source lang=python>
input = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=in)
output = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=out)
</source>

Jesteśmy już gotowi do wywołania naszego kernela:

<source lang=python>
event = mod.Sum(queue, (np.int32(300), np.int32(200), 1), (1,1,1), input, output, add)
# jak widać nasz kernel możemy wywołać jako metodę zmiennej mod. przyjmuje ona na wejściu potok (?) obliczeniowy, całkowite wymiary grup roboczych, wymiary ???? oraz wskaźniki do buforów)
event.wait()
</source>

Po wywołaniu kernela pozostaje nam tylko przenieść z powrotem dane wynikowe do pamięci hosta (tak aby były one dostępne do dalszych, standardowych obliczeń z poziomu Pythona).

<source lang=python>
event = cl.enqueue_copy(queue, out, output)
event.wait()
</source>

Proszę sprawdzić dane wynikowe i porównać je z analogicznymi danymi otrzymanymi z obliczeń liniowych:

<source lang=python>
out_liniowy = in[:,:,0] + in[:,:,1]
</source>

===Zadanie===
Proszę stworzyć kernel, który będzie dostawać dwie tablice liczb zmiennoprzecinkowych, a następnie liczyć sufit (ang. ceil) z wartości bezwzględnych elementów obu tablic i zwracać największy wspólny dzielnik elementów tych tablic. Porównać działanie (wyniki) programu z poniższą liniową implementacją:

<source lang=python>
...
</source>

==Profilowanie==
Główną przyczyną naszego zainteresowania obliczeniami na GPU jest potencjalny spadek czasu wykonywania wymaganych obliczeń. Z tego punktu widzenia istotna jest możliwość dokładnego określenia czasu obliczeń i porównania algorytmów liniowych z równoległymi. Praktyczny czas trwania obliczeń jest zależny od implementacji sprzętowej, parametrów komputera itp. Nie jest to pojęcie tożsame ze [[wikipedia:Złożoność obliczeniowa | złożonością obliczeniową]]. 
Do mierzenia czasu obliczeń liniowych wykorzystamy standardową bibliotekę time. Aby wykonać pojedynczy pomiar czasu działania funkcji y(x) należy wywołać kod poniższej postaci

<source lang=python>
from time import clock
start = clock() #wskazanie zegara w danej chwili
y(x)
end = clock()
t=end-start #zmierzony czas
</source>

Oczywiście, pomiar taki jest obarczony błędem pomiarowym. Z tego względu rozsądne szacunki należy wykonywać w oparciu o dużą liczbę pomiarów. 
Biblioteka time niestety nie współpracuje dobrze z obliczeniami równoległymi - zawyżając wartość efektywnego czasu obliczeń. Z tego względu, OpenCL dostarcza własnych narzędzi do profilowania. Aby możliwe było profilowanie danej kolejki należy deklarując użyć odpowiedniego parametru:

<source lang=python>
queue=cl.CommandQueue(ctx, device=my_gpu_device[0],properties=cl.command_queue_properties.PROFILING_ENABLE)
</source>

Aby zmierzyć czas wykonywania naszego kernela Sum z wcześniejszej części tekstu, wykonujemy

<source lang=python>
event = mod.Sum(queue, (np.int32(300), np.int32(200), 1), (1,1,1), input, output, add)
event.wait()
t=(event.profile.end-event.profile.start)*1e-9 #zmierzony czas
</source>

Zwróćmy uwagę na to, że wynik dotyczy wyłącznie samych obliczeń - nie uwzględnia on czasu potrzebnego na przeniesienie danych z pamięci hosta do pamięci urządzenia i z pamięci urządzenia do hosta. Czas komunikacji między urządzeniami jest względnie długi - sprawia to, że obliczenia równoległe są opłacalne dopiero przy odpowiednio dużych danych wejściowych. W oparciu o poniższy kod możemy zmierzyć, ile czasu trwa skopiowanie zawartości zmiennej in do bufora input w pamięci urządzenia:

<source lang=python>
start_event=cl.enqueue_marker(queue)
input = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=in)
t=(event.profile.end-start_event.profile.start)*1e-9 #zmierzony czas
</source>

Analogicznie zmierzyć możemy czas potrzebny na przeniesienie wyników z pamięci urządzenia do pamięci hosta.

===Zadanie===
#Dla kilku wybranych rozmiarów tablicy danych wejściowych zmierzyć (liczba powtórzeń >10000) czas potrzebny na:
#* przeniesienie danych z hosta do urządzenia;
#* wykonanie kernela Sum;
#* przeniesienie danych z urządzenia do hosta.
#:Czy czasy te rosną liniowo z rozmiarem tablicy?
# Zmierzyć (liczba powtórzeń >10000) czas trwania obliczeń funkcji np.add(x,y) dla kilku wybranych rozmiarów tablic danych wejściowych (rozmiary tożsame z rozmiarami z podpunktu pierwszego).
#:Czy wielkość ta zmienia się liniowo z rozmiarem tablicy?
# Na podstawie powyższego ocenić: (dla konkretnego komputera) od jakiego rozmiaru tablicy zrównoleglenie obliczeń zaczyna być opłacalne?

Mlew:

=Metoda dopplerowska=
Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:
<source lang = python>
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]
</source>

Obrazowany jest przekrój fantomu przepływowego złożonego z rurek umieszczonych w materiale tkankopodobnym.
[[Plik:Rys_fantom_flow.png|200px|thumb|right|Zdjęcie fantomu przepływowego.]]
[[Plik:Rys_fantom_flow_przekroj.png|200px|thumb|right|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 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.
[[Plik:Rys_demodulation.png|200px|thumb|right|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 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 [[Fizyka_III/Fale_dźwiękowe#Efekt_Dopplera|częstotliwość Dopplera]]:
<math>\Delta f=\frac{2 f_{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<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 
<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 
<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. Oznaczenie "PRF" pochodzi od skrótu [[wikipedia:Pulse_repetition_frequency|Pulse Repetition Frequency]] 
====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>
 <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:
<source lang=python>
#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)
</source>
===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:
<source lang=python>
scipy.signal.medfilt2d(array, K)
</source>
====Zadanie====
Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji)

USG/Doppler

2016-07-03T11:12:39Z

Mlew:

=Metoda dopplerowska=
Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:
<source lang = python>
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]
</source>

Obrazowany jest przekrój fantomu przepływowego złożonego z rurek umieszczonych w materiale tkankopodobnym.
[[Plik:Rys_fantom_flow.png|200px|thumb|right|Zdjęcie fantomu przepływowego.]]
[[Plik:Rys_fantom_flow_przekroj.png|200px|thumb|right|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 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.
[[Plik:Rys_demodulation.png|200px|thumb|right|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 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 [[Fizyka_III/Fale_dźwiękowe#Efekt_Dopplera|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<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 
<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 
<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. Oznaczenie "PRF" pochodzi od skrótu [[wikipedia:Pulse_repetition_frequency|Pulse Repetition Frequency]] 
====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>
 <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:
<source lang=python>
#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)
</source>
===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:
<source lang=python>
scipy.signal.medfilt2d(array, K)
</source>
====Zadanie====
Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji)