Alignment of intravascular optical coherence

Alignment of intravascular optical coherence tomography movies affected by non-uniform rotation distortion G. van Soesta , J. G. Boscha and A. F. W. van der Steena,b a Dept.

of Biomedical Engineering, Thorax center, Erasmus MC, PO Box 2040, 3000CA, Rotterdam, The Netherlands; b Interuniversity Cardiology Institute of The Netherlands, Utrecht, The Netherlands. ABSTRACT

Endoscopic optical coherence tomography (OCT), and other imaging modalities that use a mechanically rotated probe, often suffer from image degradation due to non-uniform rotation distortion (NURD). In this paper we present a new method to align a sequence of images by globally optimizing the match between individual lines in subsequent frames. It uses dynamic programming to find a continuous path through a cost matrix that measures the similarity between regions of two frames being aligned. The path represents the angular mismatch corresponding to the NURD. The prime advantage of this novel approach compared to earlier work is the line-to-line continuity, which accurately captures slow intra-frame variations in rotational velocity of the probe. The algorithm is optimized using data from a clinically available intravascular OCT instrument in a realistic vessel phantom. Sensitivity of the performance to imaging and optimization parameters is investigated using a computational phantom. Finally, the algorithm’s efficacy is demonstrated on an in vivo recording inside a human coronary artery, exhibiting strong motion artifact. Keywords: intravascular imaging, NURD, motion compensation, non-uniform rotation, OCT, IVUS, image correction

1. INTRODUCTION OCT has been used for imaging of atherosclerotic plaques for more than a decade. It has become available for intracoronary clinical imaging relatively recently. Virtually all intravascular OCT1–3 scan the vessel with a rotating probe. A single imaging element is in the distal tip of a long catheter that is rotated proximally. Image quality is a critical factor, not only for diagnostic accuracy, but particularly for quantitative analysis of the acquired data. Enhanced or functional imaging techniques, such as OCT elastography,4 optical attenuation mapping,5 as well as inspection of deployed stent symmetry by OCT,6, 7 require highly accurate registration of the images that are analyzed. Among other factors, non-uniform catheter rotation (NURD) may degrade image quality. The rotation velocity of the distal imaging element is not necessarily constant, even if driven at a constant speed proximally. Variations occur particularly if tortuous bends need to be made to manoeuver the tip towards the imaged site. The number of A-lines (about 200–300) constituting a movie frame may span a rotation angle that is larger or smaller than 360◦ . We have recently8 developed a correction method for NURD in OCT, based on the retrieval of a continuous shift representing the angle deviation of the frame rotation. This is done using an algorithm based on dynamic programming: we find the globally cheapest path through a cost matrix representing the similarity between lines in subsequent frames. The parameters of the method are optimized on a computational phantom, with a known NURD, and on in vitro movies of a vessel phantom. The effectiveness of the procedure is demonstrated by the successful alignment of a movie acquired in vivo, in a human coronary artery. In commercial software, image stabilization for OCT sequences, compensating NURD and other motion artifacts, is performed using a global rotational block matching (GRBM) type technique. Send correspondence to G.v.S.; e-mail [email protected]

Coherence Domain Optical Methods and Optical Coherence Tomography in Biomedicine XII, edited by Joseph A. Izatt, James G. Fujimoto, Valery V. Tuchin, Proc. of SPIE Vol. 6847, 684721, (2008) 1605-7422/08/$18 · doi: 10.1117/12.761106 Proc. of SPIE Vol. 6847 684721-1 2008 SPIE Digital Library -- Subscriber Archive Copy

(a)

(b)

(c)

Figure 1. An SPIE logo was used as a test image for simulations. (a) Low-res logo upsampled to create gray-values; (b) frame from a sequence to which non-uniform rotation is applied; (c) with simulated speckle.

2. METHODS 2.1 NURD simulation In actual OCT movies, the NURD is unknown, so there is no ground truth. The ability of the algorithm to capture the angular deviation directly, one has to resort to simulations. An SPIE logo is used as a test image. A low-resolution version of the logo was upsampled to 512×512 pixels and subsequently converted to polar coordinates (300 angular lines × 256 radial pixels). A random deviation from linear angle sampling was generated, subject to a given standard deviation and correlation length, and applied to a sequence of 50 images. Five realizations of each configuration were tested to investigate the variation in optimal parameters. The influence of speckle on algorithm performance is studied by adding a random phase to every pixel in the test image sequence, followed by bandpass filtering. The absolute value of the resulting array accurately models speckle as encountered in coherent detection techniques.9 The sequence has to be treated line-by-line, in order to ensure uncorrelated speckle, as in intravascular OCT. See Fig. 1 for an example.

2.2 OCT Imaging and Materials OCT imaging is performed using a commercially available, clinically approved system (M2CV) with an ImageWire 2 catheter (both from Lightlab Imaging, Inc.). In vitro images were acquired in a polyvinyl alcohol cryogel vessel phantom, using Nlines = 312 lines per frame, Npix = 760 samples per line. The frame rate was 10 frames per second (fps). In vivo OCT movies were recorded in a male patient, age 51, in a healthy section of coronary artery. After inflation of an occlusion balloon the blood was cleared from the vessel using a saline flush. The frame rate was 15.6 fps and the frame size Nlines = 200, Npix = 744. The speckle composing the OCT image is undersampled laterally by a factor 2–3. Adjacent lines are therefore uncorrelated on the intensity level. We tested and validated the algorithm on several OCT movies, consisting of 12–51 frames.

2.3 NURD Correction Method The NURD correction procedure is illustrated in Fig. 2. It consists of six steps: 1. concatenation of frames to an Nframes · Nlines series of lines Aj ; 2. construction of a cost matrix; 3. resampling the cost matrix, to reduce noise and ensure proper connectivity constraints; 4. DP solution of the optimal path representing the NURD; 5. restoration of the path to the original dimensions and scaling of the problem; and 6. applying the correction.

Proc. of SPIE Vol. 6847 684721-2

(a)

j'-n

frame sequence Ai,j'

j'

j'+n

i-1

i

i+1

i+2

L2 norm

cost matrix

(b)

(c)

resample 1/m, s

resample m scale 1/s

i-3

(d)

i-2 angle (actual)

i-1 i

accum

+

P'i

align

Ai,j' ← Ai,j'+P'ij'

angle (assumed)

Figure 2. Schematic diagram presenting the NURD correction algorithm, see text for details. (a) Frame sequence in polar coordinates; a cost matrix is based on the L2 norm of the differences between lines in subsequent frames. (b) The optimal path (dotted line) is traced through the cost using dynamic programming. (c) Sketch of actual angle (unknown) vs. assumed angle (regularly spaced). The path reflects the deviation from a constant rotation velocity. (d) To align the data, the path is recursed and accumulated, resulting in the recipe for the line shifts in the gray box.

The cost matrix Ck,j is taken to be the L2 distance between two image lines j and j +Nlines +k; k = −n . . .+n. As the rotation deviation across a frame is assumed to be a small fraction of 360◦, it is desirable to be able to capture shifts of smaller than one line. This is achieved by resampling the cost function, in the sense that it is stretched (upsampled using a low-pass filtered interpolation algorithm) by a factor s in the vertical (shift) direction and compressed (downsampled) by a factor m in the horizontal (time) direction. The product r = m · s is the inverse of the effective fractional step size, which must be adjusted to match the dynamics of the problem. If r is too small, the algorithm takes relatively large steps and becomes very sensitive to noise. If r is large, the steps taken are too small and the NURD cannot be followed by the optimization. The relation between m and s on one hand and the NURD characteristics on the other will be studied with simulations. In this formulation, a continuous, horizontal path is determined from left to right through the cost matrix C, using a method called windowed dynamic time warping.10, 11 Briefly, an optimal path is defined as the one with the lowest total cost, the total cost of a path being sum of the costs of all steps. If we know the total cost gk,j−1 of the optimal partial paths up to each element of column j − 1 of the matrix, we can compute the optimal costs gk,j for all elements of the next column j: gk,j = Ck,j +

min

=k−1,k,k+1

{g,j−1 } .

(1)

Given that g1 = C1 , we can calculate cumulative costs for each following column, while keeping track of the optimal predecessor of each element. The element with the minimal cumulative cost in the last column is the end point of the optimal path. Backtracking, using the archived steps [k − 1, k, k + 1], the full path can be constructed. The continuous path P represents the slowly varying angle deviation induced by the slow speed variations; in fact Nlines (j+Nlines − j ) , (2) Pj = 360◦


shift (# lines)

20 10 0 -10 -20

5

10

15 frame number

20

25

30

Figure 3. Cost matrix after resampling at m = 5; s = 8 (grayscale values) and the optimized NURD correction path through this matrix (line). 8 6

angular shift (# lines)

4 2 0 −2 −4 −6 −8

simulated NURD correction 0

1000

2000

3000

4000 5000 line number

6000

7000

8000

9000

Figure 4. Simulated frame to frame shift in terms of line numbers (1 line corresponds to a 1.2◦ angle) (black line), compared to the optimal shifts retrieved by the NURD correction algorithm (gray dashed line). Near the top of the graph the difference between the two curves is shown.

if j denotes the deviation from the assumed angle at line j. An example of a resampled cost matrix from an actual OCT movie, recorded in vitro, with the optimized path, is shown in Fig. 3.

2.4 Analysis and parameter optimization To evaluate the performance of the NURD correction algorithm, we use the mean normalized standard deviation of the aligned data. In an actual experiment, no guaranteed NURD-free measurement is available, so there is no ground truth to compare with. Provided the movie is stationary, i.e. every image contains the same structure, the standard deviation of the aligned frames, normalized to the mean aligned frame, can be used as a measure for the quality of the alignment. The smaller the standard deviation is, the better the alignment. This assumption is tested using the simulated data. We compare the optimum found by the minimized normalized standard deviation, with the root mean square (rms) difference between optimized path Pj and the simulated input.

3. RESULTS 3.1 Simulation results The optimal resampling parameters for NURD correction are studied, depending on the characteristics of the NURD: angle standard deviation σ and correlation length Γ . An example of how the simulated deviation can be retrieved by the correction algorithm is shown in Fig. 4. The smallest root mean square (rms) difference between the simulated and optimized traces occurs for similar m · s as the smallest mean normalized standard deviation, validating the latter as a measure for alignment.


100 Γε = 300 Γ = 500

90

ε

Γε = 750 Γε = 1000

80

Γ = 1500 ε

70

〈m⋅s〉

60 50 40 30 20 10 0

0

50

100

150

200

250 2π/σε

300

350

400

450

500

Figure 5. Optimized product m · s, averaged over 5 realizations, as a function of inverse angle standard deviation, for 5 different autocorrelation lengths of the NURD. m · s is the inverse step size of the optimization process. The graph confirms the expected positive relation between step size and NURD magnitude. 100 σε/2π = 0.0025 90

σ /2π = 0.005

80

σε/2π = 0.0075

70

σ /2π = 0.015

ε

σ /2π = 0.01 ε ε

〈m⋅s〉

60 50 40 30 20 10 0

0

200

400

600

800 Γ

1000

1200

1400

1600

ε

Figure 6. Optimized product m · s, averaged over 5 realizations, plotted as a function of correlation length. A larger correlation length of the NURD is seen to require smaller steps.

3.1.1 NURD characteristics For every combination of σ and Γ , 5 realizations are simulated. The effective step size 1/m · s is expected to vary with the maximum frame-to-frame deviation, requiring larger steps (smaller m · s) for larger σ . The data in Fig. 5 confirm this, although there is a considerable amount of scatter. Similarly, Fig. 6 shows that a larger correlation length requires large m · s, or small steps. 3.1.2 Influence of speckle The NURD correction algorithm was tested on images with simulated speckle (the NURD parameters were σ /2π = 0.05; Γ = 1000). Correction of the NURD was still possible, as is demonstrated in Fig. 7. The rms difference between simulated NURD and the retrieved correction is two times larger than without speckle. The extra noise in the trace can be assessed by looking at the difference curve shown at the top of the graph. The panel on the bottom right shows that the speckle (individual frames look like Fig. 1(c)) can be averaged out quite well after correction.


8

6

angular shift (# lines )

4

2

0

2

4

6

8

0

1000

2000

3000

4000 5000 line number

6000

7000

8000

9000

Figure 7. Simulated frame to frame shift in a test sequence with speckle, compared to the optimal shifts as found by the algorithm (s = 8; m = 6). The images on the right show the mean of all frames before (top) and after (bottom) correction, demonstrating that the speckle can now be averaged out.

(a)

(b)

Figure 8. (Color in pdf) (a) RGB image of three subsequent frames (converted to cartesian coordinates) from an OCT movie in a vessel phantom with a plaque. Each frame was mapped to a separate color channel. An apparent counterclockwise rotation from red to blue can be seen. (b) The same frames after NURD correction; colors have merged to white because the frames are aligned. The bar corresponds to 0.5 mm.

3.2 In vitro imaging The effect of our NURD correction algorithm is most effectively visualized with actual movies. Movie clips demonstrating the method, both in vitro and in vivo, may be found on our web site, http://www.erasmusmc. nl/ThoraxcenterBME/research/additional/nurdcorr.htm. Here, we present a quantitative analysis of the performance of the technique. The effect of the NURD correction algorithm is demonstrated in Fig. 8. There is a larger apparent difference on the right hand side of the plaque than on the left between frames in the red, green and blue channels in Fig. 8(a). As a result, a global correction method would not have worked in this case. In Fig. 8(b), our new algorithm is seen to align the frames very well. The data was a stationary movie of a vessel phantom with an eccentric plaque. The performance of the NURD correction algorithm as measured by the mean normalized standard deviation over the aligned frame σ(a) in vitro was 0.285–0.295 for different sequences. For uncorrected sequences, this number was 0.44–0.47.

3.3 In vivo imaging We applied the correction algorithm to a stationary frame sequence, acquired in vivo. In Fig. 9 the alignment of three frames with strong NURD is shown. It shows the capability of the method to stabilize the image.


(a)

(b)

Figure 9. (Color in pdf) (a) RGB image of three subsequent frames from an intracoronary OCT movie in vivo, each mapped to a separate color channel. (b) The same frames after NURD correction. The bar corresponds to 0.5 mm.

Upon visual inspection of the corrected movie, it becomes apparent that the longitudinal motion of the catheter — arising because of cardiac motion — did not prevent the algorithm from correctly aligning the sequence. Features, such as the small side branch seen at 7 o’clock in Fig. 9(b), moved in and out of view with the cardiac cycle without changing their angular position. Neither did a small in-plane translation of the catheter that could be seen reduce the quality of the alignment.

4. CONCLUSIONS We have presented a new scheme to correct the effects of variations of the rotation speed in OCT image sequences, known as NURD. It optimizes the similarity between lines in subsequent frames to find a continuously varying correction shift to align the frames. The optimal sampling of the cost function, that represents the dynamics of the system most efficiently, was investigated as a function of the parameters characterizing the NURD. The technique is unique in its strictly continuous global optimization. Our NURD correction performs better than conventional GRBM schemes. We demonstrated the utility of the scheme for in vivo imaging.

REFERENCES 1. G. J. Tearney, S. A. Boppart, B. E. Bouma, M. E. Brezinski, N. J. Weissman, J. F. Southern, and J. G. Fujimoto, “Scanning single-mode fiber optic catheter-endoscope for optical coherence tomography,” Opt. Lett. 21(7), pp. 543–545, 1996. 2. J. G. Fujimoto, S. A. Boppart, G. J. Tearney, B. E. Bouma, C. Pitris, and M. E. Brezinski, “High resolution in vivo intra-arterial imaging with optical coherence tomography,” Heart 82(2), pp. 128–133, 1999. 3. I.-K. Jang, G. J. Tearney, B. D. MacNeill, M. Takano, F. Moselewski, N. Iftimia, M. Shishkov, S. Houser, H. T. Aretz, E. F. Halpern, and B. E. Bouma, “In vivo characterization of coronary atherosclerotic plaque by use of optical coherence tomography,” Circulation 111(12), pp. 1551–1555, 2005. 4. J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Optics Express 3(6), pp. 199–211, 1998. Article 123LU English Times Cited:7 Cited References Count:13. 5. F. J. van der Meer, D. J. Faber, D. M. Baraznji 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), pp. 1369–76, 2005. 6. B. E. Bouma, G. J. Tearney, H. Yabushita, M. Shishkov, C. R. Kauffman, D. D. Gauthier, B. D. MacNeill, S. L. Houser, H. T. Aretz, E. F. Halpern, and I.-K. Jang, “Evaluation of intracoronary stenting by intravascular optical coherence tomography,” Heart 89(3), pp. 317–320, 2003. 7. C. Schulz, R. A. Herrmann, C. Beilharz, J. Pasquantonio, and E. Alt, “Coronary stent symmetry and vascular injury determine experimental restenosis,” Heart 83(4), pp. 462–467, 2000. 8. G. van Soest, J. G. Bosch, and A. F. W. van der Steen, “Azimuthal registration of image sequences affected by non-uniform rotation distortion,” IEEE Trans. Inform. Tech. Biomed. (accepted; doi 10.1109/TITB.2007.908000), 2007. 9. D. L. Marks, T. S. Ralston, and S. A. Boppart, “Speckle reduction by I-divergence regularization in optical coherence tomography,” J. Opt. Soc. Am. A 22(11), pp. 2366–2371, 2005.


10. D. J. Berndt and J. Clifford, “Using dynamic time warping to find patterns in time series,” in AAAI-94 Workshop on Knowledge Discovery in Databases, pp. 359–370, 1994. 11. M. Sonka, V. Hlavac, and R. Boyle, Image processing, analysis, and machine vision, PWS Publishing, Pacific Grove, CA, 2nd ed. ed., 1999.
