===== Lab 07 : Ultrasound Project =====

  
===== HW07 Assignment =====

  * **[1 pt]** [[#US Příklady ]] Vyřešte označené příklady.
  * **[4 pts]** [[#Zpracování videosekvence US]] Analýza srdečního tepu z US videa.

==== US Příklady ====

Použité rovnice:

  * šíření (akustického) signálu v látce $c = \frac{1}{\sqrt{\rho K}}$
  * akustická impedance $Z_t = \rho c$
  * odrazivost na rozhraní $R = \frac{Z_k - Z_t}{Z_k + Z_t}$

** Šíření signálu** Nechť máme ultrazvukové vlnění s rycholstí šíření $1540 m/s$ v látce s hustotou $1000\:kg \cdot m^{-3}$. Určete (Youngův) modul pružnosti. Jaká je specifická akustícká impedance dané látky?

// Řešení // Youngův modulus (elasticita) $E$ je převrácená hodnota stlačitelnosti $E = \frac{1}{K}$. Ze vztahu pro šíření signálu získáme
$$ c = \frac{1}{\sqrt{\rho K}} \quad \to \quad K = \frac{1}{c^2 \rho} \quad \text{tím pádem je} \quad E = c^2 \rho = 1540^2 * 1200 = 2.85 \cdot 10^9 Pa = 2.85 GPa $$

** Odrazivost ** Nechť máme měkkou tkáň s ak. impedancí $800\;\text{kRayl}$ a dále kost s ak. impedancí $6\;\text{MRayl}$. Jaká je energie (!) uzv. vlny odražené na rozhraní měkká tkáň -- kost?

// Řešení // Pro jednoduchost budeme uvažovat úplný odraz (tj. vlna se šíří kolmo na rozhraní). Impedanci známe: $Z_t = 0.8 \cdot 10^6 Rayl, Z_k = 6 \cdot 10^6 Rayl$. Koeficient odrazivosti (//poměr amplitudy dopadajícího a odraženého vlnění//) je tím pádem $$R = \frac{Z_k - Z_t}{Z_k + Z_t} = \frac{6 - 0.8}{6 + 0.8} = 0.765$$. Energie vlny je proporcionální s kvadrátem amplitudy, tím pádem je odražená energie $R_e = R^2 = 0.585$

** [HW] Odrazivost tkáně II ** Mějme tkáň s hustotou $900\;\text{kg} \cdot m^{-3}$, ultrazvuk s rychlostí šíření $1540 m \cdot s^{-1}$ a dále kost s ak. impedancí $6.75 \cdot 10^6 Rayl$. Jaká je odrazivost na rozhraní tkáň -- kost?

** [HW] Vzorkovací frekvence** Typický (užitný) dosah ultrazvuku v lidském těle je 40 cm, rychost US je $1540\;m\cdot s^{-1}$. S jakou maximální obrazovou frekvencí (//frame rate//) jsme schopni snímat v této hloubce, pokud na každý obrázek potřebujeme 200 paprsků?


==== Zpracování videosekvence US ====
<note>Protože jste nakonec neměli možnost nasnímat si vlastní data, upravil jsem úlohu pouze na zpracování dat.</note>
Jako vstupní data máme videosekvenci (//odkaz je v MS Teams//), která obsahuje záznam snímání krční tepny ultrazvukem. V rámci domácího cvičení chceme za použití dvou technik analýzy časových signálů -- **rychlé Fourierovy transformace** (FFT) a **autokorelace** určit z této sekvence frekvenci srdečního tepu (údery/min). Podrobný popis cvičení nejprve představuje způsoby, jak získát časově-závislý signál a pak popisuje způsob analýzy.

** Načtení dat**:
  - Načtěte video: <code matlab>
%% Load video
v_file='C:\Echo Images\krkavice_ax.avi';

obj = VideoReader(v_file);
frame0 = obj.readFrame();
</code>
  - Nastavte oblast ořezu (interaktivně) tak, aby obsahovala průřez krkavicí<code matlab>
[fr0, c_region] = imcrop(frame);
</code>
  - Použijte pozici a rozměry uložené oblasti ořezu na zbývající framy<code matlab>
%% Apply cropping to other frames
i = 1;
clear X;
while hasFrame(obj)
    c_frame = imcrop( obj.readFrame(), c_region);
    X(:, :, i) = c_frame(:, :, 1);
    i = i + 1;
end
</code>

** Definice signálu [2 pt]**

Pro každý frame spočteme hodnotu signálu $R(t)$, který bude mít podobnou (ideálně stejnou) charakteristiku jako srdeční tep a ten pak v další části vyhodnotíme. Představíme si dva možné způsoby definice $R(t)$
  - Budeme definovat $R(t)$ jako průměrnou intenzitu ve výřezové oblasti -- tu si zvolíme tak, aby se dostatečná část tepenní stěny dostala při systole **mimo** ořezovou oblast, pak bude při systole průměrná intentita nižší než při diastole, kdy bude celá stěna tepny uvnitř ořezové oblasti a tudíž bude průměrná intenzita vyšší.
  - Plocha uvniř tepny je pro náš úkol nejlepší signál, ale v praxi nás bude limitovat charakteristika data. {{ :courses:zsl:coronal_rg.png?400 |}} Zde si zkusíme použít algoritmus ''regiongrowing'' ({{ :courses:zsl:regiongrowing.zip |}} který pro zadaný vstupní obrázek a zadanou prahovou hranici ''thr'' začne od zadané souřadnice počátečního bodu (např. střed ořezové oblasti, pokud je průřez tepnou ve výřezu vycentrované) přidávat do segmentace sousedící pixely. Sousední pixel se do segmentace zařadí, pokud se od centrálního pixelu liší o méně než je zadaný práh. <code matlab>
J = regiongrowing(double(c_frame), int32(size(c_frame, 1)/2), int32(size(c_frame, 2)/2), threshold);
</code>. Pro stabilnější segmentaci je dobré před spuštěním regiongrowing vstupní frame vyhladit Gaussovským filtrem.
  - Vymyslete si vlastní kritérium.
  
Proveďte extrakci a následnou analýzu signálu pro **dva způsoby** spočtení $R(t)$. 

** Analýza signálu **
Jako vstup máme signál $R(t)$, který svým průběhem kopíruje srdeční tep.

  - Odečtením průběžného průměru spočítáme signál $j(t)$ (metodou //running-mean subtraction//) {{ :courses:zsl:regionsize_over_time.png?nolink&400|}}
     * pro každé $t$ spočteme průměr $R$ přes interval $[t-w, t+w]$ a odečteme od $R(t)$ \begin{equation*}
j(t) = R(t) - \frac{1}{2 \cdot w} \sum_{j=t-w}^{t+w} R(j)\end{equation*}
     * vykreslete průběh $R(t)$ a $j(t)$ do zprávy

  - **[1 pt]** FFT analýza 
     * decompose $j(t)$ with Fourier transform, compute its spectrum $J(f) = \mathcal{F}_1[j(t)]$
     * detect the dominant frequency (if needed, skip physiologically nonsensical frequencies)
     * Plot $|J(t)|$ with $x$-axis showing the heartbeat frequency [beats/min]

  - **[1 pt]** Analýza pomocí autokorelace
     * spočítejte auto-korelaci signálu $j(t)$ <code matlab> [acorr, lag] = xcorr(j(t))</code> 
     * najděte její lokální maxima (//peaks//) <code matlab>[peakval, peakpos] = findpeaks(...);</code>
     * vzdálenost lokálních maxim odpovídá srdečnímu tepu
     * vykresletete spočtenou auto-korelaci signálu $j(t)$
 

(//pozor na jednotky -- signál $j(t)$ měříme v $t$=frames, ale potřebujeme převést na srdeční tep (úderů za minutu). K převodu použijeme hodnotu ''FrameRate''//)


** Užitečné funkce **
  - Pro vizualizaci segmentace jako polo-průhlendnou vrstvu použijeme alpha-masking<code matlab>
imshow(c_frame, 'InitialMag', 'fit');

% Make a truecolor all-green image.
green = cat(3, zeros(size(fIm)), ones(size(fIm)), zeros(size(fIm)));
hold on;
h = imshow(green);
hold off;

% Use the logical values of the region to paint over the input image
set(h, 'AlphaData', J)
</code>