Open Access
Add to BibSonomy
Abstract
Catheter-based intravascular optical coherence tomography (IVOCT) is a powerful imaging modality for visualization of atherosclerosis with high resolution. Quantitative characterization of various tissue types by attenuation coefficient (AC) extraction has been proven to be a potentially significant application of OCT attenuation imaging. However, existing methods for AC extraction from OCT suffer from the challenge of variability in complex tissue types in IVOCT pullback data such as healthy vessel wall, mixed atherosclerotic plaques, plaques with a single component and stent struts, etc. This challenge leads to the ineffectiveness in the tissue differentiation by AC representation based on single scattering model of OCT signal. In this paper, we propose a novel method based on multiple scattering model for parametric imaging of optical attenuation by AC retrieval from IVOCT images conventionally acquired during cardiac catheterization. The OCT signal characterized by the AC is physically modeled by Monte Carlo simulation. Then, the pixel-wise AC retrieval is achieved by iteratively minimizing an error function regarding the modeled and measured backscattered signal. This method provides a general scheme for AC extraction from IVOCT without the premise of complete attenuation of the incident light through the imaging depths. Results of computer-simulated and clinical images demonstrate that the method can avoid overestimation at the end of the depth profile in comparison with the approaches based on the depth-resolved (DR) model. The energy error depth and structural similarity are improved by about 30% and 10% respectively compared with DR. It provides a useful way to differentiate and characterize arterial tissue types in IVOCT images.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
1.1 Backgrounds
Intravascular optical coherence tomography (IVOCT) plays an important role in the diagnosis of coronary arterial atherosclerosis and the guidance of stent implantation. It is a catheter-based imaging modality capturing real-time pullback sequences of transverse and volumetric scanning images of arterial vessels with high resolution [1]. It provides detailed information regarding morphological structures and microstructures of vessel wall and atherosclerotic plaques based on backscattering of the incident near-infrared (NIR) laser in a spectral range of 1250-1350 nm [2]. Figure 1 shows the physical overview of a spectral-domain IVOCT system.
Interpretation of tomographic structural images allows for qualification of dissimilar tissue compositions. However, a comparison mechanism of physiological information in addition to morphological structures is required to obtain parameters with clinical diagnostic value that aid the differentiation of tissue components along with determination of plaque vulnerability. Quantification of tissue optical properties has been proven to be able to assist OCT image interpretation. The attenuation coefficient (AC) has been shown a significant indicator of optical properties describing the extinction of the detected OCT signal with the depth due to optical absorption and scattering [3]. The ACs extracted from IVOCT B-scans or three-dimensional (3-D) pullback data can be used to differentiate thin fibrous cap, calcification and lipid rich plaque and characterize the formation of neointima after stent implantation [3,4]. The combination of OCT attenuation with texture features and geometric structures segmented from structural images has been demonstrated to improve the accuracy of tissue characterization and vulnerable plaque identification [5,6].
1.2 Related work
Based on the fact that the backscattered light is exponentially attenuated with the detection depth in tissues, the AC is usually estimated from the measured backscattered signal by curve fitting (CF) through single A-line analysis [7–9] or 3-D analysis [10–12]. The estimation can be optimized by extending the fitting window in those regions with homogeneously distributed optical properties [13]. Moreover, CF has been combined with machine learning to improve the accuracy of tissue characterization [14,15]. The CF methods assume a constant attenuation along a part of the A-scan. Hence, they are most suitable for use on the tissue layers with well-defined layer structure, which subsequently reduces the automation of the methods. Besides, for thin layer structures, AC extraction by CF is unstable and sensitive to the noise in the measured signal. The small variations in the ACs cannot be accurately identified, limiting the OCT potential of high resolution.
In recent years, the depth-resolved (DR) technique has been an efficient tool for AC extraction from OCT with pixel-level resolution and high automation. Vermeer et al. [16] designed a DR method for the first time which assumes complete attenuation of all incident light in the range of the imaging depth, the absence of absorption and a constant scattering anisotropy. However, overestimation is observed in the areas where the backscattered light is not fully attenuated, such as eyelid area or hair fiber. Beside this, Vermeer’s method takes the maximum detection depth as the upper limit of the integral of the backscattered signal, leading to errors near the end of an A-line due to incorrect background subtraction. Furthermore, it assumes a constant confocal function, resulting in inaccurate estimates when the focal plane is within the sample or on the sample surface.
Amaral et al. [17] extended Vermeer’s DR model to a general model by quantifying the attenuation of the NIR light with the transmittance. Their results show that the ACs of various tissue types are estimated by inputting an appropriate transmittance, even for those thin tissue layers with the thickness of 100 µm. However, this method does not take into account the multiple scattered light, leading to wrong estimates from the weak OCT signals at the tissue boundaries.
Liu et al. [18] developed a DR-CF approach by combining the conventional DR model with CF to determine the cutoff depth. It preserves the CF’s advantage of high tolerance to noise and avoids wrong estimates caused by the weak backscattered signals at the tissue boundaries. However, it does not take into account the confocal effect and requires the focal plane to be located above the sample, reducing the signal-to-noise ratio (SNR) of the measurements and making some certain applications impractical.
Smith et al. [19] designed a depth-resolved confocal (DRC) method which allows for accurate estimation of the ACs when the focusing plane is located within the sample by taking into account both the confocal function and the sensitivity fall-off. But, overestimation is still observed in the presence of incomplete attenuation of the incident light in deep tissues. To address Vermeer’s assumption of full attenuation, Liu et al. [20] developed an optimized depth-resolved estimation (ODRE) method to implement AC extraction in any depth range regardless of whether the light is completely attenuated or not. ODRE is inapplicable to arterial vessels in IVOCT images because OCT signals attenuate severely in arterial tissues with large ACs such as lipid plaques, leading to an inaccurate threshold of the SNR.
Both DRC and ORDE require the prior knowledge of the confocal function parameters of the OCT setup employed (i.e., the focal plane position and apparent Rayleigh range), which makes some applications impractical. To address this issue, Dwork et al. [21] proposed a fully automated DRC method called autoConfocal method. It obtains the confocal function parameters automatically from the differences in the features between two images of the same sample captured at slightly different view angles. However, increasing errors are observed when the focal plane is far away from the sample because of the nearly constant confocal function. Besides, this method is unsuitable for use on IVOCT because it is infeasible to image the same vascular cross section from different views by using a single catheter pullback.
For clarity, Table 1 presents the overview on the major methods of OCT-based AC extraction.
Table 1. Overview of Major Methods of OCT-Based AC Extraction
1.3 Objective of The Work
In this paper, a novel model-based method is proposed for parametric imaging of attenuation by extracting ACs in arterial cross sections from transverse IVOCT scanning images. The given data is the conventionally acquired IVOCT gray scale images satisfying the DICOM digitalization standard as well as the configuration parameters provided by the imaging system including the incident light power, central wavelength and full width at half maximum (FWHM). Pixel-wise AC estimates are achieved by iteratively minimizing an error function regarding the measured and modeled backscattered signal which is the output of a forward model based on Monte Carlo (MC) simulation. This method provides a general scheme for AC estimates without the premise of complete attenuation of the incident light in tissues, as opposed to the DR-based approaches. Experiments of numerical simulations and clinical images are conducted to verify its validity. The influence of initial plans of iteration is discussed. The experiment of comparing with the DRC method is also conducted.
2. Method
Our method aims to find the AC distribution in a vascular cross section from an IVOCT gray scale image by solving the forward problem iteratively, updating the ACs at each iteration, until the output of the solver matches the measured OCT signal in a framework of least-squares minimization. That is, the squared error between a measured OCT signal and the output of a forward model is minimized by adjusting the AC. This section gives the detailed steps.
2.1 Forward modeling of IVOCT signal
AC extraction from OCT requires a quantitative model to describe the OCT signal. Two general models have been proposed, i.e., single-scattering (SS) model [3,4] and multiple-scattering (MS) model [22,23]. Model selection depends on the tissue type of interest. The SS model assumes that the incident light propagates forward in the medium and is attenuated by absorption and scattering satisfying the Lambert-Beer's law until it is reflected back to the source by backscattering. It is based on the premise of single backscattering of photons, thus it is suitable for the AC extraction of a weakly scattering medium or the superficial thin layer of the highly scattering tissue. The MS model considers multiple backscattering events and large detection depth. Till now, there have been three approaches taking MS into account in OCT. In the order of increasing complexity and number of parameters, they are stochastic approach based on MC simulation [24–26], analytic approach based on extended Huygens–Fresnel model [27] and full electromagnetic wave modeling based on Maxwell’s equations [28]. In comparison with the MS model, the SS model has a key advantage of simplicity [3], hence it is the most commonly used in AC extraction such as CF and DR-based methods. However, clinical observations have demonstrated that multiple scattering can be expected in highly forward scattering tissues such as blood and atherosclerotic plaques [3,29].
In this work, we employ MC simulation to construct the forward model of OCT signal. The model includes three parts: Gaussian beam generation, photon propagation in the medium and backscattered signal simulation. For simplicity, we ignore the interface reflection between the optical elements as well as the phase change during light propagation [30]. As illustrated in Fig. 2, in an imaging plane, the catheter remains in the image center, although it does not stay fixed in the lumen center along the lumen length. The lumen, vessel wall intima/media, plaques and adventitia are surrounding the catheter center from inside to outside. A rectangular coordinate system XOY and a polar coordinate system θ-l are established in the imaging plane where the origin O is located at the image center.
2.1.1 Gaussian beam generation
The position and direction of the photon emission in the Gaussian beam is simulated by MC method. Then, the propagation trajectory of the photon in the focusing system is traced generating the distribution function on the back surface of the lens. More details on the principle and code of MC simulation of Gaussian beam generation can be found in Ref. [31].
2.1.2 Photon propagation in the tissues
On a microscopic scale, each vessel wall layer presents uniform optical property, while their mirostructures are regarded as heterogeneous turbid medium. MC simulation of photon propagation depicts the propagation trajectory of photons in a weakly absorbing and strongly scattering media by regarding photons as typical particles and ignoring the properties of polarization and phase of light. As illustrated in Fig. 3, photons are launched perpendicular to the tissue surface. After multiple scattering and absorption, partial photons escape from the medium and are finally received by the detector.
Fig. 3. Schematic diagram of scattering events of photons in tissue. θ and φ denote the deflection angle and azimuthal angle, respectively.
As illustrated in Fig. 4, the detailed steps are as follows.
- (1) Photon initialization
A single photon of weight W is launched vertically in the medium surface at an initial position (0,0,0). Its direction cosine and weight are initialized as (0,0,1) and 1, respectively.
- (2) Sampling of photon step-size
Assume the free path of a photon between two successive elastic scattering events is determined by the Poisson probability density function with respect to the AC and the photon free path [30]. The free path-length, i.e., the step-size of photon propagation is given by [31]
where ξ is a computer-generated random number uniformly distributed in the interval [0,1] and µt is the AC describing the decay of the incident light due to the tissue optical properties. The pathlength Li of the i-th photon is obtained by the product of its total steps and the refractive index of the medium. - (3) Photon energy absorption
The partial energy is absorbed by the medium once a photon moves by the distance s to the scattering site, resulting in the reduction in the photon weight W, i.e., ΔW=µaW/µt, where µa is the optical absorption coefficient [33].
- (4) Photon scattering
After absorption the photon undergoes the scattering event causing it is redirected. The azimuthal angle φ is uniformly distributed from 0 to 2π and is obtained by direct sampling, i.e., φ=2πξ. The cosine of the deflection angle θ∈[0,π] follows the Henyey-Greenstein probability distribution function [32] as follow,
Direct sampling gives [33]
$$\cos \theta = \left\{ {\begin{array}{ll} \frac{1}{{2g}}\left[ {1 + {g^2} - {{\left( {\frac{{1 - {g^2}}}{{1 - g + 2g\xi }}} \right)}^2}} \right], &{g > 0} \\ {2\xi - 1,}&{g = 0}\end{array}} \right., $$where g is the scattering anisotropy factor (g=1 for multiple scattering in multi-layered vessel wall tissues implying that all the light is scattered in the forward direction [22,34]).The penetration depth of the photon in the medium is given by
where n denotes the total times of a single photon scattering in the medium, sj denotes the step size of the j-th photon and µz,j is the direction cosine of the j-th scattering, - (5) Ending of photon propagation
Tracking of the photon ends when its weight is less than the threshold (10−6) or when it escapes the geometry boundaries of the medium.
2.1.3 OCT signal simulation
The MC simulation of photon propagation outputs the weights and pathlengths of the photons collected from the medium surface. Then, the OCT signal backscattered by tissues in terms of the radial depth ri is determined by [30]
Here, I0 is a constant determined by the instrumental properties of the OCT system, Nph is the total number of photons launched (typically ∼107), Wi is the weight of i-th detected photon with optical pathlength Li, 2ri is the optical pathlength in the reference arm, λc is the central wavelength of the incident radiation, and Δλ is the FWHM of the laser spectrum. Eq. (6) accounts for the speckle structure caused by the overlap of interference patterns. Speckles are eliminated in simulations by using the following expression [30] instead of Eq. (6),
2.2 AC estimation by iterative optimization
The standard output of a commercially available IVOCT system is sequential pull-back gray scale images with 8-bit gray levels overlaid with golden color maps as shown in Fig. 2. The intensity of the measured backscattered signal at each radial depth in a vascular cross section is directly proportional to the gray value of the corresponding image pixel,
where Im(ri) is the measured OCT signal and g(ri) is the 8-bit gray value (ranging from 0 to 255) of the pixel at ri in a standard IVOCT B-scan image.An error function regarding the measured and modeled OCT signal is defined as follow,
where f(·) is the error function to be minimized, µt(ri) is the AC and IOCT(ri) is the normalized theoretical OCT signal (ranging from 0 to 1). For the convenience of expression, we denote µt(ri) as µ. A nonlinear least squares (NLS) problem is defined as follow, where ${\hat{\mu }_\textrm{t}}({r_i})$ is the estimated AC at ri and F(µ) is the objective function, F(µ) is linearly approximated by the first-order Taylor expansion,As shown in Fig. 5, the iterative steps for solving Eq. (14) by Levenberg-Marquadt (L-M) optimization are detailed as follows.
Step 1. Initialization.
Initialize the involved parameters including µ0, γ0, and k=0, where µ0 is the initial AC at ri, γ0 is the initial damping factor and k is the number of iteration. µ0 can be roughly chosen by referring to the prior histological knowledge about the tissue types. γ0 is usually initialized as 0.01 [35]. In Discussion, we will discuss the dependence of the AC estimates on the initial plans of iteration.
Step 2. Calculate increment of iteration.
Solve the minimum of m(Δµ) by letting $m^{\prime}(\Delta \mu)=0$, thus
Consider
Calculate the similarity index that measures the similarity between the linear approximation and the real objective function as follow,
Step 4. Adjust the damping factor depending on ρ and determine µk+1.
- (i) ρ≤0 corresponds to a poor linear approximation, thus determine µk+1 by
Then, reduce the damping factor, i.e., γk+1=0.5γk, and proceed to Step 5.
- (ii) ρ>0 corresponds to a good linear approximation, thus determine µk+1 by
Then, adjust the damping factor as follows,
$${\gamma _{k + 1}} = \left\{ {\begin{array}{c} {\begin{array}{cc} {0.5{\gamma_k},}&{0 < \rho \le 0.25} \end{array}}\\ {\begin{array}{cc} {{\gamma_k},}&{0.25 < \rho \le 0.75} \end{array}}\\ {\begin{array}{cc} {2{\gamma_k},\begin{array}{cc} {}&{} \end{array}}&{\rho > 0.75} \end{array}} \end{array}} \right.$$and proceed to Step 5.
Step 5. Ending of the iteration.
If
proceed to Step 6; otherwise, the iteration ends and output µk+1.Step 6. Update k by k←k+1 and turn to Step 2.
The AC at each pixel position in an image is estimated by above iterative steps. After normalization and pseudocolor coding, the AC matrix is transformed into a parametric image.
2.3 Performance evaluation
We evaluated the proposed method experimentally using simulated images and clinical images. The experiments of numerical simulation aim to quantitatively evaluate the accuracy of AC extraction, discuss the influence of the initial plans of iteration and compare with other methods. The experiments on clinical pullback image data aim to demonstrate the validity of our method in retrieving ACs of arterial tissues and characterizing tissue types.
We implemented the method by programming with MATLAB (2016a, The MathWorks, Inc., Natick, MA) on a personal computer configured with a 2.6 GHz Intel Core I7-4510U as CPU, 8GB RAM and Windows 10 64 bits as the operating system. The input data involved in the method are the total number of photons launched, the initial AC, and the central wavelength and FWHM of the incident laser which are provided by the imaging probe.
2.3.1 Clinical image preparation
The clinical images were acquired using a commercially available spectral-domain OCT system (C7XR OCT system, LightLab Imaging/St. Jude Medical Inc, St Paul, MN). A 2.7F C7 Dragonfly OCT catheter was inserted into the target vascular lumen and pushed to the distal end through the culprit lesion. It was then drawn back automatically at 20 mm/s, while blood was displaced by an injection of contrast media at rates of 2 to 4 ml/s. The transversal images at an A-line rate of 120 kHz with the axial spatial resolution of 15 µm and transversal resolution of 25 µm were digitally stored for off-line analysis. Table 2 lists the operating and characteristic parameters of the setup.
Table 2. Parameters of Imaging System in Forward Simulation
2.3.2 Simulated image preparation
In simulation experiments, we constructed the cross-sectional geometry of computer-simulated arterial phantoms by referring to the histological structures of healthy coronary arterial wall and different types of atherosclerotic plaques. We employed Microsoft Visio (v.2016, Microsoft Corporation, Redmond, State of Washington) to draw the vascular cross sections. Figure 6(a) provides four examples, in which the structures of the arterial wall and plaques are referred to the histology in Fig. 6(b) [18]. The radius of the vascular cross sections is 2.0 mm.
Fig. 6. Computer-simulated coronary arterial phantoms. (a) Cross-sectional geometry of the phantoms which are numbered 1#, 2#, 3# and 4# from left to right; (b) Histological cross sections of coronary arterial vessels and atherosclerotic plaques (Hematoxylin and eosin stain) (Ref. [18]).
We simulated IVOCT images basing on the principle of OCT signal formation as described in Sec. 2.1. The instrumental parameters of the OCT system were in accordance with the characteristics of the clinical system described above. We added the Guassian random noise with the SNR of 55 dB to the simulated OCT signal to be more close to the real situation.
2.3.3 Figures of Merit
We adopted the error percentage, energy error depth (EED) [19] and structural similarity (SSIM) [36] to quantitatively evaluate AC estimates. The error percentage is defined by
3. Results
3.1 Results of simulated images
Figure 7 provides the results of the computer-simulated vessel phantoms. The simulated IVOCT images are 300×300 pixels (8 mm×8 mm) with 8-bit gray levels. For simulation of an A-scan, 1×107 photons were employed. The A-scans from each 360° rotation were arranged in a B-scan. In the iteration µ0=0.08 mm−1 and γ0=0.01. These results reveal that the AC maps display the multi-layered structures of the vessel wall and plaques with significantly improved contrast in comparison with the structural images. Table 3 lists the true and estimated AC values as well as the error percentages lower than 10%. These convincing results indicate that our method is able to estimate the ACs from IVOCT images with a high accuracy.
Fig. 7. Results of simulated demonstrations. (a) Phantom 1#; (b) Phantom 2# ; (c) Phantom 3#; (d) Phantom 4#.
Table 3. Estimated and Preset AC Values and Error Percentages for Simulated Arterial Phantoms
3.2 Results of clinical images
Figure 8 presents the AC maps retrieved from twelve typical frames of clinical IVOCT scanning images (300×300 pixels) containing healthy vessel wall, lesioned tissues and stent struts. The parameter settings were the same as the simulation experiments. The regions of interest (ROIs) including plaque areas and stent struts were manually labeled by experienced physicians. The images of AC mapping provide information regarding the optical properties of different arterial tissue types, which can be applied in the characterization and differentiation of vulnerable plaques as well as the guidance and assessment of minimally invasive coronary interventions. These results demonstrate the feasibility of using our method to quantitatively characterize and differentiate at least seven different arterial tissue types, i.e., healthy vessel wall (trilaminar structure of intima (signal-rich or high backscattering), media (signal-poor or low backscattering) and adventitia (heterogeneous and high backscattering)), lipid-rich plaque (low-intensity and irregular signal with an indistinct and irregular border, strong attenuation), calcification (signal-poor or heterogeneous region with a sharp defined border), fibrous plaque (signal-rich or high backscattering, homogeneous and weak attenuation), fibrous cap (a signal-rich tissue layer, overlying a signal-poor region), plaque rupture (fibrous-cap disruption along with cavity formation), thrombus (an irregular mass attached to the luminal surface, including red thrombus (red blood cell-rich, high backscattering and high attenuation) and white thrombus (platelet-rich, low backscattering, homogeneous and low attenuation)) [37–40]. The quantitative tissue characterization by AC extraction is helpful to the qualitative tissue classification that relies on the interpretation of image textures and structural features. Besides, the stents after implantation can be assessed from the AC maps, including neointima hyperplasia, tissue prolapse (tissue protrusion into the lumen between stent struts without apparent surface disruption), and apposition of stent struts (malapposition: the axial distance between the strut’s surface to the luminal surface is greater than the strut thickness; well apposed: the distance is less than the strut thickness [37]).
Table 4 presents the ranges of the histological Refs. [7,39] and estimated AC values of different tissue types. In the table, the ranges of the ACs were measured from the manually selected regions in the AC maps shown in Fig. 8. These results indicate that the AC estimates obtained by our method are basically consistent with the histological references.
3.3 Comparison results
In this section, the results of comparing our method with Smith’s DRC method [19] are given. In our method, the initial plan of iteration is µ0=0.08 mm−1 and γ0=0.01. DRC utilizes intensity weighted horizontal total variation (iwhTV) regularization to optimize the estimated ACs, where the regularization parameter is η=3. Figures 9–11 present the results of the simulated images and Fig. 12 provides the results of the clinical images.
Fig. 9. AC values along one radius of the vessel phantoms estimated by using our method and DRC. (a) Phantom 1# (90°); (b) Phantom 2# (105°); (c) Phantom 3# (190°); (d) Phantom 4# (75°).
Fig. 10. AC maps of the true value and estimates by our method and DRC. (a) Phantom 1#; (b) Phantom 2#; (c) Phantom 3#; (d) Phantom 4#.
Fig. 11. Comparison of figures of merit and time cost of AC estimates by our method and DRC. (a) EED; (b) SSIM; (c) Time cost.
In our simulated phantom preparation, not all the incident light is assumed to be completely attenuated at depths ranging from 0 to 2.0 mm. As indicated in Fig. 9, the overestimation at the end of the depth profile (i.e., at depths ranging from 1.7 mm to 2.0 mm) is observed in the estimates by DRC because of its premise of full attenuation. This conclusion can also be drawn from the visualization of the AC maps shown in Fig. 10. In contrast, the estimates at the end of the depth profile by using our method are closer to the true values than the DRC estimates. Figure 11 suggests that both the EED and SSIM of the estimates by our method are improved by up to 30% and 10% in comparison with those by DRC. The effective range of the AC estimation with our method is the entire radial depth, i.e., 2.0 mm, while that of DRC is about 1.7 mm. Our method has a relatively higher time cost than DRC because of its numerical iterative nature in contrast to DRC’s analytical nature. Fortunately, this computational burden can be greatly reduced by using parallel codes and a high-performance parallel computing platform that exploits graphics processing units (GPUs).
Moreover, DRC has been demonstrated to be efficient in extracting ACs in horizontal layered tissues such as bladder, skin and retina by utilizing the correlation of A-scan data to smooth tissue layers [19]. However, coronary arteries with atherosclerotic lesions usually present complicated anatomical structures and have various tissue compositions as well. Hence, a vascular cross section in a transverse structural image may involve multiple tissue types, showing complex structures instead of horizontal layered structure. Results in Fig. 12 are consistent with this conclusion, where the transverse ambiguity is observed in the DRC estimates. Our method outperforms DRC in the contrast of different tissues in the AC maps, especially at the end of the depth profile.
4. Discussion
4.1 Dependence on initial plans of iteration
We discussed the influence of the initial plans of iteration on the AC extraction based on the simulation results. First, we compared the AC estimates obtained by changing the initial AC µ0, while the damping factor remained unchanged, i.e., γ0=0.01. Figure 13 suggests small differences in the AC estimates when µ0=0.02 mm−1, 0.04 mm−1, 0.06 mm−1 and 0.08 mm−1. Figure 14 indicates that the EED and SSIM of the AC estimates when µ0=0.08 mm−1 are slightly higher than those of other cases with the least times of iteration. Second, we set the damping factor γ0=0.01, 0.02, 0.03 and 0.04 respectively, while µ0=0.08 mm−1. The results in Fig. 15 reveal that the lowest iteration times and time cost are achieved when γ0=0.01. In the L-M optimization, the damping factor is adjusted dynamically during iteration, so it has major influence on the iteration times and time cost, but minor influence on the iteration result.
Fig. 13. Estimated values of AC in the case of different µ0. (a) Phantom 1# (90°); (b) Phantom 2# (75°); (c) Phantom 3# (105°); (d) Phantom 4# (190°).
Fig. 14. Figures of merit of AC maps obtained in the case of different µ0. (a) EED; (b) SSIM; (c) Iteration times.
Fig. 15. Figures of merit of AC maps obtained in the case of different γ0. (a) Iteration times; (b) Time cost; (c) EED; (b) SSIM.
4.2 Discussion on model-based scheme
The AC estimation from OCT measurements is a distributed parameter estimation problem. It is essentially the inversion of a light transport operator. We designed a model-based minimization scheme for the recovery of AC as a function of position from IVOCT images. Two main aspects are involved in this scheme: forward model and optimization framework. We will discuss both issues in this section.
4.2.1 Model selection
The selection of the OCT signal model is critical in extracting AC from OCT data and it should be based on the chemical compositions and biological characteristics of tissue regions of interest. Considering the simplicity of the SS model, many efforts have been made to improve the accuracy and adaptability of the SS-based schemes such as CF and DR methods. For IVOCT, clinical applications have demonstrated that the effects of multiple scattering should be considered in coronary arterial vessel walls with multi-layered structures and atherosclerotic plaques [3,29]. Hence, AC measurements with higher accuracy and stability will be obtained by taking into account multiple scattering of photons in comparison with single scattering although it may be computationally intensive.
When light propagating in a turbid medium, the attenuation is mainly caused by absorption and scattering which are parameterized by the absorption coefficient µa and scattering coefficient µs of tissues, respectively. The sum of µa and µs, i.e., the AC µt, represents tissue optical properties. In OCT operating in a NIR spectrum window, the attenuation is dominated by the scattering while the absorption by those optical absorbers is very low to achieve high enough imaging depths. In the CF- and DR-based methods, the loss of the backscattered signal with the penetration depth caused by optical absorption and scattering is parameterized by a decay coefficient µOCT and then is modeled as a single exponential decay governed by the Lambert-Beer law. This simple parameterization greatly alleviates the complexity, but sacrifices the direct relationship with the optical properties of tissues (e.g., anisotropy factor g, µa and µs) to achieve robust and unique measurement [4]. Moreover, the estimated µOCT may be theoretically lower than µt in the presence of multiple scattering events because more light will be detected than expected by the SS model due to multiple scattering.
In the model-based minimization scheme of a parameter estimation problem, the forward signal model is required to be sufficiently accurate to capture the essential characteristics of the light field. Besides, it should be fast enough computationally to make its use in iterative inversion methods, where it may need to be evaluated many times [41]. In biomedical optics, the MC method represents the gold standard for modeling light transport in biological tissues, both due to its accuracy and its flexibility in modeling realistic, heterogeneous, multilayered turbid media [42]. Since both the effects of multiple scattering and dependent scattering play an important role in the quantification of optical tissue properties [43], we will take into account the dependent scattering by utilizing the well-established MC simulation software in our future work as has been done in [42]. Moreover, because of the inherently parallel nature of MC simulations, the time cost could be reduced significantly by implementing the simulator on GPUs using the compute unified device architecture (CUDA) programming language [44] as the availability of large computing clusters and clusters of GPUs is increasing quickly [31].
4.2.2 Optimization framework
We used the least-squares minimization, a common and well-developed framework to solve the inverse problem, where we solved the NLS problem with a L-M algorithm. The initial plans of iteration are required to be properly selected according to the prior histological knowledge of the tissue types and components. It may lower the automation of the method. In addition, as an iterative scheme, the computational burden that might limit the real-time applications should be considered. An efficient optimization tool such as Bregman [45] and fast iterative shrinkage thresholding algorithm (FISTA) [46], parallel codes and a high-performance parallel computing platform that exploits GPUs being used to reduce the computational burden may improve the efficiency of our method.
In recent years, deep learning (DL) has shown great potential in the field of medical imaging with the rapid growth of high-performance processors and big data technology. DL, especially convolutional neural network (CNN), has been proven an efficient tool in medical image analysis, processing and reconstruction [47]. Many efforts have also been made in DL-based approaches of quantitative parametric imaging, such as quantitative photoacoustic imaging which aims to solve an ill-posed inverse problem of recovering optical absorption coefficient and scattering coefficient from acoustic measurements [48]. Big data is crucial for DL. IVOCT has been routinely applied in minimally invasive cardiac interventions, so that a large amount of clinically acquired image data is available, which can be used to build training data sets. Hence, in the future work, we plan to develop an automated data-driven DL method for quantitative IVOCT by constructing a recovery function in the form of CNN to map an input structural image to an output parametric image of attenuation.
5. Conclusion
We proposed a novel method based on error-minimization as an efficient tool of retrieving pixel-wise AC distribution in vascular cross sections from IVOCT images to implement parametric imaging of optical attenuation. We solved the NLS problem regarding the measured and modeled backscattered signal being obtained by forward MC modeling. Results of simulated data demonstrated the accuracy of AC extraction with the error percentage less than 10%. Results of clinical images demonstrated the consistency of AC estimates of different tissue types with the histological references. The parametric imaging of AC enables the differentiation and characterization of at least seven tissue types, i.e., healthy vessel wall, lipid-rich plaque, calcified plaque, fibrous plaque, fibrous cap, plaque rupture and red/white thrombus, as well as stent apposition (well apposition, malapposition, intima hyperplasia and plaque prolapse). Comparison results validate the improvement of up to 30% in the EED and 10% in the SSIM of the parametric imaging of attenuation in comparison with the DRC approach.
One issue that must be pointed out that we evaluated the accuracy of our method quantitatively with computer-simulated vessel phantoms. However, it is evident that an AC estimation scheme should be further evaluated through validation measurements of tissue phantoms with well specified attenuation coefficients as has been done in [49]. It will be another work for us in the further study.
Funding
National Natural Science Foundation of China (62071181).
Disclosures
None.
References
1. S. Vallabhajosyula, S. C. El Hajj, M. R. Bell, A. Prasad, A. Lerman, C. S. Rihal, D. R. Holmes, and G. W. Barsness, “Intravascular ultrasound, optical coherence tomography, and fractional flow reserve use in acute myocardial infarction,” Catheter. Cardio. Inte. 96(4), E59–E66 (2020). [CrossRef]
2. M. Y. Henein, S. Vancheri, G. Bajraktari, and F. Vancheri, “Coronary atherosclerosis imaging,” Diagnostics 10(2), 65 (2020). [CrossRef]
3. S. Chang and A. K. Ellerbee, “Review of methods and applications of attenuation coefficient measurements with optical coherence tomography,” J. Biomed. Opt. 24(09), 1 (2019). [CrossRef]
4. P. Gong, M. Almasian, G. van Soest, and D. M. de Bruin, “Parametric imaging of attenuation by optical coherence tomography: review of models, methods, and clinical translation,” J. Biomed. Opt. 25(04), 1 (2020). [CrossRef]
5. F. J. van der Meer, D. J. Faber, D. M. B. Sassoon, M. C. Aalders, G. Pasterkamp, and T. G. van Leeuwen, “Localized measurement of optical attenuation coefficients of atherosclerotic plaque constituents by quantitative optical coherence tomography,” IEEE Trans. Med. Imaging 24(10), 1369–1376 (2005). [CrossRef]
6. T. Imanaka, H. Hao, K. Fujii, M. Shibuya, M. Fukunaga, K. Miki, H. Tamaru, M. Masutani, M. Ohyanagi, and T. Masuyama, “Analysis of atherosclerosis plaques by measuring attenuation coefficients in optical coherence tomography: thin-cap fibroatheroma or foam cells accumulation without necrotic core?” Eur. Heart. J. 34(suppl 1), P5482 (2013). [CrossRef]
7. G. van Soest, T. Goderie, E. Regar, S. Koljenović, G. L. J. H. van Leenders, N. Gonzalo, S. van Noorden, T. Okamura, B. E. Bouma, G. J. Tearney, J. W. Oosterhuis, P. W. Serruys, and A. F. W. van der Steen, “Atherosclerotic tissue characterization in vivo by optical coherence tomography attenuation imaging,” J. Biomed. Opt. 15(1), 011105 (2010). [CrossRef]
8. A. Boi, A. D. Jamthikar, L. Saba, D. Gupta, A. Sharma, B. Loi, J. R. Laird, N. N. Khanna, and J. S. Suri, “A survey on coronary atherosclerotic plaque tissue characterization in intravascular optical coherence tomography,” Curr. Atheroscler. Rep. 20(7), 33 (2018). [CrossRef]
9. C. Xu, J. M. Schmitt, S. G. Carlier, and R. Virmani, “Characterization of atherosclerosis plaques by measuring both backscattering and attenuation coefficients in optical coherence tomography,” J. Biomed. Opt. 13(3), 034003 (2008). [CrossRef]
10. R. Shalev, M. Gargesha, D. Prabhu, K. Tanaka, A.M. Rollins, M. Costa, H.G. Bezerra, G. Lamouche, and D.L. Wilson, “Validation of parameter estimation methods for determining optical properties of atherosclerotic tissues in intravascular OCT,” in Proceedings of SPIE International Conference on Medical Imaging 2014: Image Perception, Observer Performance, and Technology Assessment, 15 Feb. 9037, 90371D (2014).
11. M. Gargesha, R. Shalev, D. Prabhu, K. Tanaka, A. M. Rollins, M. Costa, H. G. Bezerra, and D. L. Wilson, “Parameter estimation of atherosclerotic tissue optical properties from three-dimensional intravascular optical coherence tomography,” J. Med. Imaging 2(1), 016001 (2015). [CrossRef]
12. R. Shalev, M. Gargesha, D. Prabhu, K. Tanaka, A. M. Rollins, G. Lamouche, C. E. Bisaillon, H. G. Bezerra, S. Ray, and D. L. Wilson, “Processing to determine optical parameters of atherosclerotic disease from phantom and clinical intravascular optical coherence tomography three-dimensional pullbacks,” J. Med. Imaging 3(2), 024501 (2016). [CrossRef]
13. M. Gnanadesigan, A. S. Hussain, S. White, S. Scoltock, A. Baumbach, A. F. W. van der Steen, E. Regar, T. W. Johnson, and G. van Soest, “Optical coherence tomography attenuation imaging for lipid core detection: an ex-vivo validation study,” Int. J. Cardiovasc. Imaging 33(1), 5–11 (2017). [CrossRef]
14. G. J. Ughi, T. Adriaenssens, P. Sinnaeve, W. Desmet, and J. D’hooge, “Automated tissue characterization of in vivo atherosclerotic plaques by intravascular optical coherence tomography image,” Biomed. Opt. Express 4(7), 1014–1030 (2013). [CrossRef]
15. Y. Gan, D. Tsay, S. B. Amir, C. C. Marboe, and C. P. Hendon, “Automated classification of optical coherence tomography images of human atrial tissue,” J. Biomed. Opt. 21(10), 101407 (2016). [CrossRef]
16. K. A. Vermeer, J. Mo, J. J. A. Weda, H. G. Lemij, and J. F. de Boer, “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express 5(1), 322–337 (2014). [CrossRef]
17. M. M. Amaral, D. M. Zezell, A. F. G. Monte, A. C. B. de Cara, J. C. R. Araújo, A. Antunes, and A. Z. Freitas, “General model for depth-resolved estimation of the optical attenuation coefficients in optical coherence tomography,” J. Biophotonics 12(10), e201800402 (2019). [CrossRef]
18. S. Liu, Y. Sotomi, J. Eggermont, G. Nakazawa, S. Torii, T. Ijichi, Y. Onuma, P. W. Serruys, B. P. F. Lelieveldt, and J. Dijkstra, “Tissue characterization with depth-resolved attenuation coefficient and backscatter term in intravascular optical coherence tomography images,” J. Biomed. Opt. 22(9), 1–16 (2017). [CrossRef]
19. G. T. Smith, N. Dwork, D. O’Connor, U. Sikora, K. L. Lurie, J. M. Pauly, and A. K. Ellerbee, “Automated, depth-resolved estimation of the attenuation coefficient from optical coherence tomography data,” IEEE Trans. Med. Imaging 34(12), 2592–2602 (2015). [CrossRef]
20. J. Liu, N. Ding, Y. Yu, X. Yuan, S. Luo, J. Luan, Y. Zhao, Y. Wang, and Z. Ma, “Optimized depth-resolved estimation to measure optical attenuation coefficients from optical coherence tomography and its application in cerebral damage determination,” J. Biomed. Opt. 24(3), 1–11 (2019). [CrossRef]
21. N. Dwork, G. T. Smith, T. Leng, J. M. Pauly, and A. K. Bowden, “Automatically determining the confocal parameters from OCT B-Scans for quantification of the attenuation coefficients,” IEEE Trans. Med. Imaging 38(1), 261–268 (2019). [CrossRef]
22. M. Almasian, N. Bosschaart, T. G. van Leeuwen, and D. J. Faber, “Validation of quantitative attenuation and backscattering coefficient measurements by optical coherence tomography in the concentration-dependent and multiple scattering regime,” J. Biomed. Opt. 20(12), 121314 (2015). [CrossRef]
23. B. Karamata, M. Laubscher, M. Leutenegger, S. Bourquin, and T. Lasser, “Multiple scattering in optical coherence tomography. I. Investigation and modeling,” J. Opt. Soc. Am. A 22(7), 1369–1379 (2005). [CrossRef]
24. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef]
25. J. Kalkman, A. V. Bykov, D. J. Faber, and T. G. van Leeuwen, “Multiple and dependent scattering effects in Doppler optical coherence tomography,” Opt. Express 18(4), 3883–3892 (2010). [CrossRef]
26. M. Y. Kirillin, A. V. Priezzhev, and R. Myllylä, “Role of multiple scattering in formation of OCT skin images,” Quantum Electron. 38(5), 486–490 (2008). [CrossRef]
27. L. Thrane, H. T. Yura, and P. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens-Fresnel principle,” J. Opt. Soc. Am. A 17(3), 484–490 (2000). [CrossRef]
28. P. R. T. Munro, “Three-dimensional full wave model of image formation in optical coherence tomography,” Opt. Express 24(23), 27016–27031 (2016). [CrossRef]
29. D. Levitz, L. Thrane, M. H. Frosz, and P. E. Andersen, “Determination of optical scattering properties of highly-scattering media in optical coherence tomography images,” Opt. Express 12(2), 249–259 (2004). [CrossRef]
30. M. Kirillin, I. Meglinski, V. Kuzmin, E. Sergeeva, and R. Myllylä, “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express 18(21), 21714–21724 (2010). [CrossRef]
31. S. Malektaji, I. T. Lima Jr, and S. S. Sherif, “Monte Carlo simulation of optical coherence tomography for turbid media with arbitrary spatial distributions,” J. Biomed. Opt. 19(4), 046001 (2014). [CrossRef]
32. V. V. Tuchin, Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis (Third edition) (SPIE Press, 2015.
33. Y. Zhang, “Research on Monte Carlo simulation of optical coherence tomography system,” Master Thesis. Qingdao University, China (2009).
34. V. Periyasamy and M. Pramanik, “Importance sampling based Monte Carlo simulation of time domain optical coherence tomography with embedded objects,” Appl. Opt. 55(11), 2921–2929 (2016). [CrossRef]
35. K. Madsen, H.B. Nielsen, and O. Tingleff, “Methods for non-linear least squares problems. second ed.” Department of Mathematical Modelling, Technical University of Denmark, Lyngby, Denmark. 2004.
36. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13(4), 600–612 (2004). [CrossRef]
37. G. J. Tearney, E. Regar, and T. Adriaenssens, “Consensus standards for acquisition, measurement, and reporting of intravascular optical coherence tomography studies: a report from the international working group for intravascular optical coherence tomography standardization and validation,” J. Am. Coll. Cardiol. 59(12), 1058–1072 (2012). [CrossRef]
38. H. C. Lowe, J. Narula, J. G. Fujimoto, and I. K. Jang, “Intracoronary optical diagnostics-current status, limitations, and potential,” JACC-Cardiovasc. Inte. 4(12), 1257–1270 (2011). [CrossRef]
39. Y. Zhou, “Research on catheter-based Fourier domain OCT system and imaging methods,” Master Thesis. Tianjin University, China (2017).
40. Cardiovascular Imaging group, Interventional Cardiology group, Cardiovascular Society of Chinese Medical Association, “Advise from Chinese experts on application of optical coherence tomography in interventional diagnosis and treatment of coronary heart disease,” Chin. J. Cardiol. 45(1), 5–12 (2017).
41. B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. 17(6), 061202 (2012). [CrossRef]
42. A. Varkentin, M. Otte, M. Meinhardt-Wollweber, M. Rahlves, M. Mazurenka, U. Morgner, and B. Roth, “Simple model to simulate OCT-depth signal in weakly and strongly scattering homogeneous media,” J. Opt. 18(12), 125302 (2016). [CrossRef]
43. V. D. Nguyen, D. J. Faber, E. von der Pol, T. G. van Leeuwen, and J. Kalkman, “Dependent and multiple scattering in transmission and backscattering optical coherence tomography,” Opt. Express 21(24), 29145–29156 (2013). [CrossRef]
44. A. Kalra and I. Lima Jr., S.S. Sherif, “Almost instantaneous Monte Carlo calculation of optical coherence tomography signal using graphic processing unit,” Proceedings of 2013 IEEE Photonics Conference (IPC), 8-12 Sept., 2013, Bellevue, Washington, paper MA2.3 (2013).
45. H. Gao, H. Zhao, and S. Osher, “Bregman methods in quantitative photoacoustic tomography,” CAM Report 30(6), 3043–3054 (2010).
46. B. Amir and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]
47. G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. W. M. van der Laak, B. van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” Med. Image Anal. 42(9), 60–88 (2017). [CrossRef]
48. T. Kirchner, J. Gröhl, and L. Maier-Hein, “Context encoding enables machine learning-based quantitative photoacoustics,” J. Biomed. Opt. 23(05), 1 (2018). [CrossRef]
49. B. Ghafaryasl, K. A. Vermeer, J. Kalkman, T. Callewaert, J. F. de Boer, and L. J. Van Vliet, “Analysis of attenuation coefficient estimation in Fourier-domain OCT of semi-infinite media,” Biomed. Opt. Express 11(11), 6093–6107 (2020). [CrossRef]
