O 2003 The Optical Society of Japan
ProceSS Algorlthms for Resolutlon Improvement and Contrast Enhancement
in Optical Coherence Tomography
I-Jen HSU, Chia-Wei SUN, Chih-Wei LU, Chih-Chung YANG*, Chun-Ping CHIANGl, Chii-Wann LIN2
Graduate Institute of Electro-Optical Engineering and Depa"tnrent of Electrical Engineering, National Taiwan University, J, Roosevelt Road, Section 4, Taipei, Taiwan, R.O. C.
IDepartlnent ofDentistry, National Taiwan University, 1, Roosevelt Road, Section 4, Taipei, Taiwan, R. O. C.
2Graduate Institute ofBioengineering, National Taivvan University, J, Roosevelt Road, Section 4, Taipei, Taiwan, R. O. C. (Received March 6, 2003; Accepted october 9, 2003)
We proposed two process algorithms for resolution improvement and contrast enhancement in the images of optical coherence tomography (OCT). An OCT system with a non-Gaussian light source spectrum or dispersion mismatch
usually results in sidelobes in the interference fringe envelope that may produce artifacts and reduce image contrast.
Based on the concept of deconvolution, we proposed two different process algorithms and demonstrated their effectiveness in retrieving sample structures. The effects of the process algorithms were examined by numerical simulations and real OCT scanning images. After processing with the proposed procedures, the effects of sidelobes were tremendously suppressed and the image qualities were improved.
Key words: optical coherence tomography, process algorithm, image resolution, image contrast, deconvolution
1. Introduction
Optical coherence tomography (OCT) is a noninvasive
imaging technique and has been widely studied for biomed-ical applications.1) Because the longitudinal resolution of an OCT system is inversely proportional to the spectral width of its light source, broadband light sources are required for
high-resolution OCT systems. Current developments for
improv-ing the spatial resolution of OCT systems by use of ultra-broadband light sources include the use of a photonic crystal
fiber in combination with a compact sub-10fs Ti:sapphire laser to generate a spectrum ranging from 550 to 950nm.
This broad spectrum results in submicrometer axial
resolu-tion in an OCT system 2) Also, the use of a spectral
continuum from 800 to 1400 nm generated from a photoniccrystal fiber to produce longitudinal resolution of I .3-um in
biological tissue was reported.3) Meanwhile, a continuum
generation through an ultrahigh-numerical-aperture fiber 4) and a superluminescent light source using a Ti:A1203 crystal
pumped by a frequency doubled diode-pumped laser5) were
used for high-resolution OCT systems.
Although an ultra-broadband light source can result in ultra-high resolution in an OCT system, it usually has a non-Gaussian spectrum and produces sidelobes in the interference fringe envelope. The sidelobes can be suppressed by shaping the spectrum into a Gaussian form. However, such a process will reduce the spectral width and degrade the longitudinal
resolution of the OCT system. If appropriate process
algorithms can be applied to suppress the sidelobes in the interference fringe envelope, the ultra-broadband light source
can be used without spectral shaping and the ultra-high resolution of the OCT system can be preserved. Several
signal-processing methods have been developed to increase the resolution of OCT images, such as deconvolution with the iterative restoration algorithms6~8) and reduction of the
side-9)
10bes with a spectral shaping technique.
In this paper, we propose and demonstrate two different
process algorithms to improve the resolution and to enhance the contrast of OCT images. The effects of these algorithms are demonstrated with numerical simulations and shown with real OCT scanning images. The applicability and limitation of the process algorithms are also discussed. In section 2 of
this paper, the concepts of the proposed algorithms are
discussed, and their procedures are described in section 3. The results of numerical simulations are presented in section 4. Application of the proposed algorithms to practical OCT scanning results is shown in section 5. Finally, conclusions and discussion are given in section 6.
2. Concepts of the process Algorithms
An OCT image iS a collectiOn of A-scan signals that are Obtained from the diStributiOns of interference patterns of
PhOtons backScattered by the interfaces of the sample interfering with the phOtons from the reference arm. The
detected signal in an A-scan iS
::: I ) 2 * 12)
(1( I l
I(r) ~E (t+ r) +E (t+ T)
(1)
:= 11 + 12 + 2Re ~~(E1(t + r )E (t + T))
Here, El and E2 stand fOr the wave fieldS from the sample and reference arms, respectively, and 11 and 12 are the intensities in the two arms. They are defined as
l =: I I 12)
(1 1 ,
~E (t + T)
12 =: ( E2(t) 2). (3)
*E-mail address: ccy@cc,ee.ntu.edu,tw
The time delay T is measured between the reflection at the mirrOr in the reference arm and that at the first interface of the sample. ri denotes the time delay of the photons from the ith
interface of the sample relative to those from the first
interface. The first two terms in Eq. ( 1) correspond to dc
components during scan and will be neglected in the
following discussion. For a light source with the spectral function S((D), the wave field can be expressed as the inverse-Fourier transform of the spectral function as
1 c<)
= foo
E(t) 27T S((~)) exp[jip((~))] exp(ftDt)da). (4)
Here, ip(a)) represents the frequency-dependent phase. The
one-dimensional OCT scanning signal obtained from a
single-facet reflection in the sample arm is related to the autocorrelation function oo
f
P(T) = E1(t)E~(t + T)dt
oO 1 IS(a)) 1 2 exp(-fa)T)da)= f
2 JT -OQWe assume that the dispersion between the sample and
reference arms is balanced and the spectrum is not modified in either arm. For a sample with multiple interfaces, the interference term of the detected signal will be
oo
1 i IS(a)) 12 exp[-fa)(T - Ti)]da), (6)
f
I(T) = - Re ~A
Jc ,
l -oo
where Ai is the amplitude fraction of wave field scattered
from the ith interface of the sample. If we consider a
Gaussian spectrum with the spectral density
- )J
[( 2,
IS(a))[2 = S20 exp a) - (~)o (7)
A(L)
the detected signal in Eq. (6) becomes
S2A(~) Aa)2(T - Ti)2
I(T) = fi ~Ai exp 4 (8)
x cos[a)o(T ~ 7ri)].
A one-dimensional OCT image is obtained from the
envelope of the detected signal after demodulation of the carrier frequency. The detected signal can be regarded as a superposition of interference fringe patterns with different amplitudes at different depths. When the distances between the two neighboring interfaces of the sample are larger than the widths of the individual interference fringe envelopes, the
demodulated signal is approximately a convolution of the
function describing the sample structure and the normalized interference fringe envelope of a single-facet reflection. Under this approximation, the demodulated detected signal can be expressed as
N
• ~
Il = !jfj_j' (9)
i= lHere, i and j denote the discrete coordinates that represent
the depths in the sample, and N is the pixel number. In
comparison with the undemodulated signal described in Eq. (8), I/ corresponds to the coefficient including Ai, and represents the intensity that describes the sample structure. The function f corresponds to the exponential factor in Eq.
(8) and stands for the normalized interference fringe
envelope of a single-facet reflection. fo rs defined as the peak and has the value I . The first step of our process
algorithm is to find a certain function as an assumed sample structure to be the initial guess in the procedure of iterations.
The function of the assumed sample structure is then
modified to approach the real one by comparing the
demodulated detected signal to the result of convolving the
assumed function with the normalized interference fringe
envelope of a single-facet refiection. Note that when the distance between two neighboring interfaces in the sample is less than the width of the individual interference fringe envelope, the two interference fringe patterns overlap. In this situation, the demodulated signal cannot be simply consid-ered as a convolution of the sample structure function and the normalized interference fringe envelope of a single-facet refiection. Therefore, the proposed algorithms can be used only when the separation between two neighboring structure features in a sample is not much smaller than the resolution
cell size. The applicability and the limitation of the
algorithms will be discussed in the next section.
3. Procedures of the Process Algorithms
We propose two different methods for constructing the assumed function as the initial guess of iterations. The
flowchart of the first method is shown in Fig. 1(a). With the definition in Eq. (9), the detected intensities at various pixels I - IN are constructed to form a vector I. The relevant parts of the function f are constructed to form the vectors fl ~ fN. The inner product of fi and I, namely, Ji (see Fig. 1(a)), is exactly the ith component of the convolution of functions I and f. The initial function I' is obtained by dividing the detected intensity I by r, where r is the ratio of the maximum values of Jj and li. In the second method, the initial guess
function is constructed with the local maximum of the
detected intensity I. Its flowchart is shown in Fig. I (b). The value of the initial guess function at any pixel I'i is set equal to li When li is a local maximum and set to O at other pixels.
The constructed initial guess function is then processed with the procedure described in Fig. 2. Here, I"i is the inner
product of fi and I', and is the ith component of the convolution of the initial guess function I' and f. The
difference between I" and the detected intensity I at each pixel, di, is then used to correct the initial guess function and
to check whether the iteration process converges to a
satisfactory solution. The summation of the absolute values of di Of all pixels can be regarded as a criterion of satisfactory solution. Before the criterion is reached, I'i is subtracted by cdi at each pixel to obtain a new function I"', where N is the pixel number and c is an adjustable parameter used to control the degree of correction of I'. In the operation, the choice of an appropriate c value is crucial. The iteration procedure may not converge to a satisfactory solution if c is too large. The result of a particular iteration is obtained by selecting the positive values of I"'i and setting others at O because negative intensity is unreasonable. This result is used as the guess function of the next iteration. The process is terminated when a satisfactory solution is achieved.
4. Numerical Simulations
limita-OPTICAL REVIEW Vol.
Input I I I I
l, 2' 3""' N'f_ f f_ ...,
N+1 ' J _N +2 ' N +3 ' fN_l I :: (1 1 1 IN ) 1' 2' ~""' f = (f. f. fl 3"~ fi N ) l ,1 , l~2 ' ~ ' -Ji = fi ' I + I; = Il /' Output I' Ir If ,1~ l, 7' 3"" (a) Input I I I I l' 2' ~""' N'If, I,r, I~,..., I~,
f N+1 ' f N+" f N+3-"' flv-l I'=(II I' Il I") I' " 3""' l\ f (f, I f fr 3' "' f_N ) I:'=f 'I' d =1!'-I
env' d
':= l
~1 i N i= l I:"= F - cd I:"> O I; = I;"' r oI=
r Input I I I , IN 1' 2' 3"" Output lr lr I' Il' l, 2' 3""' I > I and , i-l li > Ii+]i]
Output I;, I~, I~,•••, I~
I =0 ,
i
(b)
Fig. I . Flowcharts showing the procedures of constructing the initial guess functions with (a) the first method and (b) the second
method.
tions of the proposed algorithms, we use the iteration procedures for simulating the imaging operation of an
assumed sample structure with six interfaces, as shown in
Fig. 3(a). Here, the separations between two neighboring
interfaces range from I to 5 um. Such a sample structure is scanned with an OCT system with the spectral distribution of the light source described in Eq. (7). The center frequency a)o is assumed to be 2.355 x 1015 Hz, which corresponds to the
center wavelength of 800 nm. The full-width at half max-imum (FWHM) of the spectrum was chosen to be 100 nm by
setting Aco = 0.176 x 1015 Hz. Under these conditions, the FWHM of the interference fringe envelope of a single-facet reflection can be calculated through the equation
Az 2ln2 A~ (10)
JT AA
to obtain Az = 2.84 um. This value is conventionally
regarded as the longitudinal resolution of the OCT system. In simulation, the detected signal for the assumed sample
Fig. 2. Flowchart showing the retrieval iteration procedure.
structure can be calculated through Eq. (6) or (8) and is shown in Fig. 3(b). The corresponding one-dimensional OCT image is obtained by demodulating the detected signal and is shown in Fig. 3(c). One can see that the individual interfaces
can be well identified when the separation between two neighboring interfaces is larger than 3 um. When the
separation is smaller than 2 um, the interfering signals resulting from different interfaces overlap each other. Hence, the positions of the interfaces can no longer be correctly
identified.
The signals were processed with the proposed algorithms. Figure 3(d) shows the processed results with the first method after 10 (dotted curves) and 100 (solid curves) iterations, respectively. When the separation between two neighboring interfaces is larger than 3 um, the width of each envelope becomes smaller, corresponding to the contrast enhancement of the image. When the separation is smaller than 2 um, the original interfaces cannot be retrieved. Figure 3(e) shows the processed result with the second method after 10 iterations. In this situation, the retrieved result looks better when the interface separation is larger than 3 um. The interfaces with separations smaller than 2 um cannot be retrieved with the second method either. In both methods, the processed signals approach the convergent results after 10 iterations when the parameter c is set at 0.02. The limitation of the proposed
algorithms was estimated by considering a two-interface structure. From the raw data of OCT image, we cannot
distinguish the two interfaces or correctly identify their positions when the separation is less than 3 um. Figures 4(a) and (b) show the simulated signal and its envelope of the two-interface structure with 2 um in spacing. The interference fringe patterns from the two interfaces overlap each other and
1~]-l O-~ OO-~- O~-~Q - Oi5-~~ (!~ c 04-(D c: - O~}-O J) O IS Ifl- O~-:~ (Q - OJ]-;~ (J1) L: *OS-(D 2: - -1J~--1 5 5 lo i5 Depth (um) 20 25 15 I~-:~ OS ~Q - OJ]-~~ (!) c ~Oj;-(D ~ - -1 O--1 5 O 1 5-O 5
10 '
15 Depth (hLm) 20 2s~ ~ I b I ~ i ~ 2b 2~
~ ~ l'4 l~ Depth [~m) :i IJ]-(Q ~~ ~) ~: O; G) 4~ ~: l 5-(b) :~ IJD-CQ ~~ (!) c OS (L) c: O~ o 5 10 Is Depth (~tm) 20 2so o lb I~ l~ 20 22
d ~ ~ ~ ~ 1'4 1'6 Depth (um) O 03 -:~ 0.02-(Q~
cr) OOl c (1) c (c)o oo rb l~ 1~ 2~
d ~ ~ ~ ~ 1'4 1'5 2b Depth (urr]) O Os O J:)4 -~ OJD3 ~O ~~vl OJ]2-c (1> ~i 0~il O J~0 (d) O s lo 15 Depth (~Lm) 20 2sFig. 4. (a) The simulated interfering signal and (b) its fringe envelope of a two-interface structure with 2 um spacing. (c) The processed results after 10 (dotted curves) and 100 (solid curves) iterations by the first method.
1 S-:~ I~-~Q ~~ (j' OS c: Q) c OO (e) O s 10 15 Depth (~rn) 20 25
Fig. 3. (a) Assumed sample structure with six interfaces. (b) The simulated detected signal from the assumed sample structure and (c) its fringe envelope. (d) The processed results with the first method after 10 (dotted curves) and 100 (solid curves) iterations. (e) The processed result with the second method after 10 iterations.
cannot be distinguished. The processed results after 10 and lOO iterations with the first method are shown in Fig. 4(c)
(dotted curve for 10 iterations and solid curve for 100
iterations). The two interfaces become distinguishable and their positions can be correctly identified after the retrieval process. This means that the longitudinal resolution can be improved with the proposed algorithms if the image feature separation is not much smaller than the resolution cell size. We have found that the two interfaces can be well resolved if the interface separation is larger than 0.8 um, although their positions might not be correctly identified.
5. Experimental Results
In experiments, we used a typical free-space OCT system with a mode-10cked Ti:sapphire laser (Femtosource) as its
light source to scan various samples. The light source
provides a non-Gaussian spectrum with the spectral width of
about 80nm. Figure 5(a) shows the one-dimensional OCT
image of a human buccal mucosa sample. Before the retrieval
procedure, the raw data were smoothed through a low-pass
filter to eliminate the high-frequency noise. Figure 5(b) shows the image of filtered data. This set of data was then processed by the first method. The processed image after 10 iterations is shown in Fig. 5(c). One can see that not only does
the front boundary become sharper, but also the sample
structure becomes clearer, and the image contrast is signifi-cantly enhanced.
Figure 6(a) shows a two-dimensional OCT image of the same sample. Figure 6(b) shows the same image after
low-pass filtering. The image after processing with 10 iterations of the first method is shown in Fig. 6(c). Here, we can easily distinguish the layered structures including the (i) keratin, (ii) epithelium, (iii) connective tissue and (iv) salivary glands, as indicated in the figure. The image contrast was significantly enhanced and the boundaries between the layers are clearer.
6. Conclusrons
In conclusion, we have proposed two process algorithms to retrieve the sample structure and enhance the image contrast
OPTICAL REVIEW Vol. 04 03-::; CQ_02 ~~ ~ co ~: ~~~ O1-c OO O 03
SO ISO
100 200 2so 300 3so 400 450 500Depth (um) =; 0~-CQ
~
(!) O1 C (D C: cO (b) O O1 O 50 1 OO 1 50 200 250 300 3so 400 4so 500 Depth (um) CQ ~~ oD C Q) C O OO (c)O SO IdO I~O 200 250 300 3so 400 4so soo
Depth (um)
Fig. 5. One-dimensional backward-scattered intensity distribution in the OCT image of a human buccal mucosa sample: (a) before the
low-pass filter, (b) after the low-pass filter, and (c) after 10 iterations
with the first method.
used the methods of iteration to approach the real structure of
the sample. The effects of the process algorithms were
demonstrated with numerical simulations and real OCT
scanning images. The application limitation of the algorithms
was also examined. In real OCT operation, a human buccal mucosa sample was scanned with a free-space OCT system. After processing with the proposed algorithm, the sample
structure became clearer and the image contrast was
significantly enhanced.
Acknowledgement
This research was supported by the National Health Research Institute of
(i) (ii) (iii) (iv) -max .{~; -m n
Fig. 6. Two-dimensional OCT images of the human buccal
mucosa sample (2 mm x 0.5mm): (a) before the low-pass filter,
(b) after the low-pass filter, and (c) after 10 iterations with the first
method.
The Republic of China, under the grant NHRI-EX92-9220EI.
Ref erences
1) D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson. W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito and J. G.
Fujimoto: Science 254 (1991) 1 178.
2) B.Povazay, K. Bizheva, A,Unterhuber, B.Hermann. H. Sattmann, A. F.
Fercher, W. Drexler, A. Apolonski, W. J. Wadsworth, J. C. Knight, P. St. J. Russel, M. Vetterlein and E. Scherzer: Opt. Lett. 27 (2002) 1800. 3) Y. Wang, Y. Zhao, J. S. Nelson, Z. Chen and R. S. Windeler: Opt. Lett.
28 (2003) 182.
4) D. L. Marks, A. L. Oldenburg, J. J. Reynolds and S. A. Boppart: Opt.
Lett. 27 (2002) 2010.
5) A. M. Kowalevicz, T. Ko. I. Hartl, J. G. Fujimoto, M. Pollnau and R. P. Salath~: Opt. Express 10 (2002) 349.
6) M. D. Kulkarni, C. W. Thomas and J. A. Izatt: Electron. Lett. 33 (1997) 1365.
7) M. Bashkansky, M. D. Duncan and J. Reintjes: App. Opt. 37 (1998) 8137.
8) I. J. Hsu. C. W. Sun, C. W. Lu, C. C. Yang, C. P. Chiang and C. W. Lin:
Appl. Opt. 42 (2003) 227.
9) R. Tripathi, N. Nassif. J. S. Nelson, B. H. Park and J. F. de Boer: Opt. Lett. 27 (2002) 406.
