Suppression of motion artifacts in intravascular photoacoustic image sequences

Zheng Sun; Zheng Sun; Jiejie Du; Jiejie Du

doi:10.1364/BOE.440975

1. Introduction

Photoacoustic tomography (PAT) is a rapidly developing imaging modality that allows structural and functional imaging of biological tissues with high penetration depth and high optical contrast [1]. Intravascular photoacoustic (IVPA) imaging is a typical catheter-based application of PAT. As a complementary tool to intravascular ultrasound (IVUS) for the assessment of atherosclerosis, it is capable of providing multi-scale anatomical, functional, and molecular information of vessels by combining noninvasive PAT with endoscopic detection. It enables potential applications in the diagnosis and interventional treatment of cardiovascular diseases [2,3].

Successful application of an IVPA system for in vivo intracoronary imaging is challenging because motion artifacts owing to voluntary motion and involuntary motion are inevitable. The voluntary motion such as body movement can be avoided by decreasing the image acquisition time and keeping the imaging object still during acquisition [4]. In whole-body PAT for small animals, the animal movement is significantly reduced through mechanical clamping or fixation [5–7]. This scheme is only applicable to specific scenarios because clamping is not suitable for all body parts. The involuntary motion is associated with heart beat, breathing, pulsating blood flow within the vascular lumen, and arterial vasomotion. It leads to an undesired motion between the ultrasonic transducer and vessel wall, which involves the transversal motion and longitudinal motion. The transversal motion produces the shift (translation and rotation) of the vessel structures from one slice to the other, while the longitudinal motion induces a proximal or distal displacement additional to the pullback. As a consequence, the vessel cross-sections are not equally spaced in the longitudinal direction along the pullback trajectory and the geometrical description of the vessel wall is misled. The motion artifacts are visible as saw-tooth shaped vessel wall boundaries in the volumetric images that are reconstructed from pullback recordings. The artifacts degrade image quality and subsequently make identification and interpretation of structures like vessels and atherosclerotic plaques difficult. Furthermore, the precision of quantitative assessment of tissue properties and 3-D vessel rendering is reduced.

As an imaging modality with a high frame rate, PAT enables rendering of 2-D or 3-D images via excitation of an entire volume with a single nanosecond laser pulse, which avoids motion artifacts in a single frame [8]. However, for studies of vessel morphology, plaque characterization, and other purposes requiring 3-D imagery, acquisition and analysis of multiple frames rendering multiple vascular cross-sections are required. Accordingly, motion correction is relevant in IVPA imaging of coronary arterial vessels.

Many efforts have been made to mitigate motion artifacts in PAT applications involving multi-frame data analysis. A simple and generally effective way to account for motions caused by heart beat or breathing is to gate the sequences according to an electrocardiogram (ECG) or respiratory triggering signal prospectively or retrospectively. On-line prospective gating activates data acquisition during the same phase of each cycle by employing a triggering or synchronization scheme. Off-line retrospective gating captures images or raw signals continuously in the cycles and records ECG or respiratory waveforms simultaneously. After acquisition, images or signals that are collected at the same temporal location (cardiac/respiratory phase) are selected according to the ECG or respiratory signals. This technique is essential in many established imaging systems that need to locate objects accurately [9–12]. For instance, image-based gating has been routinely applied to suppress motion artifacts present in continuous pullback IVUS and intravascular optical coherence tomography (IV-OCT) image sequences [13–17]. One frame per cycle is extracted to form a subsequence by tracking the cyclic change of the image intensity or the vascular lumen contour across the entire pullback. In experiments of photoacoustic imaging (PAI), retrospective gating has been adopted to suppress respiratory motion artifacts in a whole-body PAT system for small animals [18,19]. It is implemented in two ways, hardware or software. The hardware respiratory gating is to collect raw PA signals when the imaging object is breathing freely, simultaneously monitoring the respiratory waveform with an external equipment. The collected PA signals are aligned according to the respiratory phases. Finally, the images in a complete respiratory cycle are reconstructed based on the PA measurements in the same respiratory phases [18]. This method needs breathing training for the imaging object because quick or uneven breathing results in the reduced accuracy of triggering. The software respiratory gating is to extract a signal that implies the respiratory phases from the PA measurements. The motion is subsequently corrected based on the signal prior to image reconstruction [19]. When the pause between two breaths exceeds the length of a single breath, it is necessary to establish a standard to distinguish between a dynamic and a static frame according to the prior knowledge regarding respiratory characteristics.

Frame motion compensation (FMC) and inter-frame motion compensation (IFMC) have been used to reduce motion artifacts in frame-averaged PA images [20]. Both methods determine motion vectors between frames by block matching and three-step search. FMC corrects the error from motion artifacts by comparing past images with the most recent image as the reference image. However, large motion might not be detected because the motion vector accumulating with the time is easy to exceed the detectable range. IFMC overcomes the limitation of FMC by comparing two consecutive frames. The accuracy of both methods highly depends on the search step that is difficult to be properly selected in an automatic manner.

In spectroscopic PAT of the heart, images are acquired under the multiple wavelengths excitation at separate time points. The cardiac motion during this process leads to blurring in the images. Motion clustering has been demonstrated effective in reducing such motion blurring [21]. A clustering algorithm such as k-means is utilized to separate a sequence of single-pulse images at multiple excitation wavelengths into clusters corresponding to different stages of the cardiac cycle. The number of clusters should be selected properly based on the severity of blurring, signal-to-noise ratio (SNR), and the performance of the clustering algorithm.

In addition, the model-based scheme has been demonstrated efficient in reconstructing PAT images with high quality besides the analytical inversion formulations such as back-projection (BP) [22]. A model physically describing the forward acoustic problem is established, which outputs the theoretical acoustic pressure induced by optical absorbers in tissues. The initial pressure or optical energy deposition is recovered by iteratively minimizing the error between the measured and theoretical acoustic pressure calculated by the forward model. The scheme is capable of including all linear effects in the forward model. By incorporating the motion into the forward model, estimation of the motion parameters and reconstruction of the desired images can be implemented simultaneously [23].

In this paper, we propose a novel method for suppressing motion artifacts associated with the cardiac cycle in IVPA pullback volumes. To the best of our knowledge this is the first paper focusing on motion artifact correction for volumetric IVPA images. The raw PA signals collected by the ultrasonic detector in successive vascular cross-sections are grouped into subsets that correspond to similar phases in the cardiac cycles by clustering. The sequential images representing the initial pressure distribution in vascular cross-sections are reconstructed based on the selected frames of signals by time reversal (TR). We tested our method on simulation data. According to the results, we analyzed the influence of the objective function threshold in the clustering algorithm on image reconstruction. In addition, we conducted experiments comparing our method with the image-based gating to demonstrate the superiority of our method in suppressing IVPA motion artifacts.

The remainder of this paper is organized as follows. Sect. 2 depicts the presented method in detail. Sect. 3 provides the demonstration and comparison results. Sect. 4 gives related discussions and Sect. 5 concludes the paper with a summary.

2. Method

2.1 Principle of IVPA imaging

As illustrated in Fig. 1, the procedure of IVPA imaging is analogous to IVUS, where a catheter is inserted into the vascular lumen for imaging and pushed to the distal end under the direction of X-ray angiography. During pullback of the catheter, the probe mounted on its tip emits short pulsed laser (∼ns) that is absorbed by the surrounding tissues, leading to a temperature rise. Subsequently, the thermo-elastic effect induces a pressure rise that is proportional to the optical energy deposition. The pressure rise propagates as wide-band (∼MHz) ultrasonic waves, that is, PA waves, to the tissue surface. The ultrasonic transducers on the probe scan the surrounding tissues circumferentially, collecting the photoacoustically generated pressure along a circular trajectory parallel to the imaging plane that is perpendicular to the catheter. The ultrasonic detector is idealized as a point-like detector by ignoring its aperture effect. A full-view (360°) scanning provides the measurements in M angles, {θ₁,θ₂,…,θ_M}, where N points are sampled on each scanning radius. Finally, images representing the spatially varying optical energy deposition or initial pressure in vascular cross-sections are reconstructed from the measured PA signals by acoustic inversion [24,25]. Moreover, the solution to optical inversion enables quantitative imaging by recovering optical properties (absorption coefficient μ_a and scattering coefficient μ_s), thermoelastic coefficient (Grüneisen coefficient), and functional properties including blood oxygenation (sO₂) and the concentration of chromophores from optical energy deposition or PA measurements [26,27].

Fig. 1. Schematic diagram of IVPA imaging. (a) Longitudinal view (L-view) of a vessel segment to be imaged; (b) Transversal view of an imaging plane; (c) Stacked transversal images in a temporal order; (d) Generation of cardiac cycle-dependent motion artifacts in a time-axis view of pullback volumes.

Download Full Size | PDF

During acquisition of PA signals in successive cross-sections, the catheter remains in the center of the transversal imaging plane as sketched in Fig. 1(d). The cardiac cycle-dependent motion causes misalignment of successive slices along the sequence of cross-sectional images, that is, B-scan images in Cartesian coordinates as well as saw-tooth shaped vessel wall boundaries in time-axis views of IVPA pullback volumes. Moreover, the variations of the luminal shape in the vertical direction are prominent than in the horizontal direction.

2.2 Motion suppression by signal clustering

Suppose that a single pullback produces W slices to be assessed. Each slice is discretized into M×N sampling locations. The PA signals collected in the kth slice are recorded in a M×N-dimensional matrix,

(1)$${{\boldsymbol P}_k} = \left[ {\begin{array}{cccc} {{{\boldsymbol p}_{11}}}&{{{\boldsymbol p}_{12}}}& \cdots &{{{\boldsymbol p}_{1N}}}\\ {{{\boldsymbol p}_{21}}}&{{{\boldsymbol p}_{22}}}& \cdots &{{{\boldsymbol p}_{2N}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{{\boldsymbol p}_{M1}}}&{{{\boldsymbol p}_{M2}}}& \cdots &{{{\boldsymbol p}_{MN}}} \end{array}} \right], $$

where k = 1,2,…,W, p_ij denotes a pressure vector collected at the jth location in the ith measuring angle, i = 1,2,…,M, and j = 1,2,…,N. A matrix P constitutes a signal frame and a single pullback produces W frames in total. The correlation matrix of these signal frames is obtained by

(2)$${\boldsymbol C} = \left[ {\begin{array}{cccc} {{\rho_{11}}}&{{\rho_{12}}}& \cdots &{{\rho_{1W}}}\\ {{\rho_{21}}}&{{\rho_{22}}}& \cdots &{{\rho_{2W}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{\rho_{W1}}}&{{\rho_{W2}}}& \cdots &{{\rho_{WW}}} \end{array}} \right], $$

where

(3)$${\rho _{ij}} = \frac{{\textrm{|Cov}({{\boldsymbol P}_i},{{\boldsymbol P}_j})|}}{{\sqrt {D({{\boldsymbol P}_i})} \cdot \sqrt {D({{\boldsymbol P}_j})} }}. $$

Here, i and j range from 1 to W, P_i and P_j denote, respectively, the ith and jth frame, ρ_ij is the correlation coefficient between P_i and P_j, D(P_i) and D(P_j) are the variances of P_i and P_j, and Cov(P_i, P_j) is the covariance of P_i and P_j. The correlation coefficients in C are rearranged into a 1-D array row by row, that is, {ρ₁₁, ρ₁₂,…, ρ_1W, ρ₂₁, ρ₂₂, …, ρ_2W, …, ρ_W₁, ρ_W₂, …, ρ_WW}. For simplicity, it is denoted as a data set, F = {f₁, f₂, f₃, …, f_Q}, where Q = W×W.

The entries in F are grouped into subsets corresponding to the similar phases in a cardiac cycle by affinity propagation (AP) clustering [28,29]. The detailed steps are as follows.

Step 1. Initialization.

The iteration times is initialized as t = 0. The responsibility matrix R_t= [r_t(i, k)]_Q_×_Q and availability matrix A_t= [a_t(i, k)]_Q_×_Q are initialized as zero matrices, where r₀(i, k) = a₀(i, k) = 0. r_t(i, k) denotes the responsibility sent from f_i to candidate cluster center f_k at the tth iteration, reflecting the accumulated evidence for how well-suited f_k is to serve as the cluster center for f_i. a_t(i, k) denotes the availability sent from candidate cluster center f_k to f_i at the tth iteration, reflecting the accumulated evidence for how appropriate it would be for f_i to choose f_k as its cluster center [28].

Step 2. Calculation of the similarity matrix.

The similarity s_t(i, k) is determined by

(4)$${s_t}(i,k) ={-} {|{{f_i} - {f_k}} |^2}, $$

which indicates how well f_k is suited to be the cluster center for f_i. The similarity matrix S_t at the tth iteration is a collection of real-valued similarities among all data points in F.

Step 3. Updating of responsibility matrix and availability matrix.

The responsibility matrix and availability matrix are updated as

(5)$${r_{t + 1}}(i,k) = \left\{ {\begin{array}{ll} {{s_t}(i,k) - \mathop {\max }\limits_{_{k^{\prime}\textrm{s}\textrm{.t}\textrm{.}k^{\prime} \ne k}} \{{{a_t}(i,k^{\prime}) + {r_t}(i,k^{\prime})} \}},&i \ne k\\ {{s_t}(i,k) - \mathop {\max }\limits_{_{k^{\prime}\textrm{s}\textrm{.t}\textrm{.}k^{\prime} \ne k}} \{{{s_t}(i,k^{\prime})} \}},&i = k \end{array}} \right.$$

and

(6)$${a_{t + 1}}(i,k) = \left\{ {\begin{array}{ll} {\min \left\{ {0,{r_{t + 1}}(k,k) + \sum\limits_{i^{\prime}\textrm{s}\textrm{.t}\textrm{.}i^{\prime} \notin \{ i,k\} } {\max ({0,{r_{t + 1}}(i^{\prime},k)} )} } \right\}},&i \ne k\\ {\sum\limits_{i^{\prime}\textrm{s}\textrm{.t}\textrm{.}i^{\prime} \ne k} {\max } \{{0,{r_{t + 1}}(i^{\prime},k)} \}},&i = k \end{array}} \right., $$

where i’, k’, i, and k range from 1 to Q, and r_t₊₁(i, k) and a_t₊₁(i, k) are, respectively, the responsibility and availability between f_i and f_k at the (t+1)th iteration.

Step 4. Attenuation of the responsibilities and availabilities.

To avoid numerical instability, the responsibilities and availabilities are attenuated by a damping factor,

(7)$$\left\{ {\begin{array}{l} {{{\hat{r}}_{t + 1}}(i,k) = \lambda {r_t}(i,k) + (1 - \lambda ){r_{t + 1}}(i,k)}\\ {{{\hat{a}}_{t + 1}}(i,k) = \lambda {a_t}(i,k) + (1 - \lambda ){a_{t + 1}}(i,k)} \end{array}} \right., $$

where λ denotes the damping factor ranging in [0.5,1]. In this study, we set it as 0.5.

Step 5. Determination of clustering results.

The sum of the updated responsibility and availability is calculated,

(8)$$e = {\hat{a}_{t + 1}}(i,k) + {\hat{r}_{t + 1}}(i,k). $$

When e reaches its maximum, f_i is the cluster center of f_k if i = k; otherwise, f_k is the cluster center of f_i.

Step 6. Determination of termination of the iteration.

Whether the iteration is terminated or not is determined by the following objective function,

(9)$$J\textrm{ = }\sum\limits_{f \in {{\boldsymbol h}_i},i = 1}^H {|f - {g_i}{|^2}}, $$

where f denotes a sample in the cluster, h_i is the cluster centered on g_i, and H is the number of clusters. If J exceeds the threshold, the iteration is terminated and accordingly the clustering results are the output; otherwise, let t←t+1 and return to Step 2.

After clustering, the frames of signals collected in the same cardiac phases are selected. Finally, the images representing the initial pressure distribution in the vascular cross-sections are reconstructed from these frames of signals with the TR reconstruction approach [30]. The complete process of the proposed method is illustrated in Fig. 2.

Fig. 2. Flowchart of complete process of our method.

Download Full Size | PDF

2.3 Performance evaluation

We validated our method on simulation data by considering the lack of continuous pullback volumetric data that was acquired on in vivo dynamic vessels owing to the limitation of experimental conditions and set-up. We implemented the method by programming using MATLAB (R2018a, The MathWorks, Inc., Natick, Massachusetts) on a laptop configured with a 2.5 GHz Intel Core i5-10300H CPU, 8 GB RAM, and Windows 10 64bits as the operating system.

2.3.1 Simulated image preparation

We constructed computer-generated phantoms that mimic coronary arterial vessels having different tissue types. Figure 3 shows representative examples of vascular cross-sections. Coronary arterial vessels follow the cyclic dynamics of the heart in the cardiac cycles, resulting in the periodical variations in the cross sectional area of the vascular lumen. For each phantom, we generated successive cross-sections along the long axis of the lumen corresponding to different time points in the cycles based on the periodic change of the luminal area. We assumed that the first frame in a pullback sequence is acquired at end-diastole at which moment the luminal area reaches its minimum. Thus, we determined the luminal area at a time-point n in the sequence by [31]

(10)$$S(n) = \left\{ {\begin{array}{lr} {A\sin (\mathrm{\pi }Rn)\exp ( - \mathrm{\pi }\alpha n) + {S_0},}&{\textrm{0 < }n < 1/R} \\ {S(n - 1/R)},&{n \ge 1/R} \end{array}} \right., $$

where A is a constant controlling the luminal area, R is the heart rate in Hz (beats per second), α is a constant determining the time when the luminal area reaches its maximum, and S₀ is the minimal luminal area at end-diastole. We obtained the vessel wall contour at time-point n+1 by expanding outward or contracting inward the contour at time-point n which is represented with discrete points ${{\boldsymbol V}_{i,n}} = ({l_{i,n}},{\theta _{i,n}})$. The polar coordinate of a point in the contour at time-point n+1, ${{\boldsymbol V}_{i,n + 1}} = ({l_{i,n + 1}},{\theta _{i,n + 1}})$, is determined by

(11)$$\left\{ {\begin{array}{l} {{l_{i,n + 1}} = {l_{i,n}}\lambda \sqrt {S(n + 1)/S(n)} }\\ {{\theta_{i,n + 1}} = {\theta_{i,n}}} \end{array}} \right., $$

where l denotes the polar radius, θ denotes the polar angle, and λ is a proportional factor controlling the extent to which the contour expands or contracts. The value of λ depends on the tissue type. We set λ = 1 for the lumen, λ < 1 for the calcified plaques owing to their poorer elasticity than the lumen, and λ > 1 for the fibrosis-lipid plaques owing to their better elasticity than the lumen. Accordingly, the cross-sectional model at each time-point in the cycle is automatically generated from the one at end-diastole.

We simulated the sequential frames of spatially varying PA signals in vascular cross-sections by inputting the generated cross-sectional geometrical models into our previously developed endoscopic PAT simulation platform [31]. Table 1 provides the optical and acoustic property parameters of the vessel phantoms defined by referring to the histological findings [32,33]. The speed of sound and density of each tissue type follow Gaussian distributions based on the values shown in the table.

Table 1. Parameters of optical and acoustic properties of vessel phantoms for forward IVPA simulation

View Table | View all tables in this article

Considering that motion artifacts associated with the cardiac cycle are not visually prominent in successive transversal images, we utilized time-axis views, that is, L-views, of the pullback volumes to facilitate analyzing motion artifacts. We obtained vertical and horizontal L-views as illustrated in Fig. 4.

2.3.2 Figures of Merit

We utilized the dissimilarity matrix (DM), average dissimilarity (AD), and average inter-frame dissimilarity (AIFD) as quantitative metrics to evaluate the quality of reconstructed image sequences. For a W-frame sequence {I₁,I₂,…,I_W}, a W×W-dimensional DM, that is, D_W_×W= [d_i_,j], is constructed by pairwise comparison of the images [17],

(12)$${d_{i,j}} = 1 - \frac{{\sum\limits_{k = 1}^{Width} {\sum\limits_{l = 1}^{Height} {|{{I_i}({k,l} )- {\mu_i}} |\cdot |{{I_j}({k,l} )- {\mu_j}} |} } }}{{\sqrt {\sum\limits_{k = 1}^{Width} {\sum\limits_{l = 1}^{Height} {{{[{{I_i}({k,l} )- {\mu_i}} ]}^2}} \sum\limits_{k = 1}^{Width} {\sum\limits_{l = 1}^{Height} {{{[{{I_j}({k,l} )- {\mu_j}} ]}^2}} } } } }}, $$

where i and j range from 1 to W, I_i and I_j are, respectively, the ith and jth image frame of dimensions Width×Height, d_i_,j is the dissimilarity between I_i and I_j, I_i(k,l) and I_j(k,l) are the gray-levels of the pixel (k, l) in I_i and I_j, respectively, and μ_i and μ_j are the average gray-levels of I_i and I_j. d_i_,j ranges in the interval [0,1] and d_i_,j= d_j_,i. A smaller element in the DM represents a frame pair that differ less in their appearances than more different frames.

AD and AIFD are, respectively, defined as [17]

(13)$$D(k )= \frac{1}{{W - k}}\sum\limits_{m = 1}^{W - k} {{d_{m,m + k}}}$$

and

(14)$$D = \frac{1}{W}\sum\limits_{i = 1}^W {\sum\limits_{j = 1}^W {{d_{i,j}}} }, $$

where k = 0,1,…,W‒1, D(k) denotes the AD between two frames with the interval of k frames, D(0) = 0, and d_m_,m+k is the dissimilarity between frames I_m and I_m_+k. Lower AD and AIFD indicate higher similarity between frames.

3. Results

3.1 Results of image reconstruction

In demonstration experiments, we generated successive 300 cross-sections corresponding to different time-points in the cardiac cycles for each phantom. We set α = 0.5 and R = 1.2 Hz, thereby the length of a cardiac cycle was ${1 / R} \approx 0.83\textrm{ }\textrm{s}$. The frame rate was set as 24 fps by referring to [34], that is, 20 frames per cycle, therefore 300 frames covered about 15 cycles. Figure 5 shows the transversal images obtained by forward simulation for the four cross-sections shown in Fig. 3. Figure 6 shows examples of transversal images representing the initial pressure distribution reconstructed from the simulated PA signals by the conventional TR algorithm without motion correction. In the figure, the shifts between successive cross-sections, that is, translation and rotation of the structures from one image to the other, can be observed. Figure 7 shows the L-views of four 300-frame sequences reconstructed by the conventional TR (non-gated sequence) and our method (gated sequence), respectively. In our method, the objective function threshold of clustering was set as 0.3. In the figure, motion artifacts associated with the cardiac cycle are observed as the saw-tooth shaped appearance of the vessel wall. Moreover, the artifacts are more prominent in the vertical L-views than in the horizontal ones. Note that the L-views of the gated sequences exhibit a significant enhancement of the visualization with the smoothed vessel wall boundaries as opposed to the non-gated sequences. However, there is an apparent loss of resolution due to gating, which is a known trade-off of the gating process [11].

Fig. 3. Geometry of vascular cross-sections which are numbered as I, II, III, and IV from left to right.

Download Full Size | PDF

Fig. 4. Schematic diagram of L-views of IVPA pullback volumes.

Download Full Size | PDF

Fig. 5. Simulated transversal images of four vessel cross-sections shown in Fig. 3. (a) Images representing the normalized optical energy deposition; (b) Images representing the normalized PA signals reaching the detector.

Download Full Size | PDF

Fig. 6. Images randomly selected from the pullback sequences of transversal images representing the initial pressure distribution, which are reconstructed directly from the simulated PA signals by TR without motion suppression. (a) Phantom I; (b) Phantom II; (c) Phantom III; (d) Phantom IV.

Download Full Size | PDF

Fig. 7. The gating results of the proposed method for four IVPA sequences. (a) Vertical L-views; (b) Horizontal L-views. The gated sequences are artificially stretched to match the same physical lengths as the non-gated ones in order to facilitate comparison due to the fact that there are much more images/mm present in the non-gated sequences, while the gated sequences have much fewer images due to the gating.

Download Full Size | PDF

Figure 8, Fig. 9 and Table 2 provide the results of the evaluation metrics obtained from the non-gated and gated image sequences. In Fig. 8, the DMs are displayed as grayscale images which exhibit periodic structures. Obviously, the overall brightness of the visualization of DMs is reduced significantly, while both AD and AIFD are decreased after motion correction, suggesting reduced dissimilarity between frames in the gating subsets.

Fig. 8. Visualization of DMs obtained from the non-gated (left column) and gated (right column) image sequences. (a) Phantom I; (b) Phantom II; (c) Phantom III; (d) Phantom IV

Download Full Size | PDF

Fig. 9. AD functions with respect to frame intervals obtained from the non-gated and gated image sequences for (a) phantom I, (b) II, (c) III, and (d) IV. There are 60, 65, 80, and 72 frames in four gated sequences, respectively.

Download Full Size | PDF

Table 2. AIFDs of the image sequences before and after motion correction.

View Table | View all tables in this article

3.2 Influence of objective function threshold

The objective function threshold in Eq. (9) determines whether the iteration of clustering is terminated or not. We set it as 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8 respectively while other conditions remained unchanged to investigate its influence on the performance of the proposed method. Table 3 provides the results of phantom I, revealing that a lower threshold leads to better results of motion suppression. However, a lower threshold leads to more iterations. Therefore, a trade-off between the time performance and the reconstruction accuracy should be made.

Table 3. AIFDs of the image sequences for phantom I in the case of different thresholds of the objective function.

View Table | View all tables in this article

3.3 Results of comparison with the image-based gating

We compared the proposed method with an image-based gating approach that was developed for IVUS [17]. The results shown in Fig. 10 reveal that the motion artifacts in the image sequences processed by image-based gating are still prominent compared with those obtained by our method. This conclusion can also be drawn from the results of evaluation metrics provided in Fig. 11. In addition, we recorded the run-time of both methods as provided in Table 4. The run-time for the image-based gating includes the time cost in reconstructing the entire image sequence by TR besides the time cost in gating. The results indicate that our method is more computationally efficient than the image-based gating. This is because our method groups raw PA signals prior to image reconstruction, avoiding reconstruction and post-processing of those non-gating frames. In contrast, the image-based gating corrects motion artifacts via a post processing procedure over the complete pullback recording. It requires to reconstruct the successive images from the raw signals collected in each slice along the pullback trajectory. This procedure is computationally expensive compared with the off-line gating itself.

Fig. 10. L-views of the IVPA sequence of phantom I processed with our method and the image-based gating method. (a) Non-gated sequence; (b) Gated sequence with our method; (c) Gated sequence with the image-based method.

Download Full Size | PDF

Fig. 11. DMs and ADs of the reconstructed image sequences for phantom I. There are 60 and 55 frames in the gated sequences obtained by our method and image-based gating, respectively. (a) Visualization of the DMs obtained from non-gated sequence (left), gated sequence by our method (middle), and gated sequence by image-based gating (right); (b) AD with respect to frame interval.

Download Full Size | PDF

Table 4. Quantitative metrics of the image sequences of phantom I processed by two methods.

View Table | View all tables in this article

4. Discussion

4.1 Superiority of our method to other motion correction methods

As summarized in Introduction, the existing strategies for motion correction in PAT involving multi-frame analysis include FMC, motion clustering, model-based reconstruction, and gating. FMC or IFMC aims to reduce motion artifacts in frame-averaged PA images and enhance the SNR of deep tissue imaging which utilizes the regulative energy of a commonly-used Nd:YAG laser. Motion clustering focuses on alleviating motion blurring which is observed in spectroscopic PA images acquired under the multiple wavelengths excitation at different time-points. Both strategies are applicable to the scenarios where multiple images are acquired for the same target with motion. It is infeasible to apply them to correct motion artifacts associated with the cardiac cycle in IVPA pullback volumetric images that are acquired at different locations along the catheter pullback path. The model-based motion correction is to solve the forward problem iteratively in nature, where the forward operator is calculated repeatedly. This procedure is computationally burdensome, which hinders its application in real-time, high-resolution, and large-volume imaging including IVPA.

Gating techniques have been commonly utilized in cardiac imaging to reduce motion artifacts owing to heart beat or breathing. However, prospective gating requires a special ECG or respiratory triggering device, causing the added complexity, long setup times, and prolonged acquisition procedure compared with the continuous non-gating acquisition. The accuracy of retrospective gating depends on the gating signals that are extracted from the raw measurements or images. The image-based gating relies on post processing of images, therefore it requires to reconstruct all images in a dynamic sequence prior to gating. A single IVPA pullback sequence that is acquired on a vessel segment of 10 mm long contains about 400 frames in the case of constant pullback (0.5 mm/sec) and a frame rate of 20 fps. The overload of image data reduces the efficiency of the image-based algorithm. Moreover, the inevitable loss of information regarding the vascular structures and properties in the procedure of image reconstruction is an important issue that should be taken into account. Our method is based on a retrospective gating scheme by selecting those frames of PA signals that are acquired in the same cardiac phase. It differs from the state-of-art image-based gating schemes in that the latter subtracts cardiac dynamics based on grayscale images themselves, whereas our method suppresses motion prior to image reconstruction avoiding reconstructing the non-gating frames.

4.2 Limitations of our method

Our method enables full suppression of motion artifacts associated with the cardiac cycle with a relatively low computation cost compared with the image-based gating. However, the loss of information and resolution is a known trade-off of the gating process. Only static sequences at certain cardiac phases are remained in the procedure of resampling signals or images, leading to the loss in the useful information regarding vascular structural and functional features. Moreover, all frames between systole and diastole are generally required to be analyzed for a continuous assessment of tissue elastic properties. Further study on the data integrity is desired in the future, so as to ensure the rendered image sequences contain complete information while enhancing the image quality.

The motion artifacts in in vivo sequential intravascular images are a combined result of various factors such as heartbeat motion, pulsating blood flow in the vascular lumen, arterial vasomotion, and catheter-based motion. The first three issues are related to the cardiac cycle. As a result, the trajectory of the catheter tip does not remain parallel to the lumen axis owing to the lateral movement of the catheter tip with respect to the lumen, leading to the shift between successive slices. In addition, an irregular deformation and pulsation of the vessel wall are related to heartbeat motion. Our method aims to suppress motion artifacts associated with the cardiac cycle. However, there remains a need to eliminate catheter-based motion artifacts, which is a common issue in catheter-based imaging systems, such as catheter bending, longitudinal oscillation of the catheter, and nonuniform rotation distortion (NURD). The longitudinal oscillation of the catheter within the vascular lumen results in the repeated sampling at a acquisition position. NURD is inherently present in intravascular imaging systems with rotary-pullback catheters owing to mechanical friction between the catheter torque cable and sheath [35–37]. Recently, the distal scanning endoscopes with miniature micromotors have been designed to deal with artifacts associated with bending and NURD of proximally rotated catheters [38]. However, it is a challenging and expensive task to miniaturize high-speed motors since the relatively large size of the motor limits the catheter to access the vascular stenosis [36,39]. Consequently, the catheter-based motion artifacts might still exist in clinical scenarios, which degrade image quality and hinder identification of various tissue types as well as the quantification of tissue properties.

4.3 Future direction related to deep learning

In recent years, deep learning (DL) has been dominating in the field of medical imaging by significantly facilitating the performance of multiple tasks [40]. DL methodologies have gained interest as potential solutions for efficient processing and analysis of large datasets since they have outstanding advantages in image processing, identification, and interpretation over non-learning ones.

Chen et al. [41] reported a method for correcting motion artifacts in optical resolution photoacoustic microscopic (OR-PAM) images by convolutional neural network (CNN). To the best of our knowledge, it is the first study on motion correction by DL in the field of PAI. They constructed a CNN with which they post-processed the maximum amplitude projection (MAP) image of OR-PAM to alleviate motion artifacts. However, motion suppression for multi-frame PAT by DL has not been explored yet up to the writing of this paper. One of the core bottlenecks is the lack of reliable experimental training data. The size and quality of the data sets determine the efficacy of DL methodologies. Three sources for obtaining training data have been adopted currently, clinical in vivo data, phantom data, and simulation data. Previous study on reconstruction of a single PA image by DL mostly uses phantom data and simulation data for network training and verification. This is attributed to the lack of the ground truth information on the underlying tissue optical properties or the initial pressure distribution when acquiring experimental measurements [42]. Motion artifacts in IVPA imaging are unobservable in a single cross-sectional image, therefore multi-frame analysis is essential, which further improves the requirements to the amount of the training data. In addition, the underlying dynamics of the arterial vessels are generally unknown in in vivo experimental settings. Another challenging issue is that it is a laborious task to manufacture a dynamic phantom mimicking the arterial vessel in a complex clinical scenario. Computer simulation enables generation of a large amount of simulation data in a flexible way. However, it is still difficult to meet the requirements of the amount and variety of training data used in motion artifacts correction by DL. Specifically, the algorithms trained by the simulation data might fail in a clinical scenario because of the fact that the simulation data suffer from several shortcomings compared with the data distribution in reality, such as domain gap, sparsity or a selection bias [42]. By considering the similarity of IVUS and IVPA in their imaging principles and image contents, transfer learning may be a possible solution in the future. Namely, the motion correction algorithm is trained on a large data set of in vivo IVUS pullback studies, whereas a smaller experimental PA data set is utilized to finely adjust the neural network to the experimental data distribution.

5. Conclusion

In summary, motion artifacts removal is essential for stabilizing IVPA image sequences that suffer from cardiac motion. Our method is driven by raw PA signals alone that are acquired in successive slices by ultrasonic transducers during continuous pullback of the catheter. We grouped PA signals by clustering to select the signal frames that are collected in the same cardiac phase. The results demonstrate an improvement in the visualization of the L-mode views, where saw-tooth shaped vessel wall boundaries are considerably smoothed. The quantitative evaluation metrics including AD and AIFD can be increased by up to 50% after motion correction, indicating the reduction in the misalignment between vascular cross-sections. In addition, our method outperforms the image-based gating in correcting motion artifacts with a low computational burden.

Funding

National Natural Science Foundation of China (62071181).

Acknowledgments

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. W. Choi, D. Oh, and C. Kim, “Practical photoacoustic tomography: realistic limitations and technical solutions,” J. Appl. Phys. 127(23), 230903 (2020). [CrossRef]

2. S.S.S. Choi and A. Mandelis, “Review of the state of the art in cardiovascular endoscopy imaging of atherosclerosis using photoacoustic techniques with pulsed and continuous-wave optical excitations,” J. Biomed. Opt. 24(8), 080902 (2019). [CrossRef]

3. Y. Li, J. Chen, and Z. Chen, “Multimodal intravascular imaging technology for characterization of atherosclerosis,” J. Innov. Opt. Health Sci. 13(01), 2030001 (2020). [CrossRef]

4. H. Zhao, N. Chen, T. Li, J. Zhang, R. Lin, X. Gong, L. Song, Z. Liu, and C. Liu, “Motion correction in optical resolution photoacoustic microscopy,” IEEE Trans. Med. Imaging 38(9), 2139–2150 (2019). [CrossRef]

5. H.P. Brecht, R. Su, M. Fronheiser, S.A. Ermilov, A. Conjusteau, and A.A. Oraevsky, “Whole-body three-dimensional optoacoustic tomography system for small animals,” J. Biomed. Opt. 14(6), 064007 (2009). [CrossRef]

6. R.B. Lam, R.A. Kruger, D.R. Reinecke, S.P. Delrio, M.M. Thornton, P. Picot, and T.G. Morgan, “Dynamic optical angiography of mouse anatomy using radial projections,” Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing, 23 Feb., vol. 7564, pp. 756405 (2010).

7. R. Lin, J. Chen, H. Wang, M. Yan, W. Zheng, and L. Song, “Longitudinal label-free optical-resolution photoacoustic microscopy of tumor angiogenesis in vivo,” Quant. Imag. Med. Surg. 5(1), 23–29 (2015). [CrossRef]

8. X. Deán-Ben, S. Gottschalk, B. Mc Larney, S. Shoham, and D. Razansky, “Advanced optoacoustic methods for multiscale imaging of in vivo dynamics,” Chem. Soc. Rev. 46(8), 2158–2198 (2017). [CrossRef]

9. M.H. Moghari, A. Barthur, M.E. Amaral, T. Geva, and A.J. Powell, “Free-breathing whole-heart 3D cine magnetic resonance imaging with prospective respiratory motion compensation,” Magn. Reson. Med. 80(1), 181–189 (2018). [CrossRef]

10. Y. Lu, K. Fontaine, T. Mulnix, J.A. Onofrey, S. Ren, V. Panin, J. Jones, M.E. Casey, R. Barnett, P. Kench, and R. Fulton, “Respiratory motion compensation for PET/CT with motion information derived from matched attenuation-corrected gated PET data,” J. Nucl. Med. 59(9), 1480–1486 (2018). [CrossRef]

11. N. Torbati, A. Ayatollahi, and P. Sadeghipour, “Image-based gating of intravascular ultrasound sequences using the phase information of dual-tree complex wavelet transform coefficients,” IEEE Trans. Med. Imaging 38(12), 2785–2795 (2019). [CrossRef]

12. R. Bajaj, X. Huang, Y. Kilic, A. Jain, A. Ramasamy, R. Torii, J. Moon, T. Koh, T. Crake, M.K. Parker, V. Tufaro, P.W. Serruys, F. Pugliese, A. Mathur, A. Baumbach, J. Dijkstra, Q. Zhang, and C.V. Bourantas, “A deep learning methodology for the automated detection of end-diastolic frames in intravascular ultrasound images,” Int J Cardiovasc Imaging 37(6), 1825–1837 (2021). [CrossRef]

13. Z. Sun and Q. Yan, “An off-line gating method for suppressing motion artifacts in ICUS sequence,” Comput. Biol. Med. 40(11-12), 860–868 (2010). [CrossRef]

14. Z. Sun and M. Li, “Suppression of cardiac motion artifacts in sequential intracoronary optical coherence images,” J. Med. Imag. Health In. 6(7), 1787–1793 (2016). [CrossRef]

15. A. Hernàndez-Sabaté, D. Gil, J. Garcia-Barnés, and E. Martí, “Image-based cardiac phase retrieval in intravascular ultrasound sequences,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 58(1), 60–72 (2011). [CrossRef]

16. S.K. Nadkarni, D.R. Boughner, and A. Fenster, “Image-based cardiac gating for three-dimensional intravascular ultrasound,” Ultrasound in Medicine & Biology 31(1), 53–63 (2005). [CrossRef]

17. S.M. O’Malley, J.F. Granada, S. Carlier, M. Naghavi, and I.A. Kakadiaris, “Image-based gating of intravascular ultrasound pullback sequences,” IEEE Trans. Inform. Technol. Biomed. 12(3), 299–306 (2008). [CrossRef]

18. J. Xia, W.Y. Chen, K. Maslov, M.A. Anastasio, and L.V. Wang, “Retrospective respiration-gated whole-body photoacoustic computed tomography of mice,” J. Biomed. Opt. 19(1), 16003 (2014). [CrossRef]

19. A. Ron, N. Davoudi, X.L. Deán-Ben, and D. Razansky, “Self-gated respiratory motion rejection for optoacoustic tomography,” Appl. Sci. 9(13), 2737 (2019). [CrossRef]

20. M. Kim, J. Kang, J.H. Chang, T.K. Song, and Y. Yoo, “Image quality improvement based on inter-frame motion compensation for photoacoustic imaging: a preliminary study,” Proceedings of 2013 IEEE International Ultrasonics Symposium (IUS), Prague, Czech Republic, 21-25 July, pp. 1527-1531 (2013).

21. A. Taruttis, J. Claussen, D. Razansky, and V. Ntziachristos, “Motion clustering for deblurring multispectral optoacoustic tomography images of the mouse heart,” J. Biomed. Opt. 17(1), 016009 (2012). [CrossRef]

22. M. Xu and L.V. Wang, “Time-domain reconstruction for thermoacoustic tomography in a spherical geometry,” IEEE Trans. Med. Imaging 21(7), 814–822 (2002). [CrossRef]

23. J. Chung and L. Nguyen, “Motion estimation and correction in photoacoustic tomographic reconstruction,” SIAM J. Imaging Sci. 10(1), 216–242 (2017). [CrossRef]

24. J. Poudel, L. Yang, and M.A. Anastasio, “A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography,” Phys. Med. Biol. 64(14), 14TR01 (2019). [CrossRef]

25. R. Manwar, M. Zafar, and Q. Xu, “Signal and image processing in biomedical photoacoustic imaging: a review,” Optics 2(1), 1–24 (2021). [CrossRef]

26. Y. Liu, M. Sun, T. Liu, Y. Ma, D. Hu, C. Li, and N. Feng, “Quantitative reconstruction of absorption coefficients for photoacoustic tomography,” Appl. Sci. 9(6), 1187 (2019). [CrossRef]

27. T. Chen, T. Lu, S. Song, S. Miao, F. Gao, and J. Li, “A deep learning method based on U-Net for quantitative photoacoustic imaging,” Proceedings of SPIE International Conference on Photons Plus Ultrasound: Imaging and Sensing, vol.11240, pp.112403 V (2020).

28. S. Zhou and Z. Xu, “Automatic grayscale image segmentation based on affinity propagation clustering,” Pattern Anal Applic 23(1), 331–348 (2020). [CrossRef]

29. C.X. Sun, Y. Yang, H. Wang, and W. Wang, “A clustering approach for motif discovery in ChIP-Seq dataset,” Entropy 21(8), 802 (2019). [CrossRef]

30. Z. Sun, D. Han, and Y. Yuan, “2-D image reconstruction of photoacoustic endoscopic imaging based on time-reversal,” Comput. Biol. Med. 76, 60–68 (2016). [CrossRef]

31. Z. Sun, Y. Yuan, and D. Han, “A computer-based simulator for intravascular photoacoustic images,” Comput. Biol. Med. 81, 176–187 (2017). [CrossRef]

32. S.L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]

33. Y. Xu, H. Xue, and G. Hu, “Diameter measurements of coronary artery segments based on image analysis of X-ray angiograms,” Chinese J. Biomed. Eng. 26(6), 874–878 (2007).

34. D. Vanderlaan, A.B. Karpiouk, D. Yeager, and S. Emelianov, “Real-time intravascular ultrasound and photoacoustic imaging,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr 64(1), 141–149 (2017). [CrossRef]

35. J. Mavadia-Shukla, J. Zhang, K. Li, and X. Li, “Stick-slip non-uniform rotation distortion correction in distal scanning optical coherence tomography catheters,” J. Innov. Opt. Heal. Sci. 13(06), 2050030 (2020). [CrossRef]

36. E. Abouei, A.M.D. Lee, H. Pahlevaninezhad, G. Hohert, M. Cua, P. Lane, S. Lam, and C. MacAulay, “Correction of motion artifacts in endoscopic optical coherence tomography and autofluorescence images based on azimuthal en face image registration,” J. Biomed. Opt. 23(01), 1 (2018). [CrossRef]

37. O.O. Ahsen, H.C. Lee, M.G. Giacomelli, Z. Wang, K. Liang, T.H. Tsai, B. Potsaid, H. Mashimo, and J.G. Fujimoto, “Correction of rotational distortion for catheter-based en face OCT and OCT angiography,” Opt. Lett. 39(20), 5973–5976 (2014). [CrossRef]

38. J. Peng, L. Ma, X. Li, H. Tang, Y. Li, and S. Chen, “A novel synchronous micro motor for intravascular ultrasound imaging,” IEEE Trans. Biomed. Eng. 66(3), 802–809 (2019). [CrossRef]

39. F. Griese, S. Latus, M. Schluter, M. Graeser, M. Lutz, A. Schlaefer, and T. Knopp, “In-vitro MPI-guided IVOCT catheter tracking in real time for motion artifact compensation,” PLoS One 15(3), e0230821 (2020). [CrossRef]

40. H. Deng, H. Qiao, Q. Dai, and C. Ma, “Deep learning in photoacoustic imaging: a review,” J. Biomed. Opt. 26(04), 040901 (2021). [CrossRef]

41. X. Chen, W. Qi, and L. Xi, “Deep-learning-based motion-correction algorithm in optical resolution photoacoustic microscopy,” Vis. Comput. Ind. Biomed. Art 2(1), 12 (2019). [CrossRef]

42. J. Grohl, M. Schellenberg, K. Dreher, and L. Maier-Hein, “Deep learning for biomedical photoacoustic imaging: A review,” Photoacoustics 22, 100241 (2021). [CrossRef]

Tissue name	Tissue component	Refractive index	Reflectance	μ_a (cm^-1)	μ_s (cm^-1)	Anisotropy factor	Sound speed (m/s)	Density (kg/L)	Average radial thickness (mm)
Adventitia	Connective tissue	1.41	0.08	0.22	18	0.82	1610	1.14	0.25-0.45
Intima/ Media	Muscular tissue	1.42	0.05	0.18	18	0.82	1600	1.12	0.24 (Phantom I and II) 0.4 (Phantom III and IV)
Fibrous cap	Fibrous tissue	1.46	0.24	0.5	100	0.8	1630	1.2	0.16
Calcified plaque	Calcium	1.46	0.34	0.6	560	0.82	1650	1.24	0.2 (Phantom I and II) 0.45 (Phantom IV)
Lipid pool	Lipid	1.46	0.42	0.65	520	0.82	1650	1.22	0.33 (Phantom II) 0.15 (Phantom III)
Lipid-fibro plaque	Lipid- fibrosis	1.46	0.34	0.45	560	0.83	1560	1.17	0.3
Lumen	Blood	1.34	0.12	0.99	750	0.999	1540	1	3 (Phantom I) 2.5 (Phantom II and III) 2 (Phantom IV)

Vessel phantom	AIFD
Vessel phantom	Non-gated sequence	Gated sequence
I	0.101	0.010
II	0.028	0.009
III	0.020	0.011
IV	0.027	0.019

Threshold of objective function	Approximation of objective function	AIFD
Threshold of objective function	Approximation of objective function	Non-gated sequence	Gated sequence
1.8	1.681	0.101	0.101
1.5	1.437		0.085
1.2	1.098		0.070
0.9	0.834		0.053
0.6	0.562		0.050
0.3	0.289		0.010

Method	Length of gated sequence in frames	Run-time in seconds	AIFD
Method	Length of gated sequence in frames	Run-time in seconds	Non-gated sequence	Gated sequence
Our method	60	227.1	0.101	0.010
Image-based gating	55	539.5	0.101	0.036

Tissue name	Tissue component	Refractive index	Reflectance	μ_a (cm^-1)	μ_s (cm^-1)	Anisotropy factor	Sound speed (m/s)	Density (kg/L)	Average radial thickness (mm)
Adventitia	Connective tissue	1.41	0.08	0.22	18	0.82	1610	1.14	0.25-0.45
Intima/ Media	Muscular tissue	1.42	0.05	0.18	18	0.82	1600	1.12	0.24 (Phantom I and II) 0.4 (Phantom III and IV)
Fibrous cap	Fibrous tissue	1.46	0.24	0.5	100	0.8	1630	1.2	0.16
Calcified plaque	Calcium	1.46	0.34	0.6	560	0.82	1650	1.24	0.2 (Phantom I and II) 0.45 (Phantom IV)
Lipid pool	Lipid	1.46	0.42	0.65	520	0.82	1650	1.22	0.33 (Phantom II) 0.15 (Phantom III)
Lipid-fibro plaque	Lipid- fibrosis	1.46	0.34	0.45	560	0.83	1560	1.17	0.3
Lumen	Blood	1.34	0.12	0.99	750	0.999	1540	1	3 (Phantom I) 2.5 (Phantom II and III) 2 (Phantom IV)

Suppression of motion artifacts in intravascular photoacoustic image sequences

Abstract

1. Introduction

2. Method

2.1 Principle of IVPA imaging

2.2 Motion suppression by signal clustering

2.3 Performance evaluation

2.3.1 Simulated image preparation

2.3.2 Figures of Merit

3. Results

3.1 Results of image reconstruction

3.2 Influence of objective function threshold

3.3 Results of comparison with the image-based gating

4. Discussion

4.1 Superiority of our method to other motion correction methods

4.2 Limitations of our method

4.3 Future direction related to deep learning

5. Conclusion

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (11)

Tables (4)

Equations (14)

Biomedical Optics Express