Global Tracking of Myocardial Motion in Ultrasound Sequence Images: A Feasibility Study

1.   Introduction

Cardiovascular diseases (CVDs) threaten human's health and lives severely. The diagnostic and prognostic mechanical changes of myocardium caused by CVDs could be better investigated by assessing myocardial motion ^[1]. Mangner et al. reported that obese children had unfavorable alterations in myocardial geometry and function compared with nonobese children ^[2]. The assessment of myocardial shortening that occurs after end-systole (post-systolic shortening, PSS) is greatly useful for the diagnosis of acute cardiac ischaemia ^[3] and have a great potential in providing prognostic information in patients with acute coronary syndrome (ACS) following percutaneous coronary intervention ^[4]. Huang et al. discovered that different patterns of layer-specific strain distribution of myocardium in left ventricle (LV) might result from different etiologies of myocardial hypertrophy such as aortic stenosis (AS) and hypertrophic cardiomyopathy ^[5]. Shiino et al. suggested that multilayer global longitudinal strain (GLS) could be used to detect the early improvement in LV systolic function in patients with transcatheter aortic valve implantation ^[6].

In order to evaluate the cardiac mechanical changes caused by CVDs, many imaging techniques have developed. Echocardiogram is widely used in clinic to evaluate the myocardial function due to its real-time feedback, noninvasiveness and portability. This approach provides M-mode ultrasound image but focuses on the evaluation of systolic function of myocardial tissues at the targeted scan line. A color Doppler velocity imaging was developed by McDicken et al. to estimate the myocardial motion ^[7]. However, the Doppler-based methods mainly estimate the velocity in axial direction, which is parallel to the ultrasound beam. It is difficult to obtain accurate evaluation of myocardial motion. Based on the premises that malignant masses are stiffer than benign masses and normal tissue, ultrasound elastography was proposed to quantitatively evaluate the kinematic parameters of soft tissue in 1991 by Ophir et al. ^[8]. Elastography is a useful technique to obtain the response (i.e. deformation) of tissues to an external applied compression. Different from other static organs, heart is a muscular organ that contracts rhythmically. Therefore, myocardial elastography (ME) tracking the deformation of the myocardium was proposed by Konofagou et al. in 2002 ^[9]. The crucial step of ME is the estimation of motion or time-delay.

Numerous approaches have been developed for motion estimation such as block matching ^{[10,11,12,13,14,15,16,17,18,19]}, optical flow ^[20], nonrigid registration ^[21] and dynamic programming ^[22] methods. Block tracking methods are the dominant methods used to calculate the displacement. The main idea of block tracking algorithm is to calculate similarity metrics including sum of absolute (SAD) ^[10], sum of squared differences (SSD) ^[11,12], cross correlation (CC), and normalized cross correlation (NCC) ^{[9,13,14,15,16]}. NCC algorithm has been commonly regarded as the gold standard of motion estimation ^[11] and successfully used to estimate the mechanical property of the tissue such as breast ^[14], myocardium ^[9,13,15,16]. However, due to the overlap of the successive searching blocks and the exhaustive search algorithm, the computational complexity of the motion estimation remains a crucial drawback in the use of NCC. This problem becomes severer when the ultrasound radio-frequency (RF) signals are used as source data. RF frame data have more information including magnitude and phase of the signals in comparison with B-mode grey scale image. Therefore, the RF signal-based methods display higher accuracy in motion estimation. However, the large amount of RF signals increases the computational complexity of the motion estimation.

Therefore, many methods have been developed to reduce the computational complexity of motion estimation using NCC. One strategy is the use of sum tables. According to the equation of NCC, sum tables were constructed in the purpose of avoiding repeated calculation of NCC between each pair of matching and candidate blocks. The sum-table-based NCC could be a very useful and flexible real-time motion estimator ^[15]. A coarse-to-fine strategy was employed to develop a two-step approach for robust estimation of displacement in real-time strain imaging for freehand ultrasound elastography ^[18]. Another useful way to reduce the computational cost is to use a guide-search strategy in the calculation of NCC instead of exhaustive search ^[23]. However, that strategy tends to be more suitable to estimate the motion of static organs such as breast. As heart is a dynamic organ with intrinsic contractility, an exhaustive search is necessary due to myocardial activity. Therefore, a fast motion estimation is required for ME. The full-field motion map of the displacements should be obtained block by block using the NCC method. It means the similar displacement calculation procedure has to repeat for all the blocks. In order to reduce the computational cost, a non-NCC based global time-delay estimation algorithm called GLobal Ultrasound Elastography (GLUE) was established by optimizing a nonlinear cost function formulated by values of all samples of RF data, and was verified to be a novel technique for calculating both axial and lateral displacement fields between two frames of RF data ^[24]. However, GLUE was more suitable for the estimation of static soft tissues for example liver, while the feasibility and accuracy in evaluating the motion of inherently dynamic organs for example heart remains unverified. ME is proposed as a cross correlation based method to assess the local and global mechanical functions of myocardium ^[9]. NCC methods had been successfully used in many myocardial mechanical studies ^[9,13,15,16].

ME provides mechanical property of myocardium to district normal and pathological tissues. However, there is no anatomic structure of heart clearly shown in the estimated displacement map. An automatic tracking algorithm was proposed by Luo et al. to track the myocardial contour of the entire LV throughout a sequence images of heart during a full cardiac cycle ^[25]. Therefore, the pseudo-color coded displacement map of the target tissues fuses with B-mode ultrasound gray-scale image for a better visualization of ME imaging.

This study proposes a novel framework of the motion estimation to evaluate the dynamic mechanical property of myocardium, including displacement estimation, myocardium extraction, and image fusion. First, a NCC based global motion estimator is proposed as vectorized normalized cross-correlation (VNCC) to estimate the displacement values of the myocardial tissue. Then, the myocardium is semi-automatically extracted in a sequence of ultrasonic images based on the displacement values. Finally, an improved fast mapping algorithm (IFMA) is developed to obtain the fusion image in which the anatomic structure and the mechanical information are displayed simultaneously. In this study, RF signals are acquired from murine heart in left ventricular parasternal long-axis (PLAX) view and used to verify the feasibility of the proposed framework.

2.   Materials and methods

2.1.   Data acquisition

This study selected one nude mouse and five C-57 mice as experimental subjects for a feasibility study on dynamically global tracking of myocardial motion in ultrasound sequence images. The animal was first anesthetized with 40 mg/kg intraperitoneal injection of pentobarbital sodium (Solarbio Inc., Beijing, China), and then its limbs were fixed to a physiological monitoring platform in supine position. The platform is part of the ultrasound imaging system (Vevo2100, Visual Sonics Inc., Canada), which can monitor electrocardiogram (ECG) and remain the body's temperature when scanning. Hair at the mouse's chest was removed using depilatory cream. A probe (MS-550S, Visual Sonics Inc., Canada) with a center frequency of 40 MHz, a geometric focus of 6mm and a bandwidth of 32–55 MHz was used to acquire the PLAX view of heart of nude mouse. A probe (MS-400, Visual Sonics Inc., Canada) with a lower center frequency of 30 MHz, a geometric focus of 7 mm and a bandwidth of 22–55 MHz was used to acquire the PLAX view of heart of C-57 mouse. All experiments were performed according to the Guidelines for the Care of Laboratory Animals of the National Institutes of Health. This study was approved by the Institutional Ethics Committee at Southern Medical University.

In this study, 32 bit In-phase quadrature (IQ) data were exported for off-line data processing. There are 300 frames that were recorded at a sample rate of 128 MHz. Then, RF data were reconstructed from IQ data. For the RF frame acquired using MS-550S transducer, the size of one pixel in the axial and lateral direction is 0.02 and 0.055 mm, respectively. For the RF frame acquired using MS-400 transducer, the size of one pixel in the axial and lateral direction is 0.02 and 0.06 mm, respectively.

2.2.   Data processing

A total number of 39 frames recorded in a whole cardiac cycle were processed in this study. Figure 1 shows the framework of myocardial motion estimation in this study. The RF data as the source data were used to reconstruct gray-scale B-mode images using Hilbert transformation and log expression for myocardium extraction and image fusion procedures. A vectorized normalized cross-correlation (VNCC) method was proposed to calculate the displacement distribution between each two adjacent RF frames. Then, the area of myocardium in the sequences of B-mode images were semi-automatically extracted. Finally, an improved fast mapping algorithm (IFMA) was proposed to obtain the fusion images, in which the anatomic structure and mechanical information of the myocardium can be dynamically displayed simultaneously.

Figure 1.  Flowchart of the framework of myocardial motion estimation. The framework includes B-mode image reconstruction from RF data, displacement estimation of two successive RF frames, myocardium extraction, and image fusion. In B-mode ultrasound image, A = apex, C = cavity, AW = anterior wall, PW = posterior wall, AO = Aorta.

DownLoad: Full-Size Img PowerPoint

2.2.1.   B-mode image reconstruction

B-mode image was reconstructed from RF data. Hilbert transformation and log compression were performed on RF data (Eq (1)):

where ${\mathit{\boldsymbol{Bmod}}}{{\mathit{\boldsymbol{e}}}_{temp}}$ in fractional form represents amplitude values of RF signals. Then, ${\mathit{\boldsymbol{Bmod}}}{{\mathit{\boldsymbol{e}}}_{temp}}$ values were quantized to integers between 0 and 255 to be brightness values in gray-scale B-mode ultrasound image.

2.2.2.   Displacement estimation

(1) Normalized Cross-Correlation (NCC)

As shown in Figure 2, ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ and ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$ represent two adjacent RF frames of PLAX view acquired from ultrasound scanning. Assume that the sizes of all RF frames are $m \times n$. The NCC value between the matching block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ and a candidate block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$ is commonly calculated according to Eq (2) ^[26].

(2)

Figure 2.  Diagram of displacement estimation using NCC. Yellow block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ is matching block. NCC values are calculated between matching block and candidate blocks (gray areas) in a search range (green area) defined in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$. The shifting values ${\tau _x}$ and ${\tau _y}$ of candidate block with maximum NCC are the axial and lateral displacement of the matching block, respectively.

DownLoad: Full-Size Img PowerPoint

where $u$ and $v$ represent the location of sample of interest in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$, $0 \leqslant u \leqslant n$ and $0 \leqslant v \leqslant m$. ${\tau _x}$ and ${\tau _y}$ are the shifts between the matching block and the candidate block, $ - a \leqslant {\tau _x} \leqslant a$ ($a \geqslant 0$) and $ - b \leqslant {\tau _y} \leqslant b$ ($b \geqslant 0$). $a$ and $b$ define the search range for candidate blocks in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$ in axial direction and lateral direction, respectively. The positive value of ${\tau _x}$ represents a downwards movement, while the negative value means an upwards movement. The positive value of ${\tau _y}$ represents a rightwards movement, while the negative value means a leftwards movement. The size of all the blocks is defined as $({B_x} + 1) \times ({B_y} + 1)$. $\overline {{f_{xy}}} $ and $\overline {{g_{xy}}} $ are the average values of the samples in the region of matching and candidate block, respectively. ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}(x, y)$ and ${{\mathit{\boldsymbol{F}}}_2}(x, y)$ represent the value of RF sample located at ${x^{th}}$ row and ${y^{th}}$ column in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ and ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$, respectively.

In the process of calculating displacements of one matching block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$, NCC values are calculated between matching block and candidate blocks in a search range in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$. The candidate block with the maximum NCC value is regarded to have the highest similarity to the matching block. The shifts ${\tau _x}$ and ${\tau _y}$ between these two blocks represent the axial and lateral displacement, respectively (Figure 2). Then, the matching block shifts in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$, this process of displacement estimation repeats until the displacement distribution map of the entire frame ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ is obtained.

(2) Vectorized Normalized Cross-Correlation (VNCC)

VNCC proposed in this study is improved from NCC to simplify the calculation of NCC and then to speed up the procedure of motion estimation using NCC. Rather than one NCC calculation of each block using traditional NCC approach, the NCC values of all samples/blocks in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ could be calculated simultaneously using VNCC.

(3) Vectorized Sum-Table (VST) for NCC Calculation

VNCC proposed in this study is improved from NCC to simplify the calculation of NCC and then to speed up the procedure of motion estimation using NCC. Rather than one NCC calculation of each block using traditional NCC approach, the NCC values of all samples/blocks in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ could be calculated simultaneously using VNCC.

(3) Vectorized Sum-Table (VST) for NCC Calculation

Firstly, Eq (2) of NCC is transformed to avoid calculation of $\overline {{f_{xy}}} $ and $\overline {{g_{xy}}} $. The denominator of Eq (2) is rewritten as two terms (Eqs (3) and (4)). The numerator of Eq (2) is calculated by Eq (5).

(5)

Then, the NCC value between the matching block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ and each candidate block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$ can be calculated only by five terms of summation, that is (1) the sum of matching block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$:$\sum\nolimits_{x = u - \frac{{{B_x}}}{2}}^{u + \frac{{{B_x}}}{2}} {\sum\nolimits_{y = v - \frac{{{B_y}}}{2}}^{v + \frac{{{B_y}}}{2}} {{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}(x, y)} } $, (2) the energy of matching block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$: $\sum\nolimits_{x = u - \frac{{{B_x}}}{2}}^{u + \frac{{{B_x}}}{2}} {\sum\nolimits_{y = v - \frac{{{B_y}}}{2}}^{v + \frac{{{B_y}}}{2}} {{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}^{\mathit{\boldsymbol{2}}}(x, y)} } $, (3) the sum of candidate block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$: $\sum\nolimits_{x = u - \frac{{{B_x}}}{2}}^{u + \frac{{{B_x}}}{2}} {\sum\nolimits_{y = v - \frac{{{B_y}}}{2}}^{v + \frac{{{B_y}}}{2}} {{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}(x + {\tau _x}, y + {\tau _y}} })$, (4) the energy of candidate block in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$: $\sum\nolimits_{x = u - \frac{{{B_x}}}{2}}^{u + \frac{{{B_x}}}{2}} {\sum\nolimits_{y = v - \frac{{{B_y}}}{2}}^{v + \frac{{{B_y}}}{2}} {{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}^{\mathit{\boldsymbol{2}}}(x + {\tau _x}, y + {\tau _y}} })$, and (5) the Cross-Correlation (CC) calculation between the two blocks: .

For further reduction of the computational complexity, NCC values between all the matching blocks in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ and its corresponding candidate blocks in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}$ are calculated simultaneously by proposing Vectorized Sum-Table (VST) for the five terms (Eq (6)):

(6)

where $u$, $v$, ${\tau _x}$ and ${\tau _y}$ are defined as same as in Eq (2). ${\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}^{\mathit{\boldsymbol{2}}}(x, y)$ represents the squared value of RF sample located at ${x^{th}}$ row and ${y^{th}}$ column in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$, while ${\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}^{\mathit{\boldsymbol{2}}}(x + {\tau _x}, y + {\tau _y})$ represents the squared value of RF sample located at ${(x + {\tau _x})^{th}}$ row and ${(y + {\tau _y})^{th}}$ column in ${{\mathit{\boldsymbol{F}}}_2}$. ${\mathit{\boldsymbol{VS}}}{{\mathit{\boldsymbol{T}}}_{{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}}}(u, v)$ and ${\mathit{\boldsymbol{VS}}}{{\mathit{\boldsymbol{T}}}_{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}^{\mathit{\boldsymbol{2}}}}}(u, v)$ represent summation of values and squared values of samples in matching block with a center point located at ${u^{th}}$ row and ${v^{th}}$ column in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ respectively, ${\mathit{\boldsymbol{VS}}}{{\mathit{\boldsymbol{T}}}_{{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}}}(u, v)$ and ${\mathit{\boldsymbol{VS}}}{{\mathit{\boldsymbol{T}}}_{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}^{\mathit{\boldsymbol{2}}}}}(u, v)$ represent summation of values and squared values of samples in candidate block with a center point located at ${(u + {\tau _x})^{th}}$ row and ${(v + {\tau _y})^{th}}$ column in ${{\mathit{\boldsymbol{F}}}_2}$ respectively, ${\mathit{\boldsymbol{VS}}}{{\mathit{\boldsymbol{T}}}_{{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}{\mathit{\boldsymbol{, }}}{{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{2}}}}}}(u, v)$ represents the sum of Hadamard product of matching and candidate block described previously.

Finally, from Eq (6), the VNCC values of the entire RF frame are formulated as Eq (7), where $ * $ represents the Hadamard product of two matrices.

(4) Displacement Calculation Using VNCC

The displacement distribution in ${{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{1}}}}$ could be obtained by Eqs (6) and (7), Initially ${\tau _x} = - a$ and ${\tau _y} = - b$. With increasing ${\tau _x}$ and ${\tau _y}$ by 1 to their upper bound ($a$ and $b$, respectively), the current matrix ${\mathit{\boldsymbol{VNC}}}{{\mathit{\boldsymbol{C}}}_{m \times n}}$ was obtained and compared with the previously obtained element by element. The ${\mathit{\boldsymbol{VNC}}}{{\mathit{\boldsymbol{C}}}_{m \times n}}$ element values and the values of ${\tau _x}$ and ${\tau _y}$ were updated with the larger current ${\mathit{\boldsymbol{VNC}}}{{\mathit{\boldsymbol{C}}}_{m \times n}}$ values and its corresponding values of ${\tau _x}$ and ${\tau _y}$, respectively. After ${\mathit{\boldsymbol{VNC}}}{{\mathit{\boldsymbol{C}}}_{m \times n}}$ was calculated in the whole search range, displacement matrices could be obtained. A median filter was employed to reduce the noise in displacement map. This study only focus on the axial displacement due to the lower accuracy of lateral displacement.

To validate the proposed algorithm, a sum-table based normalized cross-correlation (STNCC) calculation method ^[15] is employed as reference algorithm in this study. Using VNCC and STNCC, the computation time of displacement calculation and the global tracking of myocardial motion of anterior wall (AW) and posterior wall (PW) during one cardiac cycle are calculated and compared.

2.2.3.   Myocardium extraction

In this study, the region of interest (ROI), i.e. the myocardial tissue, was manually extracted in the first frame of the B-mode image sequence. The ROI in the next frame was updated according to the corresponding axial displacement values. A mask image was created by assigning 1 and 0 to the ROI and background area, respectively. The mask image was processed by image dilation and erosion to reduce the noise caused by myocardial movement or myocardium thickening/thinning. Then the aforementioned process was repeated for all the images in the sequence. The myocardial tissues were automatically extracted for displaying myocardial motion in the ultrasound sequence images.

2.2.4.   Image fusion

To clearly display both anatomic and mechanical properties of myocardium, a weight factor $p$ $(0 < p < 1)$ was introduced to make the pseudo-color-coded ROI (${\mathit{\boldsymbol{I}}}$, obtained by a image fusion method of FMA ADDIN EN.CITE.DATA ^[27]) transparently overlapped on the gray-scale B-mode image (${{\mathit{\boldsymbol{I}}}_{BMODE}}$). The final fusion image ${{\mathit{\boldsymbol{I}}}_{final}}$ was obtained by Eq (8).

In the process of displacement calculation, the block size was set to $27 \times 9$ samples and the search range was $ \pm 40$ samples in axial direction and $ \pm 4$ RF lines in lateral direction. To compare the global tracking of myocardial motion in ultrasound sequence images using VNCC with STNCC, Matlab (version 9.0, R2016b, Mathworks Inc., Natick, MA) was used to program all the codes of image reconstruction, displacement estimation, myocardium extraction, and image fusion on a PC workstation (Intel Core i5-6500 CPU, 3.2 GHz, 8GB RAM).

3.   Results

3.1.   Validation results

The deformation of myocardium during one cardiac cycle was similarly tracked by VNCC as that of STNCC. The positive and negative values of displacement indicate myocardium moving downwards and upwards, respectively. As shown in Figure 3, myocardial motion of AW and PW could be tracked using both VNCC and STNCC. Compared with AW, the motion of PW was clearly shown in an opposite direction. The averaged absolute error in displacement between the two methods ranges from 1 to 3 μm.

Figure 3.  The average displacements of AW and PW calculated in one cardiac cycle using VNCC and STNCC. (a, b) results of nude mouse heart; (c, d) typical results of one C-57 mouse heart. ∆ represents the averaged absolute error in displacement values between the two methods.

DownLoad: Full-Size Img PowerPoint

Meanwhile, in comparison with STNCC, the computation time reduced approximate 10s using VNCC (Table 1). However, more efforts are needed for the application of VNCC in real-time myocardial elastography. It was noticed that the computation time of displacement estimation for the nude mouse heart is longer than that of C-57 mice using both VNCC and STNCC. This is because of the larger size of RF frame of the nude mouse heart scanned by a higher frequency ultrasound transducer.

Table 1.  Computational time (s) of VNCC and STNCC.

	VNCC	STNCC
nude mouse	73.26 ± 1.34	85.17 ± 0.78
C-57 mice	55.73 ± 1.52	64.88 ± 1.64

---

 | Show Table

DownLoad: CSV

3.2.   Improvements in ME imaging

In this study, the myocardial tissue could be dynamically extracted with high quality. It was noticed that the myocardium could not be extracted smoothly if the ROI of the myocardial tissue was updated only based on the corresponding axial displacement values. As shown in Figure 4a, some speckles are found in the resulted myocardium area indicating that the myocardial tissues in those regions are not extracted properly. The results of myocardium extraction improved greatly with process of image dilation and erosion (Figure 4b). The myocardial tissues in the sequence images were automatically extracted. These results provided a good base for motion estimation of myocardium.

Figure 4.  The resulted mask images of myocardial extraction with (b) and without (a) image dilation and erosion. Red arrows in (a) indicate the myocardial tissues that are not extracted.

DownLoad: Full-Size Img PowerPoint

The quality of the fused images improved greatly with the application of IFMA (Figure 5). The pseudo-color-coded ROI transparently overlapped on the gray-scale B-mode image with introduction of the weight factor. This resulted in clear visualization of both anatomic and mechanical properties of myocardium, which will benefit the analysis of the abnormal motion of myocardium.

Figure 5.  The fused images using fast mapping algorithm (FMA) (a) and improved fast mapping algorithm (IFMA) (b).

DownLoad: Full-Size Img PowerPoint

3.3.   Estimation of myocardial motion

With VNCC, the displacement map of myocardium was obtained from tracking the shift between two adjacent RF frames of PLAX view.

In this study, the frames were chosen according to the ventricular systolic phase of ECG. The displacement map could clearly show the amplitude and direction of myocardial movement during contraction period. As shown in Figure 6, the AW moves downwards while the PW moves upwards showing that the heart was contracting. Furthermore, the motion of the myocardium near endocardium was larger than that of the myocardium near epicardium. The similar results were obtained for C-57 mice.

Figure 6.  The axial displacement map at one moment in ventricular systolic phases (red spot labelled on ECG signal) of nude mouse. Positive (red) and negative (blue) values indicate myocardium moving downwards and upwards, respectively. The color bar shows the amplitude of displacement. A, C, AW, PW and AO represent as same anatomic sites of left ventricle as those in Figure 1.

DownLoad: Full-Size Img PowerPoint

This study dynamically estimated the movement of myocardium during cardiac cycle. Figure 7 shows three ME images at begin, middle and end of the ventricular systolic phases of nude mouse. The myocardial motion was visualized in the sequence ultrasound images. It was found that the myocardium grown thicker with the contraction of left ventricle. Additionally, the area of myocardium with the largest movement propagated from apex to base of the posterior wall, and a reversed trend (from base to apex) was found in anterior wall. The obtained ME images visualized the morphological and mechanical changes during the contraction period of left ventricle, which could estimate the contraction function of the myocardium in the systolic phase. The similar results were obtained for C-57 mice.

Figure 7.  Axial displacement maps obtained by the proposed VNCC algorithm at begin (a), middle (b) and end (c) of the ventricular systolic phases of nude mouse heart. Positive (red) and negative (blue) values in color bar represent moving upwards and downwards, respectively.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

This study introduced a novel framework for global tracking of myocardial motion in ultrasound sequence images, which consists of four parts: B-mode image reconstruction from RF data, global displacement estimation, myocardium extraction and image fusion. The algorithm of global displacement estimation plays an important role in myocardial motion tracking.

An exhaustive search and high overlap between blocks were essential to maintain the high-quality motion estimation of myocardium. Under this condition, the proposed VNCC reduced computational time by simultaneously calculating the displacements of all samples in a RF frame. A previous study proposed sum-table-based NCC approach to reduce the computational complexity of NCC, and had been verified by calculating the displacements of a human abdominal aorta and a human left-ventricle in a parasternal long-axis view, but this approach is still unable to calculate all samples simultaneously ^[15]. Dynamic programming (DP) approaches ^[28] were also proposed and had been verified in static soft tissues. However, the accuracy of myocardial motion estimation using DP still needs to be verified.

The extraction of myocardium from background tissue was necessary for the better visualization and interpretation of myocardial motion. This study proposed a semi-supervised automatic extraction method, which is automatic extraction in the sequence images based on the manual segmentation of the ROI in the first frame. Similar methods had been proposed in other study ^[25], but they focused to extract the contours of ROI instead of the region of ROI and interpolation was employed to overcome the discontinuity of the traced contours. Similarly, in our method, the speckles (discontinuity) were found in the extracted myocardial tissue, which probably were caused by myocardial movement or myocardium thickening/thinning. We applied image dilation and erosion to obtain the intact target tissue.

The effective tracking of myocardial motion in the ultrasound sequence images could be useful to evaluate the systolic function of myocardium. As shown in Figure 5, the motion of the myocardium close to endocardium is larger than that at epicardium. This phenomenon is consistent with previous study on transmural heterogeneity of myocardium ^[29]. Furthermore, previous study reported that the mechanical activation sequence of heart was physiologically along the direction of apex-to-base ^[30]. In this study, cardiac contractility was dynamically visualized in Figure 6, showing a similar propagation of contractility of myocardium at posterior wall. This finding was confirmed by the results reported in a recent study ^[31]. However, a base-to-apex propagation was also found in anterior wall in this study. The possible reason might be the complexity of heart movements, including myocardial rotation, systole, and diastole. The robustness of VNCC was demonstrated by the effective myocardial motion tracking of different mouse species with different transducers.

The proposed VNCC method in this study achieved a high processing speed at expense of memory space consumption. 2-D matrixes i.e. vectorized sum-tables were established to calculate the NCC values. This method resulted in the decrease of computational complexity. The large memory space consumption appears to be acceptable due to advanced equipment of large computer workstation.

Although this study could efficiently and globally calculate the myocardial displacement using VNCC and informatively present the mechanical property and anatomical structure of myocardium, there are also limitations of this study. First, the sub-sample displacement of each sample was not estimated. In order to compensate this weakness, a median filter was simply employed in this study. A novel method for sub-sample displacement evaluation should be emphasized in the future work to improve precision ^[32] and continuity of the myocardial motion. Secondly, this study only verified the feasibility of the proposed framework using RF data of murine heart. The mechanical difference between healthy and pathological myocardium will be quantitatively evaluated in future animal experiments. Thirdly, myocardial motion was only estimated in two dimensional plane in this study. However, the movement of the heart is three dimensional. It can lead to artifacts in ME images due to the out-of-plane motion like elevational tissue strains ^[33]. The development of three dimensional (3D) ultrasound technique ^[34,35] lays a solid foundation in the performance of 3D myocardial elastography including the estimation of out-of-plane motion. 3D motion estimation of the heart will be conducted in future work.

5.   Conclusion

In this study, a novel framework of global tracking of myocardial motion was proposed to visualize contractility of myocardium in ultrasound sequence images in ventricular systole. To speed up the processing speed, a VNCC approach was proposed to globally estimate the displacements. An improved fast mapping algorithm (IFMA) was used to display the bimodal fusion images. The results demonstrate that the proposed framework is an efficient tool for the estimation of myocardial motion reflecting systolic function of murine heart. The visualization of myocardial motion is useful to evaluate cardiac function and will benefit diagnosis and prognosis of CVDs.

Acknowledgments

This work is supported by National Natural Science Foundation of China (No. 81371560), Key R & D Programs in Guangdong (No. 2018B030333001), and Guangdong Science and Technology Program (No. 2016A020216017 and No. 2013B021800039).

Conflict of interest

All authors declare no conflicts of interest in this paper.