|Home | About | Journals | Submit | Contact Us | Français|
Intravascular ultrasound elastography (IVUSe) could improve the diagnosis of cardiovascular disease by revealing vulnerable plaques through their mechanical tissue properties. To improve the performance of IVUSe, we developed and implemented a non-rigid image-registration method to visualize the radial and circumferential component of strain within vascular tissues. We evaluated the algorithm’s performance with four initialization schemes using simulated and experimentally acquired ultrasound images. Applying the registration method to radio-frequency (RF) echo frames improved the accuracy of displacements compared to when B-mode images were employed. However, strain elastograms measured from RF echo frames produce erroneous results when both the zero-initialization method and the mesh-refinement scheme were employed. For most strain levels, the cross-correlation-initialization method produced the best performance. The simulation study predicted that elastograms obtained from vessels with average strains in the range of 3%–5% should have high elastographic signal-to-noise ratio (SNRe)–on the order of 4.5 and 7.5 for the radial and circumferential components of strain, respectively. The preliminary in vivo validation study (phantom and an atherosclerotic rabbit) demonstrated that the non-rigid registration method could produce useful radial and circumferential strain elastograms under realistic physiologic conditions. The results of this investigation were sufficiently encouraging to warrant a more comprehensive in vivo validation.
Intravascular ultrasound (IVUS) devices offer a minimally invasive, high-resolution imaging technique to visualize blood vessels from within the vessel lumen. IVUS provides real-time two-dimensional (2-D) cross-sectional images of the arterial walls. IVUS imaging has been used to improve the diagnosis of cardiovascular pathologies, such as atherosclerosis, by identifying arterial wall morphology, fibrous plaque locations and thicknesses, and internal lumen area (Gussenhoven et al. 1989). However, the use of IVUS for diagnosing vulnerable plaques has been shown to have low sensitivity in identifying specific plaque components that may indicate a propensity to rupture, such as lipid-filled plaques or fibrous caps (Palmer et al. 1999). In addition to plaque components, mechanical markers, such as circumferential stress, have been suggested as possible indicators of a plaque’s likelihood to rupture (Lendon et al. 1991; Loree et al. 1992).
Intravascular ultrasound elastography (IVUSe) visualizes the mechanical properties of vascular tissue. It has been suggested that IVUSe may provide additional markers that can improve the prediction of plaque rupture. These markers include peak strain, mechanical modulus and peak stress (Baldewsing et al. 2004, 2005; de Korte et al. 2000; Maurice et al. 2008a; Richards and Doyley 2011; Wan et al. 2001).
The general principles of elastography can be described as follows. First, a particular tissue or organ is deformed under an applied force, either externally or internally. Second, the tissue is imaged multiple times as it deforms and the images at different times are compared to provide a relative measure of deformation. Lastly, the deformation of the tissue is used to infer the mechanical properties of the tissue, either through the creation of strain images or through further processing steps. For intravascular ultrasound elastography, vascular tissue is imaged with an IVUS catheter transducer as the vessel expands and contracts with the internal cardiac pressure. Strain images are then calculated from the gradient of the radial displacement field (de Korte et al. 2000) or modulus elastograms are computed by solving the inverse elasticity problem (Baldewsing et al. 2006; Le Floc’h et al. 2009; Richards and Doyley 2011). IVUSe has been shown to accurately determine plaque composition and potentially improve the identification of rupture-prone plaques (de Korte et al. 2000, 2002; Fromageau et al. 2008; Liang et al. 2009). To obtain strain images with high signal-to-noise ratios (SNR) in elastographic image processing, an accurate estimation of tissue displacement is an essential step. Furthermore, to reconstruct the mechanical parameters of vascular tissue, accurate displacement vectors (with two or three components) are also necessary, to limit the “ill-posedness” of the inverse elasticity problem.
In this article, we hypothesized that applying a non-rigid image-registration (IR) method to radio-frequency (RF) echo frames would improve the performance of IVUSe. Furthermore, we suggest that the use of the RF image signal will improve the accuracy of the measurement over the use of B-mode images. However, we recognize that to register RF images, improved techniques for algorithm initialization must be considered. The registration algorithm developed here uses a bi-linear, finite-element-based discretization of the displacement field. The algorithm also includes a novel, spatially modulated regularization of the strain magnitudes.
To corroborate our hypothesis, the algorithm was tested on simulated IVUS ultrasound images created using finite-element modeled displacement fields. Using these images, a direct comparison of the accuracy of B-mode image registration and RF image registration was performed. We performed studies to determine the impact of several initialization schemes for RF image registration and employed a mesh-refinement scheme similar to that used by Liang et al. (2008). The results are presented and compared for both RF and B-mode measurements. We also performed an experiment on a vessel phantom to validate the simulation study. Our simulation studies showed that, indeed, registration of RF images outperformed that of B-mode images if the proper initialization techniques were used. Relatively high radial strain SNR were reported (4.9 ± 0.1) for average strain values of 1% (3.2% peak strain). Strain elastograms obtained from vessel phantoms suggested that use of RF images would lead to more accurate measurements.
The most prevalent feature-matching algorithm for measuring displacements in ultrasound (US) images is known as localized cross-correlation (CC) echo tracking (Chen et al. 1994; Hall et al. 2003, Lubinski et al. 1999; Ophir et al. 1999; Walker and Trahey 1995). This method is computationally inexpensive, easy to implement in real time (Hall et al. 2003) and can accurately estimate the component of the displacement in the direction of the ultrasound beam. Application of CC techniques in IVUS elastography (IVUSe) has been widely studied in vitro and in vivo (de Korte et al. 1997, 2002; Ryan and Foster 1997). Many groups have expanded on the ideas of localized CC techniques to further improve the performance of the axial, lateral and, in limited cases, the elevational displacement estimates. For example, temporal stretching has been used to re-correlate the RF echo frames prior to performing the standard correlation technique to reduce signal decorrelation (Alam and Ophir 1997; Cespedes and Ophir 1993). This technique uses a constant strain value for the entire image in the beam (or axial) direction only, which assumes that the tissue deformation is only axial and is constant within the image. However, the tissue itself deforms axially, laterally and in shear, and is rarely constant in the images of interest. For IVUSe, the assumption of constant radial strain is impractical because of the radial decay of the strain in a deforming vessel.Chaturvedi et al. (1998) used a similar global stretching of the post-compression image in two dimensions to reduce decorrelation noise attributable to the lateral tissue motion and, thus, improved the displacement variance. The stretching of the windows in this article was both axial and lateral; however, the lateral stretching was determined from an estimated axial stretching and the assumption of tissue incompressibility (i.e., a Poisson’s ratio of 0.5).Konofagou et al. (1998) developed a similar approach, in which both global and local stretching was performed, in conjunction with a weighted lateral interpolation, to improve lateral displacement estimates. This was done to overcome the poor precision of lateral displacement estimates attributable to low lateral sampling. Maurice et al. (2005, 2008b) developed the Lagrangian speckle model estimator for use in IVUSe. Here, a 2-D CC-based rigid registration was used to compensate for the constant shifts of window regions, both radially and angularly. Then, the strain was calculated in this window via an optimization algorithm designed to best fit the speckle of a pre-deformation image and a post-deformation image deformed by affine transformation. Here, the strains were measured intrinsically at each window location, rather than in postprocessing as a displacement gradient.
Several groups have investigated the use of regularization to penalize unrealistically large spatial variations of the displacements. Pellot-Barakat et al. (2004, 2006) developed a regularized optical flow technique that limited the square of the finite difference of neighboring displacements in 2-D, which was used to improve the axial displacement estimates. A similar approach was used by Rivaz et al. (2008) in the context of dynamic programming to improve the axial strain estimates in a real-time evaluation. Brusseau et al. (2008) used a spatially varying regularization technique to strongly penalize high variances in the displacement and strain in areas of poor CC coefficient values, using a discrete strain-estimation algorithm.
Deformable meshes have also been used to improve the quality of elastographic measurements (Basarab et al. 2008; Yeung et al. 1998; Zhu et al. 1999). This type of displacement discretization allows for higher-order representations (linear and second-order variations) of the local motion, which intrinsically compensate for image decorrelation in two and three dimensions. Unlike the previous methods, however, higher-order deformable meshes require an initial guess of the displacement field. The accuracy of the initialization can have a significant effect on the measurement accuracy.
Liang et al. (2008) have developed an image-registration based technique specifically for IVUS B-mode images of arteries. Their technique also uses a deformable mesh and cubic B-spline basis functions as the underlying displacement fields. The discretization of cubic B-splines, although common in many medical image-registration problems, implicitly imposes a smoothness constraint on the strain values. This may lead to inaccuracies in the measured strains if material property boundaries are present, leading to discontinuities in the strains. In addition, no work has been done to compare measurement accuracy using third-order displacement interpolation compared with second-order. Their algorithm includes regularization that penalizes the second-order spatial derivatives of the measured displacement field. Liang et al. (2008) also use a multi-resolution mesh-refinement strategy to avoid issues of local trapping or local minima. However, it is unclear as to whether such an initialization strategy would continue to work at high strain levels or if it is feasible for registration on the higher-frequency RF images.
In this work, we modeled the vascular tissue as a continuous material body with a deformation that can be described as:
Here, xo is the initial position of the material, xt is the position of the material at some time t > 0 and u(xo, t) is the time-dependent displacement field. We also assume that changes in the US signal are attributable only to the motion of the underlying material being imaged (i.e., u(xo, t)). Thus, image sequences (I) of this material taken in time can be related via the approximation:
The displacement field acts as a nonlinear scaling of the position vector that defines the image intensities. Equation (2) allows us to create an intensity-matching cost function that can be minimized, using a nonlinear optimization technique, to find the underlying displacement (u(xo, t)). The optimization function that we employed, including the cost function and regularization, has the form:
Here, Ω is the image domain for which the displacement measurement was calculated, I1 is the undeformed reference image, and I2 is the post-deformed image. To simplify the above expression, we excluded the explicit time dependence of the images and the displacement. Thus, xo = x. The first term of eqn (3) is the image-matching functional, πI. The second term of eqn (3) is the regularization, πR. We use an H1 semi-norm regularization to penalize large strains (usym = (uT + u)/2). The “:” in πR denotes the matrix inner product and α is an algorithmic constant that determines the weighting of the regularization, relative to the image-matching functional. The spatial weighting function, G(x), was used to spatially modulate the strain regularization. A spatially varying regularization can be used in situations where the a priori expectation of the regularized function is not constant to improve measurement accuracy (Pogue et al. 1999). In this algorithm, we set this function to G(x) = (dct)4, where dct is the radial distance from the centroid of the vessel’s inner lumen. This form of the spatial modulation was chosen based on the a priori expectation that the vessel’s strain field should have an approximate (dct)−2 decay from the inner lumen, as discussed in Timoshenko and Goodier (1970).
The functional (eqn ) minimization was performed using a Gauss-Newton method. Using this method, updates were computed iteratively as follows:
For simplicity, the x dependence in eqn 4 is implied (i.e., u(x) = u and I2(x + u(x)) = I2(u)). Here, w is an arbitrary, admissible variation of u. The subscript “v” indicated on the displacement field denotes the iteration number. Thus uv represents the current solution of u, at iteration v, and δu is the update to uv such that:
In this algorithm, bi-linear interpolation functions, with four-node, quadrilateral finite elements, were used to approximate u, δu and w. A linear representation of the element displacements reduces the effect of image decorrelation by allowing the image warping within each finite element to compensate for strain. This is similar to the Lagrangian speckle estimator method (LSME) described in Maurice et al. (2008b). The LSME method minimizes intensity differences, subject to a localized affine transformation within an US image to assess tissue strains. In our work, the deformation within each element is bi-linear. While this model adds an increased complexity, this level of displacement interpolation is similar to those found in most finite-element-based deformation models (Hughes 1987) and more realistically matches the true deformation. In addition, the minimization for a given measurement window of the LSME method is independent of its neighboring measurement. The minimization is global and the displacements across each element are continuous, thus creating an implicit smoothness to the displacement field in the formulation.
The discretized measurements were defined at the nodal coordinates. A rectangle rule was used to integrate the image matching terms of eqn (4), where the number of integration points was greater than or equal to the number of image pixels within an element. The images were interpolated at each integration point using cubic Lagrange polynomials. A 3 × 3 point Gaussian integration of the regularization terms was used. Finite-element shape functions were used for the interpolation of the displacement field and its variants (u, w and δu).
We implemented the non-rigid registration method in Fortran 90 and compiled it on a 16-core Intel Xeon Server (Intel, Santa Clara, CA, USA) that was operating at 2.93 GHz under the Centros 5.6 (64-bits) operating system. We used a parallelized linear solver as described in Schenk et al. (2008) to solve eqn (4) and the nodal values of δuv. The registration algorithm typically converged with 20 iterations (approximately 11–13 s) depending on the initialization method employed. More specifically, registration performed using the mesh refinement initialization took up to four times longer to converge than the other initialization schemes.
To evaluate the performance of the non-rigid displacement estimator, we performed computer simulations with strain levels ranging from 0.1% to 10%. Radio-frequency (RF) echo frames were synthesized by combining a finite-element mechanical model with the Field II (Jensen 1991) acoustic model using a four-step process. First, the point spread function (PSF) of the simulated IVUS system was computed using Field II. The simulated IVUS system consisted of a 40 MHz single-element unfocused transducer that had a 30% fractional bandwidth. Second, the acoustic response of the pre-compressed vessel was simulated by randomly distributing point scatters in two dimensions between a concentric 1.5 mm inner radius and a 4 mm outer radius and convolving the simulated PSF with these scatters. Third, the acoustic response of the post-compressed vessel was computed by redistributing the point scatters of the pre-compressed tissue using displacements that had been computed using a commercially available finite-element (FE) package (Abaqus; Dassault Systemes, Velizy-villacaublay, France) and convolving the simulated PSF with these scatters.
The commercial software package (Abaqus; Dassault Systemes) used to create the scatterer deformations modeled the vessel as a 2-D, plane-strain section of a nearly incompressible (v = 0.495) cylindrical vessel. The vessel contained a background material (22.5 kPa) with a 1.5 mm radius central hole, mimicking a vessel lumen. Adjacent to the lumen hole, a soft region (7.5 kPa) was made to simulate a lipid plaque. To simulate infinite boundary conditions, we constructed a FE model with an outer diameter of 50 mm. The circumferential component of displacement was assumed to be zero on four points on the periphery of the simulated vessel. Motion was induced in the simulate vessel by applying a distributed pressure to the inner lumen. Displacements from a smaller radial distance (4 mm) were used in the simulation. We have demonstrated in a previously reported study (Richards and Doyley ) that displacements computed by solving the forward FE problem with infinite boundary conditions were comparable to those measured from calibrated vessel phantoms. An average deformation of 1.0% strain was induced within the vessel by applying a uniformly distributed pressure (305 Pa) on the inner lumen. We scaled the computed displacements to simulate vessels with an average radial strain of 0.1%, 0.5%, 1%, 2%, 4%, 6%, 8% and 10%. Radial strain decayed from its peak value with the inverse square of radial distance in vessels; therefore, the reader should realize that although average strains were used to quantify the deformation, all elastograms contained strains that were larger and smaller than the average strain. Therefore, unless otherwise stated, the percentage strain and strain level refers to the average radial strain within the vessel. Figure 1 shows the simulated vessel model and the applied boundary conditions. The dotted line in Figure 1 represents the region of tissue imaged by the simulated IVUS system.
To create 2-D IVUS images, we rotated the scatters between 0° and 360° in steps of 1.41° for each rotation (256 A-lines) and the PSF was convolved with the scatters. The acoustic signal was sampled at 200 MHz axially. Additive white, Gaussian noise was added to the RF images with SNR of 8 dB, 12 dB and 17 dB. The speed of sound and attenuation were assumed to be 1540 m/s and zero, respectively. B-mode images were computed from the RF signals by complex demodulation using the analytic signals as described in (Tuthill et al. 1988).
The registration was performed using a uniform mesh. The radial and angular grid spacings were 0.094 mm and 5.63° (0.098 radians), respectively. The optimal magnitude of the regularization constant (α) was dependent on the strain magnitude (εavg) and the image type used for the regularization (i.e., B-mode or RF) (see Appendix). For image registrations using B-mode images, the regularization constant was chosen based on the equation α = 1 × 1018/ε avg. For image registrations using the RF images, the regularization constant was chosen using α = 5 × 1018/ε avg.
As stated earlier, the optimization technique required an initial guess of the displacement field to begin the minimization process. To compare the performance of different initialization techniques, four independent methods of creating the initial displacement values were used:
The initialization schemes listed above were compared for registrations performed on both B-mode and RF images.
Experimental studies were conducted to evaluate the performance of the registration method under more realistic conditions. A vessel phantom (3.16 mm inner diameter by 14 mm long) was constructed, containing a soft, crescent-shaped inclusion constructed from polyvinyl alcohol (PVA), as described in Richards et al. (2011). The vessel wall was fabricated from 10% by weight PVA solution in a cylindrical mold. The vessel background material and plaque regions were subjected to five and three freeze-thaw cycles, respectively. After the thermal cycling, the phantom was removed from the mold and stored at room temperature in water.
All ultrasonic imaging was performed at room temperature (approximately 20°C). IVUS imaging was performed using a commercially available ILab™ (Boston Scientific/Scimed, Natick, MA, USA) intravascular ultrasound scanner that was equipped with a 40 MHz Atlantis Pro imaging catheter (Boston Scientific, Natick, MA, USA). We used a PCI bus data-acquisition card (Compuscope 14200-1GB; Gage Applied, Lockport, IL, USA) to stream RF echo frames from the IVUS scanner at the full frame rate (30 fps) to a high-performance computer workstation. All RF echo signals were digitized to 14 bits at 200 MHz. Figure 2a shows a schematic of the experimental setup used to acquire the IVUS images. The vessel was slowly pressurized using the syringe pump and 20 RF echo frames were acquired during pressurization.Figure 2b shows an example B-mode image of the vessel phantom in the undeformed state.
The registration mesh was constructed using a three-step process. First, the inner lumen of the vessel phantom was manually segmented (see Fig. 2b). Then the elements were evenly spaced radially, for a given angular position, between the segmented inner lumen and approximately 0.5 mm from the image’s outer radius, with respect to the catheter position. At each angular position, 24 radial elements were used, with an average element size of 0.11 mm. Elements were evenly spaced angularly, with respect to the catheter position, with angular element size of 5.63° (0.098 radians) for a total of 64 elements. The image center did not coincide with the vessel center; thus, the spatial modulation function was set to G = q4, where q was the distance between the image point and the vessel’s centroid, calculated from the segmented inner lumen (see Fig. 2b). Seven registration sets were performed, where the first image frame taken in the experimental protocol was registered to the following seven image frames. We computed the average radial strain from displacements measured with the cross-correlation echo tracking method, which we interpolated onto the registration mesh. The average strain between a pair of RF echo frames was approximately 2.0%. Both the B-mode images and corresponding RF images were registered at each measured frame, for a total of 14 image registrations.
In the experimental study, the 2-D block matching, cross-correlation-based tracking algorithm was used to compute the radial displacement component from the two IVUS image frames of interest. All correlation-based measurements were performed using a measurement window size of 0.385 mm radially and 14.1° angularly. Each measurement overlapped 80% radially and 50% angularly. The measured radial displacements were then linearly interpolated at the nodal locations of the registration mesh and used as the initial displacement field. The initial angular displacement component was assumed to be zero at all the nodal locations. We set the regularization constant to α = 1 × 1018/ε avg for all registrations using B-mode images. Similarly, when using the RF images, the regularization constant was α = 5 × 1018/εavg (see Appendix). In this case, ε avg was estimated from the initialized cross-correlation displacement values.
IVUS RF images of artificially induced, abdominal aortic atherosclerosis were previously collected using a similar protocol as described above. Atherosclerosis was induced in the aorta of a New Zealand White (NZW) rabbit (2.5 months old, weight 1.5 kg) using a combination of high cholesterol diet (2% Cholesterol Research Diets, New Brunswick, NJ, USA) and balloon de-endothelization. The rabbit was maintained on a high cholesterol diet for 2 weeks before the arterial wall was injured (from the femoral artery to the iliac bifurcation) with a 4F Fogarty catheter. The rabbit was maintained on a high cholesterol diet after balloon injury, and once the atherosclerosis was developed (approximately six week post injury), the rabbit was sedated with ketamine and xylazine and a sheath was introduced in the carotid artery to facilitate echo imaging. After imaging, the rabbit was sacrificed and opening the pericardium exposed the heart. A small incision was made in the right atrium for draining of blood and perfusion fluids after they traveled through the circulatory system, and another was made at the apex of the left ventricle. A cannula was inserted into the left ventricle with the saline flowing, then moved through the left ventricle into the aorta. Saline perfusion was followed with para-formaldehyde (1%–4%), then with 10% then 30% sucrose. We used anatomic landmarks (identified with angiography) and radio-opaque rulers to register the histologic cross-sections and the ultrasound images. All procedures were performed in accordance with the University’s Institutional Animal Care and Use Committee.
We captured a total for 125 RF echo frames from the scanner at full frame rate (30 fps) from multiple cardiac cycles. The negative derivative of the intracoronary pressure signal was used to select RF echo frames (n = 4 for each cycle) from the diastolic phase of the cardiac. The decision to limit the processing to these frames was based on observation that the contraction of the heart and, therefore, catheter motion is relatively small in this phase of the cardiac cycle. Additionally, in this phase of the cardiac cycle, the intraluminal pressure decays relatively slowly and approximately constantly compared with other regions (i.e., systole), which allowed us to cope with structural decorrelation noise.
We performed a series of simulation experiments and a phantom study to assess the performance of an image-registration-based, displacement-estimation technique for IVUS elastography. In the analysis of the measured displacement values, accuracy measurements were calculated based on the strain values. In this work, the linearized Cauchy strains were calculated (i.e., ε = usym = (uT + u)/2). Finite differences were used to find the displacement derivatives at the center of each finite element. The results are shown with respect to the catheter-centric, polar coordinate strains. The strains (εrr, εθθ) were calculated as a function of the radial (ur) and angular displacements (uθ) as follows
The true strains were calculated from the true displacements, as determined by first interpolating the displacement nodal values of the finite-element simulation mesh onto the registration mesh. Figure 3a and b show an example (at 2% strain) of the true, finite-element simulated radial and angular displacement fields, respectively, used to deform the point scatterers in the simulation. Figure 3c shows the corresponding, true radial strain image calculated from the radial displacement field shown in Figure 3a.
The elastographic signal-to-noise ratio (SNRe) performance metric was used to quantify the simulation results, which was calculated as follows:
In this equation, εt is a single component of the true mechanical response (strain or displacement) and εm is the corresponding component of the measured response. The SNRe was calculated separately for each component of the measured response. The index n and number N correspond to the element number and total number of elements, respectively. In this work, the true and measured strains were calculated from the true and measured displacement fields, respectively, using the bilinear, finite-element shape functions at the centroid of each element (i.e., finite difference strain calculations). The calculations of the SNRe measurements, at each strain level and each noise level, were repeated for five independent realizations of the electronic noise and then averaged.
We have presented the results of the phantom registration experiments as a quantification of the strain elastographic contrast-to-noise ratio (CNRe) in addition to the corresponding strain images. The contrast was determined by manually segmenting the region of the simulated plaque (see Fig. 2b). The magnitude of the CNRe was calculated as (Techavipoo and Varghese 2005):
Here s1 is the mean strain in the simulated plaque and s2 is the mean strain in the background region. σ1 is the standard deviation of the strain in the inclusion and σ2 is the standard deviation in the background region. The CNRe analysis requires accurate segmentation of lumen/intima and media/adventitia boundaries, which can be challenging when IVUS imaging is performed at 40 MHz during clinical or animal studies. This issue was not encountered in either the phantom or animal study; however, we plan to use semiautomated segmentation methods in future studies. A comprehensive review of the most promising IVUS segmentation approaches are discussed in (Katouzian et al. 2012).
In Figure 4 the SNRe computed from radial strain elastograms estimated with the cross-correlation echo tracking method is plotted as a function of the strain. This figure reveals the characteristic trade-off between electronic noise and structural decorrelation reported inVarghese and Ophir 1997. More specifically, electronic noise degrades the performance of radial strain elastograms obtained at small strain levels (<1%). Structural decorrelation reduces the performance of elastograms obtained at higher strain levels.
Figure 5a–d show the radial strain signal-to-noise ratio, plotted as a function of the applied strain (εavg) when the registration algorithm was applied to B-mode images. Figure 5e–h show the radial strain SNRe when registering RF images using initialization procedures (1)–(4), respectively. The first initialization scheme uses the true displacement values as the initial guess. Thus these results represent the upper limit of the accuracy of this technique, with respect to the initialization. The peak SNRe achieved when elastograms were computed by applying the registration procedure to B-mode images was similar to that obtained using a simple cross-correlation method. The quality of both sets of elastograms degraded with increasing applied strains, but the rate of degradation was much slower for elastograms computed with the registration method. However, the peak SNRe of elastograms computed by applying the registration method to RF images was twice that obtained using the cross-correlation method.
The performance of the registration method when it was applied to RF echo frames improved with increasing applied SNRe when initialization scheme 1 was employed. SNRe reached its peak value (~5) at 2% applied strain. At higher strains, the accuracy decreased to approximately 3.5 SNRe. The B-mode measurement SNRe showed a similar increase until about 2% strain, then a dip followed by a gradual increase to larger strains. Figure 5b and f show similar SNRe plots when initialization scheme 2 was employed. Here the B-mode measurements followed a similar pattern as the B-mode measurements using scheme 1 for smaller strains (<2%). However, for (>2%) higher strains, the accuracy of these measurements dropped much lower than that of measurements from scheme 1. Figure 5c and g show the resulting radial strain accuracy as a function of the applied strain images, using the B-mode and RF images, respectively, under initialization scheme 3. These plots show that, for both B-mode and RF images, the resulting SNRe was greatly improved over initialization scheme 2. No appreciable difference in the accuracy was seen for US images of varying electronic noise. Figure 5d and h show the resulting radial strain accuracy as a function of the applied strain images, using the B-mode and RF images, respectively, under initialization scheme 4.
Figure 6 shows the circumferential strain SNRe when B-mode images were used in the registration procedure, for the initialization procedures (1)–(4), respectively. Figure 6 displayed a similar trend to that observed in Figure 5. However, the absolute values of the circumferential strain SNRe were higher because the circumferential strain was calculated by summing the angular derivative of the angular displacements and the radial displacements (eqn ).
In Figures 7 and and88 the elastographic signal-to-noise ratio computed from the estimated radial and circumferential components of displacement—the output of the registration algorithm—are plotted as function of the applied strain. These figures revealed that initialization scheme 1 gave the best performance, and initialization method 3 produced more accurate displacement estimates (an order of magnitude) for strain levels from 1%–5% when the registration method was applied to RF echo data. We expected higher SNRe values, which suggests that the meshing and interpolation process used in the registration process introduced errors. The amplitude of the SNRe plot shown in Figure 7 was considerably higher (a factor of 50) than those observed in Figure 5, which was not surprising since differentiation is a noise amplifier. However, the amplitude of the SNRe plot shown in Figure 8 was lower than that observed inFigure 6. We have to examine eqn (6) to understand this behavior. More specifically, circumferential strain was computed by summing a noisy term (the derivative of the circumferential component of displacement) and less noisy term (radial component of displacements) that was more dominate. Since Figure 8 only shows the behavior of the noisy circumferential component of displacements with applied strain, it wasn’t surprising that SNRe values were lower.
Figure 9 shows representative examples of radial and circumferential strain elastograms obtained when the registration procedure was applied to noisy B-mode and RF echo frames (i.e., images SNR of 12 dB) at 2% strain. Figure 3c shows the true strain elastograms.Figure 10a shows an example of radial displacement profiles that was obtained from the displacements used to compute the elastograms shown in Figure 9a and c.Figure 10b shows the corresponding radial strain profiles.Figure 11 shows a representative example of radial strain profiles obtained from elastograms computed from B-mode and RF images with different the initialization schemes.
Figure 12a and b show the example of radial and angular displacements computed from RF echo frames obtained from the vessel phantom. Figure 12c shows a plot of CNRe as a function of the average strain. CNRe was computed from strain elastograms measured from B-mode and RF images. CNRe increased from a very low value, peaked at approximately 1% and stayed relatively level. The radial strain CNRe of the RF image measurements increased initially from a low value, and then steadily declined beyond approximately 0.25% strain. Figure 13a and b show the radial and circumferential strain images measured from RF images obtained the vessel phantom at an approximate 0.4% strain. Figure 13c and d show the radial and circumferential strain images, respectively, measured using the phantom RF images at an approximate 2.0% strain.
Figure 14a shows a representative example of a sonogram obtained from the rabbit aorta. The sonogram reveals that the plaque was eccentric. Figure 14b–d show examples of composite radial strain elastograms computed near end diastole for three different cardiac cycles. These elastograms were computed from the average of four elastograms. Each elastogram was computed by first applying the digitized RF echo frames to a cross- correlation echo tracking algorithm and the using the estimated radial displacements to initialize the non-rigid registration method. The elastograms obtained from different cardiac cycles contained some local variations. Nonetheless, a localized region of high strain was observed between 10 and 11 o’clock in all elastograms. Figure 14e shows the corresponding histology image obtained from the corresponding cross section of the aorta. The histologic report of the rabbit aorta revealed that the high strain region between 10 and 11 o’clock corresponded to a mural thrombus.
Medical imaging techniques, such as intravascular ultrasound, are the diagnostic tools used to reduce deaths associated with cardiovascular disease. IVUS elastography is a promising adjunct technique that may aid in identifying vulnerable plaques by offering information about their mechanical properties. As the process of measuring material deformation from images is inherently ill posed, much of the work in the field of elastography has focused on improving the accuracy of measurements of tissue-displacement fields and subsequent strain fields. This work presents an image-registration technique that incorporates, in a general framework, methods such as local 2-D companding and spatially modulated regularization designed to improve measurement accuracy. We have shown that the use of our algorithm and the high-frequency (40 MHz) RF IVUS image information further improves the accuracy of the measurements compared with the use of lower-frequency B-mode images. However, one drawback of such a technique is the need for an initialization scheme to avoid local minima, or local trapping, caused by the localized search required by optimization techniques. It is expected that because the images are so different in their frequency characteristics, any initialization scheme would behave differently depending on the image type and the magnitude of the underlying displacements.
A similar registration-based, displacement-measurement technique was developed in Liang et al. (2008) and further used in Liang et al. (2009). The algorithm developed here differs from their implementation in two distinct ways. First, the displacement interpolation used in this work was linear, rather than cubic. A finite-element-based, linear interpolation was chosen here to facilitate the incorporation of model constraints in any future work. In addition, the use of regularization and cubic interpolation may lead to over-smoothing of the strain field. Second, the regularization used in our work penalizes strain magnitude, as opposed to strain gradients, and includes a spatially modulated regularization parameter. Since our regularization penalizes strain magnitudes, and we know a priori that the strain magnitude in a blood vessel will always be larger at the inner lumen, decaying with the inverse of the distance from the lumen squared, we can utilize this knowledge with a spatially varying regularization.
The algorithm evaluation performed here also differs from the evaluation of the algorithm in Liang et al. (2008) in several important ways. First, our algorithm was tested on both simulated B-mode and RF images, with varying signal-to-noise ratios, for a large range of strains. The simulated images created for our analysis mimicked both the image-formation process, which would result in more realistic image decorrelation with strain, and a vessel deformation with some variation in the material properties. This leads to strain fields that were inherently not smooth at material boundaries. In addition, several different initialization schemes were tested when using both B-mode and RF images. One initialization provided the true displacement field as a starting point to the algorithm to establish a measurement upper limit. Another provided a zero guess as an initialization. The third used a pre-registration, CC-based technique to estimate the radial (axial) displacements and the fourth used a mesh-refinement strategy designed to perform multiple registrations, varying the mesh resolution from low to high.
Figure 5a and e demonstrate that RF images had the potential to yield more accurate displacement and strain measurements than B-mode images, assuming that the issues of local minima could be overcome. When the RF images were used (Fig. 5e), very little difference was appreciated between the images with varying degrees of electronic noise. The plots of Figure 5b and f suggest that a zero initialization caused the algorithm to get trapped in local minima and resulted in decreased measurement accuracy. The differences between these plots show that the RF image measurements suffered from the problem of local trapping to a much larger degree than did the B-mode measurements. Indeed, the use of a zero initialization scheme would seem impractical for RF images.
Figure 5c and g demonstrate that the use of a cross-correlation-based measurement overcame the issue of local trapping, at least at smaller strains (<3%). The lower SNRe at larger strains suggested that the accuracy of the cross-correlation initialized measurement decreased beyond a 3% strain. In Figure 5d and h, the B-mode image measurements suggested that a mesh-refinement scheme yielded higher measurement accuracy than did zero initialization for all strain levels; however, mesh refinement yielded a lower accuracy than CC-based initialization. Thus, although a mesh-refinement initialization may indeed overcome some of the local minima issues for smaller strains, a cross-correlation scheme yields more accurate results. For the measurements using the RF images, the mesh refinement scheme is a slight improvement over a zero-initialization for very small strains; however, the cross-correlation-based initialization remains the best scheme using RF images for most strain levels considered.
Figure 6 demonstrates that the circumferential strain measurements had a higher accuracy at a given strain level. This was the case, at least in part, because the equation we used for the circumferential strain contained a radial displacement term, which had a higher accuracy than those in which a derivative is necessary. In practice, the accuracy of the circumferential strains may be largely dependent on the position of the catheter within the vessel and the particular tissue deformation. This is further shown in the strain elastograms of Figure 9.
Although the purpose of the regularization used in this work was to smooth the displacement measurement, this method can have the effect of over-penalizing the strains at the boundaries of the measurement domain. InFigure 10a, it is clear, particularly for the B-mode measurements, which the strains were relatively lower at the inner and outer radii of the images. This is further seen in Figure 10b, which shows the same lines of strain measurements. Although these boundary artifacts were always present, use of the RF images did reduce them. In addition, the application of the spatially varying regularization used here also mitigated the presence on this bias at the inner lumen by reducing the relative regularization strength.
Figure 11 demonstrates that at 2% strain, the results of the B-mode image measurements were very similar, showing only slight variations in the shapes. However, the RF image measurements showed that both the zero-initialization method and the mesh-refinement scheme completely failed to represent the expected strain profile. The cross-correlation initialization resulted in the most accurate of the image-initialization schemes, not including the true-displacement initialization. For the RF image measurements shown here, the results of true-displacement initialization and the CC- based initialization are approximately the same.
Figure 12 demonstrates the result of displacement measurements, both radially and angularly, performed on the phantom images for the average 2% strain case. The measurement of the angular displacement field shown here was slightly obscured by an imaging artifact known as the non-uniform rotational distortion, which is caused by the friction between the catheter sheath and the imaging wire as the wire rotates (Kimura et al. 1996). Although this distortion can be minimized via experimental precautions, it is difficult to remove in practice. When measuring image deformations, the consequence of such a distortion is that an angular displacement, constant radially at a given angle but different at each angle, will be added to measurements of the tissue displacements. Although this distortion might only be on the order of one to two A-lines (±3°), it can be significant relative to the tissue deformation, particularly at the low strain values that are found away from the inner lumen. This distortion was not modeled in our simulation, and it might affect the choice of regularization, since it alters the measured strain magnitude. In this work, these effects were neglected and the results were presented as is. However, in the future, we intend to investigate the possibility of using material-model constraints to distinguish and separate the tissue deformation and the distortion displacements as part of the registration algorithm.
It is clear from Figure 12 that the RF image measurements of the phantom had a much higher CNRe than those of the B-mode image measurements for the strains shown. This was further evident in the strain elastograms of Figure 13. Here the lower-magnitude radial strain and higher radial strain showed similar shapes, differing only in magnitude and apparent noisiness. The non-uniform distortion seemed to be apparent in the lower-magnitude circumferential strain, manifesting as angular fluctuations in the strain near the phantom’s inner lumen. This effect was less true in the higher-magnitude circumferential strain image.
The limited number of animals used in the in vivo pilot study does not allow us to make any conclusions regarding the clinical impact of the proposed technique. Nonetheless, Figure 14 demonstrates (1) the proposed method can produce useful strain elastograms under more realistic physiologic conditions; and (2) the strain elastograms were consistent with histologic observations.
We have developed and implemented an image-registration-based algorithm to measure the radial and circumferential motion induced within the vessel tissue using high-frequency intravascular ultrasound images. The results of computer simulations predicted that this strain estimator should provide reliable elastograms of vessels with strains ranging from 3%–5%–radial strain elastograms with SNRe as high as 4.5, and circumferential strain elastograms with SNRe as high as 7.5. The results of our initial in vivo study (phantom and an atherosclerotic rabbit) suggest that the non-rigid registration approach should provide useful radial and circumferential strain elastograms under realistic physiologic conditions and were sufficiently encouraging to justify a more extensive in vivo validation.
This work was funded by a National Heart and Lungs Research Grant R01 HL088523. The authors thank Dr Ebo de Muinck at Dartmouth College for his assistance in the rabbit study.
A series of reconstructions was performed for each image type (B-mode and RF), and each noise level (8, 12 and 17 dB), at each strain level (0.1%–10%) to determine the optimal choice of the regularization parameter for all images. Radial strain was calculated from the resulting displacement measurements and the radial strain SNRe was then calculated using eqn (7). This process was repeated for a series of varying magnitudes of the regularization parameter (α) to determine the optimal value of α (see example plot in Fig. 15a). The optimal α magnitude was then plotted for each image type (Fig. 15b and c). To obtain an objective choice for the regularization parameter used in this work, the plots were fitted with an inverse strain function, as this was the most accurate predictor of the regularization parameter. Although slight differences were apparent between the noise level and the measurement accuracy, due to its unpredictable nature, no noise dependence was modeled in the choice of α. It was found that the best-fit function for the regularization parameter of the B-mode image measurements was approximately α = 1 × 1018/ε avg; thus, it was this function that was used for the B-mode image registrations performed in the bulk of this article. The regularization parameter function of the RF image measurements was found to be approximately α = 5 × 1018/ε avg; thus, it was this function that was used for the RF registrations performed in the bulk of this article. These same equations were used in the case of the vessel phantom images, as we assumed that the simulation images closely resembled the phantom IVUS images. For the phantom images, εavg was estimated from the cross-correlation-based, initialization measurements.