Determination of confocal profile and curved focal plane for OCT mapping of the attenuation coefficient

: The attenuation coeﬃcient has proven to be a useful tool in numerous biological applications, but accurate calculation is dependent on the characterization of the confocal eﬀect. This study presents a method to precisely determine the confocal eﬀect and its focal plane within a sample by examining the ratio of two optical coherence tomography (OCT) images. The method can be employed to produce a single-value estimate, or a 2D map of the focal plane accounting for the curvature or tilt within the sample. Furthermore, this method is applicable to data obtained with both high numerical aperture (NA) and low-NA lenses, thereby furthering the applicability of the attenuation coeﬃcient to high-NA OCT data. We test and validate this method using standard samples of Intralipid 20% and 5%, improving the accuracy to 99% from 65% compared to the traditional method and preliminarily show applicability to real biological data of glioblastoma acquired in vivo in a murine model.


Introduction
Optical coherence tomography (OCT) is an imaging technology that has gained popularity especially in recent years due to its capability to produce high-resolution (on the order of micrometers), label-free, three-dimensional images in vivo [1]. OCT is an interferometric technique based on the back-scattering of near-infrared light, allowing various structures within an image to be visualized by differences in optical properties [2]. OCT thus easily lends itself as an effective tool to probe optical properties of samples, including the attenuation coefficient: a measure of the degree of attenuation of light within a sample. Several studies have demonstrated the usefulness of this measure in a host of various medical contexts, namely the quantification of ovarian collagen content [3], detection of cerebral edema [4], diagnosis of cancer [5,6], bladder biopsies [7], quantification of atheroscleretic plaque [8,9], glucose monitoring [10][11][12], etc. There is also a growing interest in tools for functional mapping during surgery such as the use of OCT intraoperatively, for example in the delineation of cancerous and healthy tissue in the brain during surgery [13,14]. The accuracy of resection during brain surgery is especially important in eloquent areas, such as in the speech and motor areas, where it is imperative to minimize surgery-inflicted deficits due to the removal of non-cancerous [15]. One recent study proposed to detect cancer infiltration by thresholding the attenuation coefficient, demonstrating the potential importance of high-resolution mapping of the absolute value of the attenuation coefficient [13]. In such applications, it becomes especially critical that the attenuation coefficient is calculated accurately and robustly. law, which is a single-scattering model that describes signal attenuation through exponential decay [16]. Although another model has been introduced to additionally encompass multiple scattering in highly scattering samples, generally the contribution of multiple scattering to OCT signal of biological tissue is fairly small [17,18]. However, an important component that needs to be considered when calculating the attenuation coefficient from OCT data is the confocal effect, which arises from focusing of Gaussian beams for transverse resolution. The modulation of OCT signal by confocal effect has been characterized in literature and depends on two factors: the position of focus within a sample, and the apparent Rayleigh range [19]. The position of focus is often simply determined by moving the sample by a known distance, assuming that the focal plane is flat; however, the robustness and accuracy of such methods are questionable as detailed below, especially when a 2D map of the attenuation coefficient is desired. Such mapping with high transverse resolution can lead to a particularly severe confocal effect, which has limited the use of the attenuation coefficient when using high numerical aperture (NA) lenses that have small Rayleigh ranges, thus making it critical to accurately specify the position of the focal plane [20]. Addressing this issue would extend applicability of the attenuation coefficient in clinical usage, where high magnification lenses are required to examine tissue microstructure.
Commonly, the focal plane position is found manually during imaging by moving the surface of the sample through the maximum field of brightness [21]. This method relies on assumptions of the behavior of the beam in the sample; when moving the lens (or the sample), the focal plane is shifted by an equal distance within the sample. This assumption, however, is not accurate when the refractive index of the sample is different from that of the air (the case in most biological samples). Regarding the confocal profile, existing studies have opted to ignore it completely [4,11,12] or otherwise employed methods that avoid explicit determination of the profile, like focusing the focal plane on the surface of the sample [3,9]. Another method involves normalizing the OCT signal by a signal of a reference sample with known attenuation coefficient, in theory effectively canceling the dependence on the confocal effect [8,[22][23][24]. However, this requires an additional reference sample and is only valid if the position of the focal plane is identical between the sample of interest and the reference sample, which is especially critical for high-NA lenses but not always possible. A recent study suggested to measure both the position of the focal plane and the confocal profile through curve-fitting with a simple model, but it still uses the above inaccurate assumption of equal shift and furthermore can only produce a single-value estimate, as such data without averaging is too noisy to allow for the focal plane to vary with (x,y) positions [8,20]. The single-value determination after averaging all A-scans of a volume into a single A-scan, based on the assumption of flat focal plane, can be problematic as the shape of the confocal profile becomes distorted due to the focal plane often lying at an angle within the sample and/or having a curvature when high-NA lenses are used. Furthermore, averaging is not appropriate if the sample is heterogeneous or consists of various structures or high-resolution 2D mapping is desired.
These issues have led us to develop a method to produce high-resolution 2D maps of the attenuation coefficient by accurately determining the positions of the focal plane within a sample while considering the confocal effect and allowing for curved focal planes, thus being applicable to OCT data taken with both low-and high-NA lenses. This method allows for either a single-valued estimate of the focal plane position or a 2D map of the focus positions to account for the curvature and tilt within the sample. We compare these methods to the standard method [20,24] by calculating attenuation coefficient maps of standard samples of Intralipid 20% and 5% using both high-NA and low-NA lenses. Lastly, we apply the newly-developed method to produce an attenuation coefficient map of a xenograft glioblastoma in mouse, demonstrating applicability to real biological data.

OCT system
All OCT measurements were collected with a commercial SD-OCT system (Thorlabs, Newton, NJ, USA). The system uses a large-bandwidth near-infrared light source with center wavelength of 1310 nm, and wavelength bandwidth of 170 nm, which leads to a high axial resolution of 3.5 µm. The system uses a high-speed 2048-pixel line-scan camera to achieve 147,000 A-scan/s with a relatively large imaging depth (2.5 mm maximum). Two objective lenses were tested: low-NA lens (LSM03, Thorlabs) with the focal length of 36 mm, Rayleigh range (z R ) of 135 µm, and the transverse resolution of 25 µm with the 1310-nm wavelength (both in air); and high-NA lens (10X Mitutoyo Plan Apo NIR Infinity Corrected Objective, Mitutoyo) with the focal length of 20 mm, Rayleigh range of 40 µm, and the resolution of 5 µm. The experimentally-determined values of the NA are 0.053 and 0.27 for the low-NA and high-NA lens, respectively.

Scanning and reconstruction of OCT data
When using the low-NA lens, we scanned samples with the lateral sampling interval of 3 µm, while a lateral sampling interval of 1.5 µm was used for the high-NA lens. This volumetric scan was repeated 2 times for later volume averaging to reduce speckle noise. The lens was moved a distance of about 0.1mm (such that the focal plane was still within the sample) and imaging process repeated. Acquired fringe data were reconstructed to 3D images of the OCT intensity following the widely-used method that involves numerical dispersion compensation, k-space interpolation, and inverse Fourier transformation [25]. The reconstructed 3D images underwent a preprocessing to reduce noise: median filtering (kernel: 3x3x3 voxels), volume averaging, and then Gaussian filtering (kernel: normalized Gaussian with 1x1x1 voxel widths in σ and 5x5x5 voxels in window size).

Standard samples
To assess the performance of the proposed and existing methods, we test the methods on solutions of Intralipid with concentrations of 20% and 5% (v/v). Intralipid is a fat emulsion and standard sample that has been used extensively within literature because it has attenuation coefficients in the same range as biological tissue [26,27]. The attenuation coefficients are well described in literature [28] and are referred to in this paper as µ r e f t , such that µ r e f t = 4.7 mm −1 for Intralipid 20% and µ r e f t = 2.3 mm −1 for Intralipid 5%. Using each lens, for each concentration of Intralipid (20% and 5%), two datasets were taken such that the datasets were related by a vertical translation of the lens of about 0.1 mm. The vertical translation was controlled by a motorized translation stage (Thorlabs, K-Cube DC Servo Motor Controller) which was integrated with the OCT system to position the OCT probe.

Xenograft glioblastoma animal models
The mouse model for glioblastoma involved the growth of human-derived tumor in immunosuppressed mice for its superior predictability in location and expected growth over genetic models as well as the selection of cells that constitute the tumor in question. These xenograft models grant said consistency in tumor development and location without greatly sacrificing fidelity to the disease as these models retain invasive growth patterns mimicking human glioblastoma. Additionally, these injections require the use of immunocompromised "nude" mice, without which the tumor would never develop due to the body's immune system recognizing and attacking the foreign cells. As such, this absence leads to the loss of an analogous immune response. However, the replicability and certainty in tumor location and development from onset grants advantages that offset these costs. Sedation of 8-week old Nu/J mice was initiated with 4% isoflurane and continued with 2% isoflurane for the duration of the intracranial injection of human glioma stem cells (hGSCs) via a burr hole in the skull several weeks prior to imaging using a stereotaxic apparatus in collaboration with the Tapinos Lab. After leveling the skull, a hole was drilled with a #72 micro drill bit (Kyocera) at coordinates -2.0 mm AP and +1.5 mm ML relative to Bregma. A 75 RN Hamilton syringe was then lowered to a depth of -2.5 mm DV at a rate of .5 mm per minute, and 200,000 primary hGSCs resuspended in a total volume of 4 µl were injected at a rate of 0.5 µl per minute using a Stoelting Quintessiential Stereotaxic Injector. To reduce backflow, the syringe rested for 5 minutes post-injection, before it was withdrawn at a rate of 0.5 mm per minute, and the cavity was immediately sealed with bonewax (Ethicon). The mice revived with additional oxygen and aged up to four weeks before terminal imaging via a cranial window.

Animal preparation and imaging
The animal was administered 2.5% isoflurane. An area of the skull overlying the somatosensory cortex was shaved, and surgical scrub was applied over incision site to prevent contamination. A 2 cm midline incision was made on the scalp and the skin and both side of the midline was retracted. A 6 mm area of the skull overlying the region of interest of the cortex was thinned with a dental burr until transparent (100 µm thickness). The thinned skull was carefully removed over the region of interest by cutting the thinned skull with scissors, keeping the dura intact. A well was created around the region of interest using dental acrylic, filled with a 1% agarose solution in artificial CSF at the room temperature, and covered with a small glass cover slip. With this window sealed by agarose and covered by the cover slip, the brain was protected from contamination and dehydration. A metal frame was glued by dental acrylic to the exposed skull, used to hold the head and prevent motion during microscopy measurements. After surgery, the animal was placed into the OCT system for in vivo imaging. The animal was kept under the isoflurane anesthesia during imaging. We used the scanning parameters as detailed above. Oxygen saturation, pulse rate, and temperature were continuously monitored with pulse oximetry and a rectal probe during the whole surgical procedure and experiment. The body temperature was maintained at 37°C and the pulse rate remained within the normal range of 250-350 pulse/min. All animal-based experimental procedures were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Brown University and Rhode Island Hospital, according to the guidelines and policies of office of laboratory animal welfare and public health service, National Institutes of Health.

Determination of focal plane: three approaches
Traditionally, single back-scattered light detected with OCT is modelled using Beer's law [16], which takes the form of exponential decay with the attenuation coefficient as the decay constant, where I(z) is the intensity of light at depth from the surface of a sample, µ t is the attenuation coefficient, and h(z) is a function known as the confocal profile which describes the effect of confocal gating on the OCT signal [19]. For Gaussian beams, the confocal effect is described as where z f is the position of the focal plane, z R is the Rayleigh range, and n is the refractive index. The factor of α is incorporated to account for the effects of scattering on the Rayleigh range within the sample, and its value α = 2 has been derived and widely used under the assumption that the interaction of the scatterers with the incident beam forms a secondary beam with new beam waist and Rayleigh range [19]. Although this theoretically estimated value can slightly differ from empirically reported values depending on various factors like the NA and incident angle [16], we used the theoretical value as it is a reasonable estimate for single-scatterers and positions of focus close to the surface of the sample [16,20]. For analyzing data of the standard samples (Intralipid), n= 1.35 was used.
To determine z f , the proposed method involves scanning two images of the sample of interest, I 1 and I 2 , such that I 1 is related to I 2 by a small vertical shift (on the order of 0.1 mm). Once the surface of the sample is identified from each image, the images are coregistered for subsequent analysis, often by vertically shifting both and such that the surfaces in the two images are aligned. Dividing the image I 1 by the I 2 image then yields I d which is simply a function of the two confocal effects, where z r is the apparent Rayleigh range (z r = αnz R ), ∆z is the vertical shift, and the coefficient β accounts for the possibility that the shift of the focal plane differs from ∆z. Most of the previous methods, either comparing two B-scans or moving the focus from the sample surface by a known distance, assumed that this focus shift is identical to the sample (or lens) shift [21,29]; however, in this paper we argue that it is not true in most biological samples with different refractive indices from the air. We find that even when the refractive index is taken into account when estimating the focal shift (even though the refractive index is another unknown in many cases), the results are not significantly different (see Section 4 for details), such that it is optimal to allow this parameter to be a variable. To compare our newly-developed approaches to the previous methods, we implemented the assumption of equal shift in Method 1 by adopting β = 1. Method 1 also adopts the strategy of averaging all A-scans (as commonly done) and then fitting Eq. To reduce noise in the raw 2D map of z f , we fitted the map to a 2D second-order polynomial (S(x, y) = const. + ax + by + cx 2 + dxy + ey 2 ) and used the smooth curved plane in further analysis. When fitting Eq. (3) to I d to determine z f , the region of fitting was determined as the one with the largest decay in the laterally-averaged I d , identically in all methods. Three methods are illustrated in Fig. 1 and Table 1. Figure 2 further demonstrates the procedure with examples and Fig. 3 shows an example of the curved plane fitting. As can be seen in the example of Fig. 2(a), Method 1 forcing β = 1 generally results in worse fitting (red line) than Method 2 with variable (green line). The fitting can be performed as described by Faber et al. [16]: suppose we are interested in fitting a model f with M adjustable parameters a j to N data points (x i , y i ). The estimates of the model parameters a j is found by minimizing the quantity χ 2 given by: For non-linear models the χ 2 -minimization is an iterative process, implemented by the Levenberg-Marquardt method, using the MATLAB function lsqcurvefit. The goodness of fit is given by R 2 , which is a measure that ranges from 0 to 1. R 2 =1 means that the model perfectly fits the data. Simple ray-optics calculation suggests that β should be approximately identical to the ratio of refractive indices of the sample to the air (about 1.35 in our standard samples), although we show that even when the refractive index is accounted for, the focal plane position can not always be accurately estimated (see Discussion). Empirically, our fitting result from Method 2 leads to β values ranging from 0.8 to 5.5, suggesting a wide range of values around the refractive index of Intralipid (n = 1.35). With employing β as a variable, the coefficient of determination of fitting has been improved overall, from R 2 = 0.50 (Method 1) to 0.98 (Method 2) on average. Fig. 1. Illustration of the pipeline performed in this paper to estimate the confocal profile using the novel methods, and the application of these methods in calculating attenuation coefficient maps.

Determination of 'true' focal plane with known attenuation coefficient
To examine the focal planes determined by the three methods, we determined a theoretical focal plane by using the known attenuation coefficients of Intralipid. First, using Eq. (1), we extracted a theoretical confocal profile h(z) by dividing the OCT image I 1 by the term exp(−2µ t z) with the true µ t value ( Figure 4). Second, theoretical h(z) was smoothened to reduce noise. Then, the theoretical focus (z true f ) was determined as the position of the peak in h(z), for each (x, y) position. We found that a curve-fitting approach was not possible here due to the high noise of each individual A-scan. This procedure allows us to determine the theoretical focal plane as a 2D map and use it as a reference to compare those determined by Methods 1-3. Fig. 4. An example of the procedure to find the theoretical focal plane with the true attenuation coefficient. For this example, the same dataset to Figures 1 and 2 was used. Figure 5 displays the focal planes determined with true attenuation coefficients (yellow), the previous Method 1 (red), our Method 2 (green), and its extension Method 3 (blue), for each combination of different scattering samples and different NAs. Method 1 based on the incorrect assumption of equal shifts not only results in wrong estimates (especially with low NA) but also cannot account for the curvature of focal planes. Although Method 2 also cannot account for the curvature, it works better in most cases. Method 3 produces the closest estimates to the theoretical planes, in all cases.

Attenuation coefficient calculation
For each method, we used the estimated z f value(s) to determine the confocal profile h(z) from Eq. (2). Then, we divided the intensity map I 1 by h(z) to minimize the confocal effect. Finally, the attenuation coefficient was determined by fitting the exponential component of Eq. (1) to each resulting A-scan, producing a 2D map of µ t . The region of fitting for each sample was chosen by examining the log intensity of the mean A-scan. A region was chosen from the surface of the sample which appeared to be relatively linear, ensuring that the region was close to the surface of the sample to reduce the effects of multiple scattering which are amplified in depth. Figure 6 illustrates examples of this process for Methods 1-3.
Repeating the above procedure over all (x,y) positions results in a 2D map of the attenuation coefficient. In Fig. 7, we demonstrate the resulting attenuation coefficient maps produced by Methods 1, 2 and 3, as well as the distribution of µ t . Method 1 produces not only inaccurate average µ t values (see the histograms) but also larger variations (more heterogeneous) in attenuation coefficient from the homogeneous standard samples (see the maps). Methods 2 and 3 produce similarly accurate average µ t values (the center of each histogram is close to the true value) with smaller variations, except for the high-NA data of Intralipid 20% sample. This result suggests that Method 3 results in more accurate µ t measurements when high-NA lenses are used for high-resolution mapping of attenuation coefficient in highly scattering samples. The mean and standard deviations of the calculated attenuation coefficients are presented in Table  2. To quantify the differences in performance between the different methods, we determined the accuracy and precision of each method in each condition,

Attenuation coefficient map of xenograft glioblastoma murine model
Having tested and validated the newly proposed methods for attenuation coefficient measurement, we apply the most accurate method, Method 3, to image the xenograft glioblastoma mouse model in vivo. Two OCT volumetric scans were performed with the low-NA lens and analyzed in the same way to the above standard-sample imaging, producing an attenuation coefficient map as presented in Figure 8. The difference in the attenuation coefficients between tumor and healthy tissue is very distinct in the 2D map. All the algorithms were implemented using MATLAB, and our scripts have been uploaded to Brown Digital Repository with the DOI of 10.7301/Z0PV6HW1 for public access [30].

Discussion
In this study, we have established and evaluated a method to robustly determine the positions of the focal plane within a sample post imaging by using two volumetric scans with a small vertical translation. By taking advantage of the fact that the attenuation coefficients are the same in the two images, the method performs a curve-fitting procedure on the ratio of the two confocal profiles to extract the focal plane position. In this course, the major improvements include that we have corrected the false assumption often made that the focus within a sample shifts in a predictable manner with the vertical translation [20,21,31], and we have allowed for a curved focal plane instead of the previously-used flat plane [8,[22][23][24]. When compared with the Method 1, our method results in focal planes much closer to the theoretical one determined with true attenuation coefficients. We further show that even when the refractive index is accounted for, it is difficult to achieve an accurate estimation of the position of focus (Table 4). It is likely that confocal effect and the behavior of light within the sample is dependent on several factors in addition to the refractive index that are not captured in theoretical models, such as the angle of incidence that the beam makes with the sample, the contribution of specular reflection, the scattering coefficient of the substance, the NA of the lens, etc. which may collectively interact to contribute to a wide spread in the β values obtained. Empirically, we have discovered that allowing beta to vary accounts for some of these effects. Finally, the method has been shown to produce more accurate 2D maps of the attenuation coefficient in the standard samples, improving the accuracy by 34 percent points and the precision by 11 percent points. Additionally, the method provides significant improvements on data taken with high-NA lenses, in which the confocal effect contributes significantly to the modulation of the depth profiles. This may have implications for clinical use, where it is desirable to use high magnification lenses to examine tissue microstructures [24].
This study suggests that the accurate determination of the focal plane position is an important factor in precise measurement of the attenuation coefficient. We found that Method 3 which uses a 2D map of the focal plane performed significantly better than Method 1, which assumes that the focal plane shifts within the sample by a predetermined distance as defined by the lens distance. Additionally, averaging of A-scans poses a problem since it is difficult to extract a single-value estimate of the position of the focal plane from the averaged A-scans of a homogeneous sample because often the focal plane is curved or slanted within the sample, such that h(z) is distorted and Eq. (1) is not an appropriate description of the data. If instead the sample were homogeneous and the focal plane perfectly flat within the sample, simple averaging would be theoretically possible. However, if the sample is heterogeneous (as most biological tissues are) simple averaging is once again not appropriate, even when the focal plane is flat, because the attenuation coefficient varies at different (x, y) positions, i.e. the second term in Eq. (1) is no longer appropriate.
By determining a 2D map of the focal plane, the attenuation coefficient maps of the Intralipid solutions are more homogeneous, as expected. We note that it appears that when lateral averaging is performed before curve-fitting to obtain a single µ t high accuracy may be achieved as demonstrated in several papers [3,6,16,28]. Despite this, the method fails if a 2D attenuation coefficient map is desired, as the results in Table 3 indicate. The results are particularly improved in the high-NA case, in which the Rayleigh range is small thus making it critical to accurately specify the position of focus.
The presented paper, however, does not estimate the apparent Rayleigh range empirically which is known to be influenced by factors such as the scattering coefficient, the angle of incidence with the sample, and multiple scattering at deeper depths [16]. For the most accurate calculations, the dependence of the effective Rayleigh range on these factors should be characterized or otherwise determined empirically. Additionally, compared to the previous methods that measure the attenuation coefficient from single OCT volume data [3,4,9,11,11,16], our method requires two scans with a small vertical shift, which increases the total scanning time twice and may increase the total data processing time. However, the vertical shift can be easily implemented using a motorized translation stage, and furthermore, compared to the other method that required a reference sample and identical focal depth between two samples [8,22], our method is much easier to implement and may result in more robust measurements.
There is a growing interest in producing and utilizing 3D maps of the attenuation coefficient [21,32]. For example, depth-resolved attenuation coefficient maps have already been employed to investigate various tissues, including atheroscletertic [33], brain [20], and esophageal [34]. These depth-resolved methods are still vulnerable to the limitations discussed previously, namely the restriction to low-NA lenses and difficulty in precisely specifying the position of focus. Once the confocal effect is accurately described using the proposed method, any method of determining the attenuation coefficient may be applied, producing more accurate estimates of the attenuation coefficient, and also allowing for the use of high-NA lenses. The method described in this paper to determine the position of the focal plane can be applied in conjunction with more complex methods of determining the attenuation coefficient, such as 3D methods or methods that include sensitivity roll-off which describe the reduced signal-to-noise ratio with depth. Sensitivity roll-off is incorporated into the model of intensity as, I(z) = S(z)h(z)exp(2µ t z) where S(z) is the sensitivity roll-off, which was described by Vermeer et al. as a Gaussian [32], and modeled more comprehensively for SD-OCT by Yun et al. as an explicit parameter of spectrometer and CCD parameters [35]. Once sensitivity roll-off is characterized for a OCT system, any images obtained may simply be divided by this factor since there is no dependence on the sample being imaged. Finally, this paper demonstrates an example of applying the developed method to real biological data, producing an attenuation coefficient map of glioblastoma in a murine xenograft model, showing improved contrast of cancerous and non-cancerous tissue compared to the log mean intensity. Considering that the attenuation coefficient is gaining popularity as a tool in various medical applications, like drug development in human tumor spheroids and intraoperative tumor resection [13,36], it is becoming increasingly important that the measurements are robust and accurate. Although this paper has thoroughly evaluated the new method with standard samples and showcased its applicability with the preliminary in vivo tumor imaging, more studies are required to translate the presented method into clinical applications.

Conclusion
This paper presents a method to accurately determine the confocal effect and its focal plane, as an either flat or curved plane, for OCT measurements of the attenuation coefficient. This is achieved by examining the ratio of confocal effects of two OCT images related by a small vertical translation. This paper also shows that more accurate and precise attenuation coefficient maps can be obtained if we avoid the assumption that the focal plane moves equally with movement of the lens, as demonstrated using standard samples of Intralipid. Lastly, this paper exemplifies how this method may be applied to real biological data of a murine glioblastoma model system.

Disclosures
The authors declare that there are no conflicts of interest related to this article.