• [1 pt] US exercises Solve denoted exercises.
• [4 pts] Processing US video Estimation of heart rate from US video.

Used formulas:

• relationship between the wavelength and frequency in a matter $\lambda = \frac{c}{f}$
• propagation of acoustic signal in a matter $c = \frac{1}{\sqrt{\rho K}}$ ($K$ is compressibility, the inverse of Young's modulus of elasticity $E$)
• acoustic impedance $Z_t = \rho c$
• intensity reflection coefficient on a boundary $R = \frac{(Z_1 - Z_2)^2}{(Z_1 + Z_2)^2}$

1. Signal propagation

1.1 The time delay between transmission and the arrival of the reflected wave of a signal using ultrasound travelling through a piece of fat tissue was $0.13 ms$. At what depth did this reflection occur? The speed of propagation of ultrasound waves in fat is $1450 m/s$

Solution The wave travels $2 \cdot d$ in $0.13 ms$, therefore $$d = \frac {c_{fat} \cdot t_d} {2} = \frac{1450 m/s \cdot 0.00013 s} {2} = 0.09425 m$$

1.2 Calculate the minimum frequency of ultrasound that will allow you to see details as small as $0.250 mm$ in human tissue. What is the effective depth to which this sound is effective as a diagnostic probe?

Solution Whenever a wave is used as a probe, it is very difficult to detect details smaller than its wavelength $\lambda$. The relationship between wavelength and the frequency is following $$\lambda = \frac{c}{f} -> f = \frac{c}{\lambda} = \frac{1540 m/s}{2.5 \cdot 10^{-4}m} = 6.16 MHz$$ There is a rule of thumb for effective ultrasound depth of penetration which is $500\lambda$, therefore the effective depth of penetration for ultrasound waves with wavelength $0.25 mm$ is $12.5 cm$

1.3 The acoustic (ultrasound) waves propagate with a speed $1540 m/s$ through a matter with a specific weight $1000 kg \cdot m^{-3}$. Determine the Young's modulus of elasticity. What is the specific acoustic impedance of the matter?

Solution Young's modulus of elasticity $E$ is the inverse of compressibility $E = \frac{1}{K}$. From the formula for the speed of signal propagation we obtain $$ c = \frac{1}{\sqrt{\rho K}} \quad \to \quad K = \frac{1}{c^2 \rho} \quad \text{therefore} \quad E = c^2 \rho = 1540^2 * 1000 = 2.37 \cdot 10^9 Pa = 2.37 GPa $$ The acoustic impedance is $$ Z = /rho \cdot c = 1000kg/m^3 \cdot 1540 m/s = 1.54 cdot 10^6 Rayl $$

2 Reflectivity

2.1 In the clinical use of ultrasound, transducers are always coupled to the skin by a thin layer of gel or oil, replacing the air that would otherwise exist between the transducer and the skin. Using the values of acoustic impedance for transducer material $30.8 \cdot 10^6 \frac{kg}{m^2 \cdot s}$, air $429 \frac{kg}{m^2 \cdot s}$ and water $1.5 \cdot 10^6 \frac{kg}{m^2 \cdot s}$. Calculate the intensity reflection coefficient between transducer material and air. Calculate the intensity reflection coefficient between transducer material and gel (assuming for this problem that its acoustic impedance is identical to that of water). Based on the results of your calculations, explain why the gel is used.

Solution The intensity reflection coefficient transducer–air $$R_{tr - air} = \frac{(30.8 \cdot 10^6 - 429)^2}{(30.8 \cdot 10^6 + 429)^2}=0.9999$$ The intensity reflection coefficient transducer–water $$R_{tr - gel} = \frac{(30.8 \cdot 10^6 - 1.5 \cdot 10^6)^2}{(30.8 \cdot 10^6 + 1.5 \cdot 10^6)^2}=0.8228$$ In the case of transducer to air boundary, almost all the ultrasound energy is reflected back, in the case of transducer to gel boundary, the intensity reflection coefficient is not great, but not terrible. In actual use, the gel acoustic impedance would be different than the acoustic impedance of water.

2.2 A soft tissue has the acoustic impedance $800\;\text{kRayl}$ and and a bone $6\;\text{MRayl}$. What is the proportion of reflected energy of the ultrasound wave reflected from the boundary of the soft tissue and the bone?

Solution For simplicity we will consider total reflection ( the wave propagates perpendicularly to the interface). We know the impedance: $Z_t = 0.8 \cdot 10^6 Rayl, Z_b = 6 \cdot 10^6 Rayl$. The coefficient of reflectivity (ratio of the energy of incoming and reflected waves) is therefore $$R = \frac{(Z_b - Z_t)^2}{(Z_b + Z_t)^2} = \frac{(6 - 0.8)^2}{(6 + 0.8)^2} = 0.585$$

3 Miscellaneous

3.1 A bat uses ultrasound to find its way among trees. If this bat can detect echoes $1.00 ms$ apart, what minimum distance between objects can it detect? Could this distance explain the difficulty that bats have finding an open door when they accidentally get into a house? Speed of ultrasound in air is $330 m/s$.

Solution The ability detecting echoes $1ms$ apart implies wavelength $\lambda = \frac{ c \cdot t}{2} = \frac{330 m/s \cdot 10^{-3}s}{2} = 0.15 m$. However this is axial resolution (alongside the ultrasound beam). It would apply if the bat for example circled in the room. The lateral resolution which seems to be more reasonable for this exercise is proportional to the ultrasound frequency, focusing and in medical ultrasound devices to the size of the transducer and distance of the object from the transducer. But who knows how bats are doing these things…

3.2 A dolphin is able to tell in the dark that the ultrasound echoes received from two sharks come from two different objects only if the sharks are separated by $3.50 m$, one being that much farther away than the other. If the ultrasound has a frequency of $100 kHz$, show this ability is not limited by its wavelength. If this ability is due to the dolphin’s ability to detect the arrival times of echoes, what is the minimum time difference the dolphin can perceive?

Solution The wavelength of ultrasound with frequency $100 kHz$ in water, where the waves travel with speed $1540 m/s$, is $ \lambda = \frac{v}{f} = \frac{1540 m/s}{10^5 1/s} = 0.0154m < 3.5m$ So the inability distinguish sharks less than $3.5m$ is not due to the parameters of the dolphin's sonar. The $3.5m$ length travels the wave in $t = \frac{2d}{v} = \frac{2\cdot3.5m}{1540 m/s} = 0.0045s$

[HW] Reflectivity II A tissue has specific weight $900\;\text{kg} \cdot m^{-3}$, ultrasound propagates through the tissue with the speed $1540 m \cdot s^{-1}$ and a bone with acoustic impedance $6.75 \cdot 10^6 Rayl$. What is the reflection coefficient of the tissue – bone boundary?

[HW] Sampling frequency Let's assume that the penetration depth of ultrasound in human body is $20 cm$, speed of ultrasound is $1540m\cdot s^{-1}$. What is the smallest possible sampling frequency (frame rate)we can use for ultrasound imaging if we need 200 rays for each image?

The input videos lab07_videos.zip, contain measurements of carotid artery with US in two different sections - alongside the artery and perpendicular to the artery. In this lab we will try different approaches to frequency analysis – fast Fourier transform (FFT), autocorrelation and power spectral density (PSD) – to estimate the heart rate in beats per minute from the carotid US videos. The following guide describes the extraction of the area of carotid artery in each frame of the video and processing of the resulting signal.

Loading data:

1. Load the video:

   %% Load video
   v_file='C:\Echo Images\us_carotid_perp.mp4';
    
   obj = VideoReader(v_file);
   frame0 = obj.readFrame();

1. Crop the frame (interactively) so that it only contains the carotid artery

   [fr0, c_region] = imcrop(frame);

2. Apply the obtained mask for cropping the consecutive frames

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

Extraction of carotid area in each frame [1 pt]

We need to extract a signal $C(t)$, which will have similar (ideally identical) characteristic as the heart pulsations. The signal of heart pulsations will be analysed in the following steps. We will try two approaches to the extraction of $C(t)$

1. In the first case the $C(t)$ will be defined as the average brightness intensity of the cropped frame – we have to crop the image so that an reasonable part of the artery wall moves in and out of the resulting image during the pulsations (systole and diastole). The average brightness intensity will then be lower in systole (the bright artery wall is not in the cropped frame) and higher during diastole (the bright artery wall is in the cropped frame).
2. In the second case we will define the $C(t)$ will be defined as the area of the artery. For this task we will use the region growing segmentation which assigns pixels to a segmentation if they don't differ less than a certain threshold value. The code for region growing from (regiongrowing.zip requires setting the starting pixel for the segmentation (in case of artery area in a cropped frame the centre of the image) nd the threshold value. The artery area should be a better signal for the estimation of the heart rate, however the data quality may be a limiting factor.

   J = regiongrowing(double(c_frame), int32(size(c_frame, 1)/2), int32(size(c_frame, 2)/2), threshold);

   We can obtain more stable segmentation if we filter the cropped image with gaussian filter (blurring)

   gauss_blurr = fspecial('gaussian',HSIZE,SIGMA);
   c_frame_blur = conv2(c_frame,gauss_blur,'valid'); 

   The HSIZE controls the size of the kernel (the bigger the more the resulting image will be blurred) and the SIGMA is standard deviation of the Gaussian distribution (the bigger the more the resulting image will be blurred).

Perform the average brightness extraction for the us_carotid_long.mp4 video and artery area extraction with region growing segmentation for the us_carotid_perp.mp4 video.

Heart pulsation signal processing

The input signal is $C(t)$ which correlates with the heart pulsation.

1. By subtracting the moving average we obtain signal $j(t)$ which will be free of slow fluctuations and an offset.
   • for each $t$ we compute average of $C$ over the interval $[t-w, t+w]$ and subtract it from $C(t)$ \begin{equation*} j(t) = C(t) - \frac{1}{2 \cdot w} \sum_{j=t-w}^{t+w} C(j)\end{equation*}
   • Visualize the signals $C(t)$ and $j(t)$ in the report from the lab.

1. [1 pt] FFT analysis
   • transform $j(t)$ by fast Fourier transform, compute the magnitude spectrum $J(f) = \mathcal{F}_1[j(t)]$
   • detect the dominant frequency (if it's necessary,skip some of the physiologically nonsensical frequencies)
   • visualize $|J(t)|$ with $x$-axis in beats per minute

1. [1 pt] Autocorrelation analysis
   • compute autocorrelationof $j(t)$

      [acorr, lag] = xcorr(j(t))

   • find the local maxima (peaks)

     [peakval, peakpos] = findpeaks(...);

   • the distance between the local maxima is proportional to the heart rate
   • visualize the computed autocorrelation of $j(t)$ (use proper units – signal $j(t)$ is measured in $t$=frames, sampling frequency is the FrameRate property of the loaded video)

1. [1 pt] Analysis with power spectral density PSD
   • both the FFT analysis and the autocorrelation function have certain frawbacks – the FFT spectrum might be noisy and detecting peaks in autocorrelation function (if you would not be able to use findpeaks function) and converting it to frequency is compicated
   • a more robust heart rate estimation can be done using PSD
   • PSD provides power spectrum estimates that eliminates various transitional and other effects in the signal – this is done by averaging out the random component in the power spectrum
   • there are many methods for estimating – PSD can be directly estimated from

     fft(xcorr(j(t)))

   • or by Welch method, where a signal $x(t)$ is split into shorter segments of length $n$ (which may overlap with an overlap $n_{overlap}$ samples), for each shorter segment $x_{segm}(t_{segm}-n/2,t_{segm}+n/2)$ of the original $x(t)$ we estimate the power spectrum by FFT as $X_{segm}(f) = \mathcal{F}_1[x_{segm}(t)]$ and we average all the power spectra $X_{segm}(f)$ over all segments. When splitting the original signal $x(t)$ it is practical to use a window function, which further improves the resulting PSD by limiting higher frequencies. In Matlab the Welch method for estimating PSD is implemented in pwelch function

     [Pxx, w] = pwelch(X,WINDOW,NOVERLAP,F,FS)

     where Pxx is the PSD estimate, w is frequency (normalized, in case we set the FS value it is the actual frequency), WINDOW window function (for example hamming(100), the length of the window is the length of the signal segment), NOVERLAP number of overlapping samples, F is a vector of frequencies for which we want to estimate PSD (can be left empty [] and the function sets the frequencies according to FFT), FS is sampling frequency.

Useful functions

1. alpha-masking can be used To visualize the segmentation as a semi-transparent layer

   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)