Multi-shaping technique reduces sidelobe magnitude in optical coherence tomography

: Shaping methods that are commonly used in Fourier-domain optical coherence tomography (FD-OCT) can suppress sidelobe artifacts in the axial direction, but they typically broaden the mainlobe of the point spread function (PSF) and reduce the axial resolution. To improve OCT image quality without this tradeoff, we have developed a multi-shaping technique that reduces the axial sidelobe magnitude dramatically and achieves better axial resolution than conventional shaping methods. This technique is robust and compatible in various FD-OCT imaging systems. Testing of multi-shaping in three experimental settings shows that it reduced the axial sidelobe contribution by more than 8 dB and improved the contrast to noise by at least 30% and up to three-fold. Multi-shaping enables accurate image analysis and is potentially useful in many OCT applications.


Introduction
Optical coherence tomography (OCT) is one of the most successful optical imaging techniques to have been translated from bench to bedside in the past two decades [1]. It has become a standard imaging tool to examine the retina in ophthalmology, providing valuable diagnostic information [2,3]. In addition to retinal imaging, OCT plays important roles in high resolution imaging of both human and animal tissues [4,5]. Because OCT offers both structural and functional information with no reliance on exogenous labels or dyes [6], it has become a powerful label-free 3D imaging technique in both clinical and research settings.
A major artifact in OCT B-scan images is the sidelobe artifacts in the axial (z) direction. These sidelobes are a natural byproduct of the non-Gaussian spectrum emitted by the broadband light sources employed in OCT. To understand the origin of the sidelobes, consider that the axial point spread function (PSF), also known as the coherence function, is proportional to the Fourier transform of the light source spectrum [7] is the free-space wave number, z  is the path difference between reference and sample arms, and   Sk is the light source spectrum. When the spectrum of light source or the spectrum captured by the detector is non-Gaussian in shape, undesirable high sidelobes will be present in B-scan images. High sidelobes lower the image contrast and confound small target signals. To suppress axial sidelobes, a number of methods have been used to shape the non-Gaussian spectra, including Gaussian fitting [7][8][9], Hamming/Hanning windowing functions [10], and correction of spectra shape directly on the light source [11,12] or in post-processing [13][14][15][16][17]. These methods suppress sidelobes, but at a cost: they typically sacrifice some axial resolution by widening the mainlobe in the axial PSF.
In synthetic aperture radar (SAR) and medical ultrasound imaging, image quality suffers from the sidelobe contributions in lateral (x-y) directions, where they arise because the Fourier transforms of the rectangular aperture result in a sinc function in the spatial domain, resulting in sidelobes as high as 13.5 dB below the peak [18]. Sidelobe magnitude is commonly suppressed by applying an amplitude weighting function on the rectangular aperture data, including Hanning, Hamming or Chebyshev functions [18][19][20]. Such "apodization" processes work at the expense of worsened resolution (a widened mainlobe). To minimize the broadening of the mainlobe, advanced techniques, such as dual-apodization and tri-apodization by using more than one weighting function [18,20,21], have been deployed.
To minimize sidelobe contribution and maximize mainlobe resolution in OCT, we have developed a novel multi-shaping technique, using multiple windowing functions to shape the interference spectrum in FD-OCT. Unlike traditional shaping methods, the multi-shaping technique reduces the axial sidelobe magnitudes significantly without a major loss in axial resolution. This technique can be implemented via different forms, offering flexibility in its application in various scenarios.

Principle of multi-shaping technique
Multi-shaping is based on a simple principle: the position of the mainlobe in the axial PSF is largely insensitive to the width or position of windowing function used to shape the interference spectra in FD-OCT; in contrast, the positions of the sidelobes are very sensitive to the details of the windowing function. This difference in position with different windowing offers a simple means to recognize and suppress the sidelobes. While detailed model fitting might offer optimal suppression, such multiple-constraint deconvolutions can be time consuming [22][23][24][25][26]. Here, we show that even simple operations, such as a rank order (median or minimum) filter, is sufficient to eliminate the sidelobe contribution without the significant broadening of the mainlobe that results from other techniques.
The following outlines the procedures of this technique using a spectrometer-based spectral-domain OCT (SD-OCT) as the example: 1) A line scan detector (1 × N pixels) in the spectrometer is used to detect B-scans of interference spectra (raw B-scans) that contain structural information of the sample.
2) Generate a cluster of rectangular shaping functions centered at pixel N/2 with different rectangular window sizes of 1/2N, 1/2N + N step , …, N-N step and N pixels. 3) Multiply each raw B-scan with these shaping functions to create a cluster of B-scans of shaped spectra (shaped B-scan cluster).

4)
Process each member of the B-scan cluster by standard SD-OCT processing procedures, including background removal [27], spectral interpolation, numerical dispersion compensation [28], zero padding [27] and Fourier transforms.

5)
For each B-scan position, generate a cluster of B-scan images (B-scan image cluster), normalizing each member of the B-scan image cluster to its own maximum intensity.
6) Within each B-scan image cluster, process the intensity values for each corresponding pixel by using a mathematical operation, such as a rank order filter for the linear intensity. Use the result as the pixel value in the processed B-scan image.

7)
Smooth each processed B-scan image by a filter, such as a median filter or a Gaussian filter.
8) Display smoothed B-scan images as improved B-scan images.
In this paper, we used simple rectangular functions, instead of Gaussian or other functions, since rectangular functions do not broaden the mainlobe as severely as Gaussian functions and therefore maintain better axial resolution. We note that with proper deconvolution, many different windowing functions could be deployed, but typically with more complicated modeling and greater computational times. The rectangular windowing functions can be expressed as: where p is the pixel position, N is the total number of pixels in the line scan detector, and L is the width of the rectangular window. Four example rectangular shaping functions for a 4096pixel spectrometer (N = 4096) are shown in Fig. 1(a). Widths of rectangular windows are 4096, 3500, 3000 and 2500 pixels, respectively. Figure 1(b) depicts the cluster of B-scan images that result from using n different shaping functions on a single raw B-scan. An improved B-scan image is then generated by processing the intensity values for each corresponding pixel within the B-scan image cluster.

Numerical simulations of multi-shaping technique
Two numerical simulations were performed to illustrate and validate the multi-shaping technique. We obtained the greatest sidelobe reduction and best mainlobe resolution by using a simple minimum rank order filter to process intensity values (step 6) for both simulations. In the first simulation, an interference spectrum is simplified to have a constant frequency and fringe magnitude over all 4096 pixels of a spectrometer, as shown in Fig. 2(a). It is defined as: where   Sp is the intensity of spectrum at pixel p , and f is assumed to be 1/40. It corresponds to a single-layer mirror reflection, or axial PSF, in the spatial domain. A common shaping approach is to multiply the spectrum by a Gaussian function to obtain a Gaussianshaped spectrum, as shown in Fig. 2(b). In the multi-shaping technique, PSFs generated from four rectangular shaping functions in Fig. 1(a) and combined into a single PSF by a minimum rank order filter are shown in Fig. 2(c), illustrating the basic operation of this technique. After Fourier transform, the original spectrum ( Fig. 2(a)) results in a PSF with a narrow mainlobe but high sidelobe levels (13.4 dB), shown as the yellow curve in Fig. 2(d). For the Gaussian-shaped spectrum ( Fig. 2(b)), the resulting PSF offers reduced sidelobe levels (34.2 dB) but the mainlobe width (at 6 dB) is 53% broader, shown as the green curve in Fig. 2

(d).
In the multi-shaping example of Fig. 2(d), the processing used a total of 42 rectangular functions, each centered at pixel 2048, the rectangular windows have widths from 2000 to 4050 in 50-pixel increments, in addition to the full width of 4096 pixels. The improved PSF, shown as the red curve in Fig. 2(d), significantly reduces the sidelobe level to 41.1 dB, which is 6.9 dB below the Gaussian-shaping function described earlier. What is more, the mainlobe width is as small as the original PSF, which indicates the best resolution. In a more demanding simulation, we considered the performance of the multi-shaping technique in identifying two nearby reflections in space. The simulated interference spectrum is defined as: A is as small as 1% of 1 A , which is almost 40 dB lower. These simulation results indicate that the multi-shaping technique is able to identify features that are very close in space, even when one feature is dramatically dimmer. A is much smaller than 1 A .

OCT imaging experiments with the multi-shaping technique
The multi-shaping technique was applied to experimental OCT and OCT microscopy (OCM) data acquired from three specimens: a human eye, a microbead phantom, and a zebrafish embryo. OCT imaging data of the human retina were acquired from a custom SD-OCT system (2048-pixel spectrometer) developed at Varocto Inc, as described in our previous work [29]. OCM imaging data of the microbead phantom and zebrafish embryo were acquired from custom spectral-domain OCM (SD-OCM) systems (4096-pixel spectrometer) described in our previous studies [30,31].
Due to the finite spectral resolutions and pixel sizes of the spectrometers, the fringe magnitude of the detected interference spectrum was uneven. Therefore, in contrast to the numerical simulations in which all shaping functions centered on the central pixel, we used a larger family of shaping functions, shifting the center of rectangular windows for processing the experimental data. Table 1 summarizes the 67 rectangular shaping functions used for OCT data sets (N = 2048), created with center positions ranging from pixel 500 to 1500, with 100pixel increments, and window widths of 100-pixel increments from 1000 to 2000. We included a full window width of 2048 pixels and removed those that would not fit within pixels 0~2047. Similarly, we created 232 rectangular shaping functions for the OCM data sets (N = 4096), with center positions ranging from pixel 1000 to 3000 at 100-pixel increments, and window widths of 100-pixel increments from 2000 to 4000 and including a 4096-pixel window width. In all data sets, we chose rectangular window widths between around N/2 to N pixels and increments of 100 pixels, because we found that window widths below N/2 and increments smaller than 100 pixels did not noticeably reduce the sidelobe noise. As with the previous numerical simulations, we generated intensity values in improved B-scan images with a minimum filter, which appeared to offer the greatest reduction of sidelobes and the best mainlobe resolution. All B-scan images were then smoothed by a 3 × 3 median filter and displayed in a 50-dB dynamic range.
For our experimental results of imaging the human retina and zebrafish embryo, the image quality was evaluated quantitatively by comparing the contrast-to-noise ratio (CNR) between sample and background regions, given by the equation below:  represent the mean and standard deviation of the background noise region [25,32,33]. These parameters were calculated from logarithmic OCT/OCM images, in which the pixel values were in dB scales within the 50-dB dynamic range. In our quantitative analysis, we used three regions (N = 3) to calculate CNR values.

Imaging of microbead phantom
The multi-shaping technique was applied on experimental data obtained from imaging a microbead phantom by our extended-focus OCM (xf-OCM) system [31]. The phantom, used to characterize the performance of the xf-OCM system, consisted of latex microbeads of 0.3μm diameter (LB-3, Sigma-Aldrich) immersed in a 0.75% phytagel solution, with a concentration of ~10 8 particles/mL. Such a concentration can provide sufficient numbers of microbeads in OCM images, while reducing their tendency to form aggregates of several beads. Microbeads with 0.3-μm diameter are sufficiently small to permit measuring the system's PSF (lateral resolution = 1.6 μm, axial resolution = 2.9 μm), while large enough to generate strong reflection signals for analysis. B-scan images of the microbead phantom with and without the multi-shaping technique are shown in Fig. 4. The comparison of Fig. 4(a) and 4(b) reveals that the multi-shaping technique reduced undesired background noise and provided better visualization of individual beads. The profile of single-bead reflections was the same as the measured axial PSF of the xf-OCM system (Fig. 4(c)). The plot shows that the sidelobe magnitude on the right side of the mainlobe of a single bead was reduced from 17.7 dB to 25.8 dB, compared with standard Gaussian shaping. There is a smaller peak on the left side of the mainlobe, which possibly came from the reflection of another bead in the neighboring area. It only decreased from 16.3 dB to 21.9 dB, suggesting that the multi-shaping technique preserved small signals that were next to bright signals, as predicted by our second numerical simulation.

In vivo imaging of zebrafish embryo
We applied the multi-shaping technique to in vivo imaging data from a 3 days-postfertilization (dpf) zebrafish embryo. The protocol for the animal experiments, such as mounting the specimen in phytagel, was described in our previous publication [30]. Images of the embryo were acquired using our SD-OCM system, with 2.8-µm lateral and axial resolution [30]. We elected to image the fish embryo in the mid-trunk region because of its variety of complex cellular features, including thin cell boundaries (~5-µm thick) inside the notochord, which are appropriate for evaluating the performance of our technique. Figure 5 shows that after applying the multi-shaping technique to identical image data sets, the background within the notochord became much cleaner, resulting in much better identification of cell boundaries inside. In the middle of the notochord, the features where two cell layers meet are much more clearly resolved in Fig. 5(b) due to reduced sidelobe intensity. The CNR was increased from 9.0 (Fig. 5a) to 23.1 (Fig. 5(b)), more than a two-fold improvement. Axial intensity profiles in Fig. 5(c) show that reflection signals from two cell layers became better identified thanks to the multi-shaping technique, as the gap between them was decreased by 5 dB, while their intensities remained roughly the same. The background inside the notochord was nearly 15 dB lower.

In vivo imaging of human retina
B-scan images of the human retina from a healthy subject (male, 38 years old), shown in Fig.  6, were acquired by a custom-built SD-OCT system developed at Varocto Inc. This system was similar to that was used in our previous OCT work [29], with a center wavelength of 835 nm and a bandwidth of 50 nm, offering 4.5-µm axial resolution in the retina (n = 1.38). The OCT scan was performed over a retinal area of 3 × 3 mm 2 .
As shown in Fig. 6(b), the B-scan image processed by the multi-shaping technique had much lower noise in both sample and background regions. The CNR of Fig. 6(b) was improved by almost 30% compared with Fig. 6(a), increasing from 5.3 to 6.8. The improved image was superior in resolving not only bright layers, such as the junction between the inner and outer photoreceptors (IS/OS), retinal pigment epithelium (RPE) and choriocapillaris (CC); but also features within relatively dim layers, such as the outer plexiform layer (OPL) and external limiting membrane (ELM). Moreover, weak signals of the choroid (C) region can be better identified, as they are not overwhelmed by the background noise in Fig. 6(a). Axial intensity profiles in Fig. 6(c) show that reflection signals from bright layers (IS/OS, RPE and CC) became sharper, and that the intensities of the hypointense regions between these layers were lowered down by more than 10 dB.

Discussion
This paper presents a simple and straightforward technique that can suppress axial sidelobe magnitude without sacrificing the axial resolution in FD-OCT or adding significant complexity to the standard FD-OCT processing procedures. Our Multi-shaping approach is not only easier to implement than the alternative approaches that involve Gaussian fit, adaptive apodization or search of optimal shaping functions [7][8][9][10][11][12][13][14][15][16][17], but also avoids the usual tradeoffs of most techniques in which good mainlobe resolution is sacrificed for good sidelobe reduction. Most shaping techniques require precise measurement of the light source spectrum; in contrast, the performance of multi-shaping is independent of such measurements. To better define these differences, we compared the multi-shaping technique with an automatic spectral reshaping method using a Wiener filter [16,17] for the zebrafish data (presented in Fig. 7) and found that its performance offered results similar to the standard Gaussian shaping (CNR = 9.0, Fig. 5(a)), with roughly the same levels of contrast (CNR = 9.2) and sidelobe artifacts. Image quality was not improved, mainly because the broadband superluminescent diode (SLD) used in this OCM system had strong spectral modulations and deviated from the Gaussian shape [30]. In contrast, multi-shaping achieved improved results ( Fig. 5(b)-5(c) and Fig. 7(b)). Multi-shaping is computationally simple, requiring the processing of multiple rectangular windows applied to the spectra to create a B-scan image cluster for each B-scan position. On our computer (Python 3.6, Intel Core i7, 2.40GHz, 16GB RAM, 64-bit Windows 10), the time required for processing a standard OCT B-scan with standard Gaussian shaping (Fig. 6) was 0.5 sec; with multi-shaping is was 5.4 sec. This could be further sped by parallel computing using a GPU or by multi-core processing, since the OCT processing by different shaping functions are independent processes.
In both the simulations and the experimental results shown here, we used a minimum rank order filter, as it was simple to implement and achieved good mainlobe resolution and strong sidelobe reduction; however, this might result in as intensity loss, raising concerns that weakly scattering structures might become even dimmer. For example, in Fig. 6(b), some parts of the structure around the ONL layer almost disappear compared with (a). If this proves problematic in some settings, less aggressive multi-shaping can be performed using a higher rank order filter, such as a median filter or a rank order filter in between the median and minimum filter. A higher rank order filter should increase the intensity of the processed image and better preserve signals of the weakly scattering structure, but with slightly worse mainlobe resolution and weaker sidelobe reduction. To illustrate different options of filters, we processed the B-scan of the human retina by multi-shaping with a median filter ( Fig. 8(a)) and a rank order filter ( Fig. 8(b)) picking the 17 th smallest value over 67 B-scans within a Bscan image cluster. The CNR of Fig. 8(a) and 8(b) was 7.5 and 6.4, both were better than Fig.  6(a). Axial intensity profiles in Fig. 8(c) show that although reflection signals resulted from the median filter and rank order (17 th smallest) filter were not as sharp as the minimum filter, they preserved weak signals similarly to Gaussian shaping while achieving axial resolution better than it and slightly worse than the minimum filter. To further evaluate the multi-shaping technique, we performed a histogram analysis of the standard and improved images in Figs. 6 and 8, using the regions selected for the CNR calculation (yellow dashed boxes, Fig. 6). As shown in Fig. 9, multi-shaping by the minimum filter reduced the averaged intensity (µ) by 8 to 10 dB relative to standard Gaussian shaping; in contrast, multi-shaping by the median filter and rank order (17 th smallest) filter are more moderate approaches, with only 0 to 2 dB increase or decrease. In addition, multi-shaping by the minimum filter increased the standard deviation (σ) of pixel values in the sample region, possibly because speckle patterns were highlighted by the minimum filter rather than being masked by the sidelobe contribution. In the background region processed by the minimum filter, majority of the pixel values centered around zero because they were very low and below the 50-dB dynamic range for visualization.  Fig. 6. From top to bottom: results processed by standard Gaussian shaping, multishaping with a minimum filter, multi-shaping with a median filter, and multi-shaping with a rank order (17 th smallest) filter. The region size is 50 × 50 pixels. The averaged intensity (µ) and the standard deviation (σ) of pixel values in these regions are presented.
Multi-shaping requires the selection of appropriate window widths and number of the rectangular shaping functions to achieve a proper balance between sidelobe reduction, intensity loss and processing time. Using more shaping functions, we are likely to achieve stronger sidelobe reduction and higher contrast, but at the price of greater intensity loss and longer processing time. In this work, we selected the widths and number empirically, based on the observation that window widths below half of the amount of spectrometer pixels and increments smaller than 100 pixels did not lead to noticeably higher sidelobe reduction. Moreover, the processing time with current parameters were reasonable for reconstructing a volume of imaging data. In our future work, we will develop a quantitative method of determining these parameters to optimize the performance of the multi-shaping technique. The method should be applicable to various OCT imaging systems of different light sources, axial PSFs and noise levels.

Conclusion
We have demonstrated a novel multi-shaping technique for reducing sidelobe artifacts in FD-OCT imaging systems. The multi-shaping technique is able to improve the image quality remarkably without the price of degrading the axial resolution. Our experimental results show that the axial sidelobe contribution was reduced by more than 8 dB and the CNR was enhanced by at least 30% to three-fold. It can be a very efficient method because of its compatibility with parallel computing. It is potentially useful for various applications, as it captures sample features in greater detail and permits more precise layer segmentation in retinal imaging. In this simple implementation, the multi-shaping technique is a nonparametric method because it does not rely on noise models to remove sidelobe artifacts in images, permitting it to be robust and applicable in a wide range of FD-OCT imaging systems.

Funding
This work was funded by the Translational Imaging Center, the Alfred E. Mann Institute for Biomedical Engineering, and the Viterbi Doctoral Fellowship at the University of Southern California.
Disclosures Y.C. and S.E.F filed a provisional patent application on the multi-shaping technology.