Generation and optimization of superpixels as image processing kernels for Jones matrix optical coherence tomography

Arata Miyazawa; Young-Joo Hong; Shuichi Makita; Deepa Kasaragod; Yoshiaki Yasuno

doi:10.1364/BOE.8.004396

1. Introduction

Optical coherence tomography (OCT) possesses high resolution and high acquisition speed, and three-dimensional volumetric imaging and video rate monitoring capabilities. OCT has been applied in scientific, industrial, and medical fields [1, 2]. Particularly in ophthalmology, OCT has been widespread and is becoming an essential tool in the diagnosis and monitoring of human retinal disease [3].

An OCT signal intensity image provides layered structures of the retina and helps in the accurate diagnosis of retinal disease. Polarization sensitive OCT (PS-OCT) provides not only structural information but also the polarization properties of the sample [4, 5]. Among the polarization properties, local birefringence (BR) is considered to provide collagen contrast [6–13]. The degree of polarization uniformity (DOPU) is another polarization property [14–17]. Low DOPU signal is considered to be an indicator of melanin [18]. OCT angiography (OCTA) is another extension of OCT, that provides vascular information through the time variation analysis of OCT signals [19–22]. Jones matrix-based PS-OCT, the so-called multifunctional Jones matrix OCT (JM-OCT), provides the following four types of OCT images: scattering intensity, BR/phase retardation, DOPU, and OCTA using a single scan [10, 12, 23–25].

Local statistics are frequently used for the image processing of JM-OCT. For example, the local mean is used to reduce noise and speckle. Similarly, DOPU is computed as a circular variance of Stokes vectors in a local region [14]. BR is also estimated using signals in a small local region [8, 13]. Conventionally, these local statistics are computed using a fixed-size rectangular kernel. However, the fixed-size kernel results in a trade-off between image sharpness and statistical accuracy, that is, a larger kernel provides better statistical accuracy but reduces the image resolution.

We introduce clusters of pixels with a flexible shape, so-called superpixels, as the kernel for computing local statistics. A superpixel is formed with image pixels that share similar signal values and possess high spatial proximity [26]. Using the superpixel as the local statistics kernel, we can preserve tissue structures and simultaneously achieve accurate statistics.

In this study, we present a superpixel method that is based on the idea of the SLIC algorithm [27, 28] but is specially designed for multifunctional JM-OCT. The SLIC algorithm generates superpixels by clustering pixels based on their spatial proximity and color similarity. Our JM-OCT superpixel method generates superpixels by clustering pixels based on their spatial proximity and optical feature similarity where the optical features include OCT intensity, BR, DOPU, and OCTA. The performance of JM-OCT superpixels is evaluated in detail for datasets obtained from in vivo human posterior eyes. Systemic methods for optimizing parameters used for superpixel generation are also presented.

2. JM-OCT system and measurement protocol

In this study, we used a multifunctional JM-OCT designed for posterior eye imaging [24]. A MEMS-based wavelength sweeping light source (Axsun Technology Inc., MA) with a center wavelength of 1.05 μm was used. The scanning rate of the light source was 100 kHz and the average output power was 30 mW. The optical power on the sample was configured to be approximately 1.15 mW to satisfy the safety standard defined by ANSI [29].

This JM-OCT multiplexed two incident polarization states using passive polarization delay, and two output polarizations were measured using a polarization diversity detector [9, 10, 12, 18, 23, 24]. Thus, it measured a set of four OCT images by a single scan. The four images formed a Jones matrix, which is a similarity transformed matrix of the round-trip Jones matrix of the sample. Additionally, four Jones matrix B-scans were repeatedly obtained at a single location on a sample. More details of this JM-OCT are described in Ref. [24].

Four types of optical features, OCT intensity, BR, DOPU, and OCTA, were then computed from the Jones matrix. OCT intensity was computed by coherently combining four entries of the Jones matrix and also combining four repeated Jones matrices [24]. BR was computed using a local Jones matrix analysis method [7] and maximum a-posteriori BR estimator [8]. DOPU was computed using a DOPU algorithm with Makita’s noise correction [16]. Here each DOPU value is computed from spatial 3 by 3 pixels times 4 frames times 2 input polarizations. The 4 frames and 2 input polarizations were accounted by the same manner described in Eq. (3) of Ref. [16]. Additionally, OCTA was obtained by complex correlation analysis with noise-correction [22]. More details of the signal processing are summarized in Refs. [10, 24].

The sensitivity of each of the four images was measured to be 91 dB, which is equivalent to 97 dB for conventional non-polarization sensitive OCT [24]. The depth resolution was measured to be 8.5 μm in air, which corresponds to 6.2 μm in tissue. The transversal resolution is around 20 μm, although it varies by the aberration of the eye. Each A-line consisted of 480 pixels, and the axial pixel separation was 4.0 μm.

Maculae and optic nerve heads (ONHs) of four healthy human subjects were measured. Horizontal cross-sections were obtained with a 6.0 mm scan range laterally and 1.9 mm in depth. Sixty-four B-scans were obtained continuously, and each B-scan consisted of 490 A-lines. Four B-scans among the 64 were used for the analysis. The study protocol adhered to the tenets of the Declaration of Helsinki, and was approved by the institutional review board of the University of Tsukuba.

3. Generation and optimization of superpixels for JM-OCT

In this section, we first explain a method to generate superpixels from multi-contrast Jones matrix OCT images (Section 3.1). Then the methods and algorithms to optimize the parameters used to generate the superpixels are described in 3.2.

3.1. Superpixel generation

3.1.1. Generation of initial superpixels

A superpixel is a cluster of image pixels with high spatial proximity and high optical feature similarity. To generate superpixels, we first created initial superpixels that had a fixed size and hexagonal shape. They were regularly distributed in space and covered the entire image region. To generate the initial superpixel, we first selected the interval of initial superpixels S in a pixel unit. The number of superpixels K in an image with N image pixels is computed as K = N/S². Additionally, the number of pixels in each initial superpixel is S². In the present study, S was set to six pixels, and hence the number of image pixels per superpixel was 36 pixels. The number of image pixels per image was 235,200 pixels (490 × 480 pixels), and hence the number of superpixels was 6,533/image.

3.1.2. Definition of distance in feature space

The initial superpixels were reshaped to increase both the optical feature similarity and spatial proximity of the image pixels in each superpixel. For this reshaping, we had to simultaneously evaluate the optical feature similarity and spatial proximity. So, we had to define an image pixel distance in six-dimensional (6-D) feature space, which consisted of four optical features (OCT intensity, OCTA, BR, and DOPU) and two spatial coordinates (lateral and axial).

The first step to define the 6-D distance was to define a distance in 4-D optical feature space. The optical feature distance between two image pixel points a and b is defined as the weighted Euclidean distance:

D_{o} (a, b) \equiv \sqrt{\sum_{i = 1}^{4} w_{i} {(a_{i} - b_{i})}^{2}},

where a and b are position vectors in the 6-D feature space, which can be expressed as a = [a₁, a₂, a₃, a₄, a₅, a₆]^T and b = [b₁, b₂, b₃, b₄, b₅, b₆]^T. The first to fourth entries of the vectors represent the OCT intensity on a linear scale, OCTA, BR, and DOPU, respectively. The fifth and sixth entries are the lateral and axial spatial positions in the pixel unit. w_i is the weight of the i-th optical feature. A method to define an optimal w_i is described in Section 3.2.2. In our implementation, the optical feature values are normalized in the range of [0, 255] and represented in floating point numbers in its software implementation. And it does not affect the final shape of the superpixels.

Similarly, the spatial distance is defined as the Euclidean distance:

D_{s} (a, b) \equiv \sqrt{\sum_{i = 5}^{6} {(a_{i} - b_{i})}^{2}},

Additionally, the total distance D_t in 6-D feature space is then defined as a weighted sum of the optical feature distance D_o and spatial distance D_s:

D_{t} (a, b) \equiv D_{o} (a, b) + \frac{m}{S} D_{s} (a, b),

where m is a weight between the optical feature distance and spatial distance. A larger m results in a larger contribution of spatial distance to total distance. As will be described in the next subsection, the superpixel is reshaped to reduce D_t among pixels within a superpixel. Hence, a larger m results in more spatially compact superpixels, and m is called as the “compactness factor” [27]. A method to determine the optimal compactness factor is described in Section 3.2.3.

3.1.3. Superpixel reshaping

After generating the initial superpixels, they were reshaped in an iterative process similar to the SLIC superpixel method [27]. The reshaping was an iterative clustering process of image pixels in the 6-D distance D_t. For each iteration, the center of gravity of all the image pixels that belonged to each superpixel was computed in the 6-D feature space. This center of gravity is denoted as the centroid of the superpixel. The centroids of all superpixels were computed. Then each image pixel was reassigned to the superpixel whose centroid was the nearest to that image pixel. In our implementation, the searching area of the nearest superpixel was limited to within a 2S × 2S area, which is double the interval of the initial superpixel, to accelerate the searching speed. The iteration continued until the superpixels converged into particular shapes. In practice, the iteration continued until the spatial distance between the previous and recomputed centroids became less than a threshold distance, which was 1 pixel in this study, or reached 10 iterations. This reassignment reshaped the superpixels to reflect the tissue structures.

After this iteration, some superpixels occasionally split into small fragments because this clustering process did not constrain spatial connectivity. Hence, after the convergence, small fragments of superpixels, that is, fragments smaller than four pixels in this study, were merged into the largest neighboring superpixel. This process is referred to as connectivity enforcement.

3.2. Optimization method of superpixel parameters

The superpixel method relies on some arbitrary defined parameters including the initial superpixel interval S, weights w_i, and compactness factor m. Among them, the initial superpixel interval is mainly defined based on the allowable computational load, while the weights and compactness factor should be optimized to make superpixels well adhere to tissue boundaries. In the following sections, we describe strategies and methods to optimize the weights (Section 3.2.2) and compactness factor (Section 3.2.3). To perform these optimizations, we requested quantitative metrics of the goodness of the superpixels. Thus, we first define the metrics in Section 3.2.1 before describing the optimization methods.

3.2.1. Intra-superpixel variance and contribution metric

We assumed that each type of tissue had a specific optical property. Thus, the measured optical feature values from the same tissue had similar values. If a superpixel adheres well to tissue boundaries and contains only a homogeneous tissue, the variance of the optical feature values within the superpixel should be low. Hence, the intra-superpixel variance (ISPV) of each optical feature can be considered as a metric to evaluate how well a superpixel adheres to tissue boundaries.

The variance of the i-th optical feature within the k-th superpixel is defined as

σ_{i, k}^{2} = \frac{1}{s_{k}} \sum_{j = 1}^{s_{k}} w_{i} {(x_{i, j, k} - {\bar{x}}_{i, k})}^{2},

where s_k is the number of pixels within the k-th superpixel; w_i is the weight of the i-th optical feature, which is the same as that previously used for the distance calculation in Eq. (1); x_i,_j,_k is the i-th optical feature value of the j-th pixel in the k-th superpixel; and

{\bar{x}}_{i, k}

is the mean of the i-th optical feature values within the k-th superpixel.

We used the mean of all the variances within each superpixel as a metric for superpixel quality, which is referred to as the ISPV. The ISPV of the i-th optical feature V_i is expressed as

V_{i} \equiv \frac{1}{K} \sum_{k = 1}^{K} σ_{i, k}^{2},

where K is the number of superpixels in the image. Superpixels that adhere well to tissue boundaries obtained by the i-th optical feature should have low V_i.

3.2.2. Optimization of weights for the optical features

We used four optical features to discriminate tissue boundaries. Frequently, the boundary of a particular tissue is clearly delineated with some optical properties; however, this does not appear in other optical property images. For example, blood vessels are clearly visible in OCTA, but are barely visible in BR. Hence, if the superpixel shape is dominated by only a few optical properties, the superpixel does not follow the shape of tissues that appear only in the non-dominant optical properties. Thus, it is a rational optimization strategy to make the four optical properties equally contribute to superpixel formation. In this section, we present a method that controls the weights of the optical properties (w_i) so that all optical properties make an equal contribution.

In our optimization method, the contribution of the i-th optical property was evaluated using the following metric:

C_{i} \equiv {(\frac{V_{i}^{r e s h a p e d}}{V_{i}^{i n i t i a l}})}^{- 1} = \frac{V_{i}^{i n i t i a l}}{V_{i}^{r e s h a p e d}},

where

V_{i}^{i n i t i a l}

and

V_{i}^{r e s h a p e d}

are the ISPVs of the i-th optical property of the initial and reshaped superpixels, respectively. If an optical feature contributes to reshaping, the ISPV of the optical feature is reduced as the reshaping iteration progresses. Thus, ISPV may be useful as a measure of the contribution. However, ISPV is not directly usable for comparison among the optical features because each optical feature has different magnitudes of variances. Although we have applied min-max scaling to optical features into [0, 255] (Section 3.1.2), it has not regulated the variances and standard deviations of the optical features. To overcome this problem, the ISPV was normalized by the ISPV of the initial superpixel shape, as in the middle part of Eq. (6). We assumed that all the optical features were equally random within the initial uniform superpixels. Under this assumption, the normalized ISPV become comparable among the optical features. Finally, a contribution metric C_i is defined as the inverse of the normalized ISPV, that is, a greater contribution from an optical feature results in a higher contribution metric of the optical feature.

To equalize the contributions from the four optical features, the variance of the four contribution metrics was used as a cost function E:

E (w) = \frac{1}{4} \sum_{i = 1}^{4} {(C_{i} (w_{i}) - \frac{1}{4} \sum_{i = 1}^{4} C_{i} (w_{i}))}^{2},

where C_i (w_i) is a contribution metric of the i-th optical feature when the weight of the optical feature is w_i and w is a vector of weights defined as w ≡ [w₁, w₂, w₃, w₄]^T. As the contribution of all the optical features becomes more equal, E(w) becomes smaller. Thus, the optimization of weights was performed to minimize E(w).

We used an iterative method to minimize the cost function. The initial weights were determined randomly and the initial weight vector $w^{(0)} = {[w_{1}^{(0)}, w_{2}^{(0)}, w_{3}^{(0)}, w_{4}^{(0)}]}^{T}$ was normalized to a unit vector. Hereafter, a superscript in brackets represents an iteration index. A new weight vector in the (l + 1)-th iteration is then defined to be

w^{(l + 1)} \propto w^{(l)} + α d^{(l)},

where α was a step length, which controls the optimization speed. It was set to be 0.315 in this study. d⁽^l⁾ is an update of the weight vector defined as

d^{(l)} \equiv - [C^{(l)} - \max {C^{(l)}}],

where C⁽^l⁾ is a vector of contribution metrics at the l-th iteration defined as

C^{(l)} = {[C_{1}^{(l)}, C_{2}^{(l)}, C_{3}^{(l)}, C_{4}^{(l)}]}^{T}

.

\max {C^{(l)}}

is the maximum entry of C⁽^l⁾. w⁽^l⁺¹⁾ is then normalized to a unit vector. Thus, Eq. (8) sufficiently specifies w⁽^l⁺¹⁾, although it is in a proportional form. The weights of the optical features with smaller contributions were increased more in the iteration. This iteration continued until E(w) become smaller than 10⁻⁶ or reached 20 iterations.

Current optimization algorithm equalizes the reduction of variances of each optical feature. However, this equalization does not work correctly if the measured tissue region does not have any structures in one or some of the optical feature images. For example, the algorithm cannot properly optimize the weights of OCTA or birefringence, if a tissue has no perfusion or birefringence, respectively. Fortunately, for the current in vivo applications, such cases are hard to be found. However, if it was found, the optical feature without structure should be excluded from the equalization.

3.2.3. Optimization strategy of the compactness factor

The compactness factor, m in Eq. (3), controls the balance between optical feature similarity and spatial proximity. If the compactness factor is too large, superpixels are not reshaped flexibly and remain almost in their initial hexagonal shape. If it is too small, the superpixels ignore spatial proximity and split into small fragments more frequently. In this case, the small fragments are absorbed into the largest neighboring superpixel in the connectivity enforcement step (as described in the last paragraph of Section 3.1.3). This connectivity enforcement ignores optical feature similarity, and hence, superpixels do not adhere to tissue boundaries well and the ISPV increases.

We had to determine an optimal compactness factor that might provide small ISPVs for all the optical features even after the connectivity enforcement. We optimized the compactness factor to minimize the normalized ISPVs of all the optical features, where the normalized ISPV is normalized by the ISPV of the initial superpixels. Recall that the normalized ISPV is the reciprocal of the contribution metric C_i, and the C_is of all the optical features are nearly the same after optimizing the features’ weights, as described in Section 3.2.2. Thus, the ISPVs of all optical features became similar to each other. We selected a compactness factor so that the average of the ISPVs of all the optical features (mean normalized-ISPV) had a minimum value. Superpixels with a lower mean normalized-ISPV can be considered to adhere to tissue boundaries better. This optimization was performed by brute force optimization, as described in Section 4.1.2.

3.2.4. Flow of optimizations

The entire optimization process is a three-level-nested iterative process, as summarized in Fig. 1. The outermost iteration optimizes the compactness factor (solid red box in Fig. 1) and corresponds to Section 3.2.3. The second-level iteration is within this outermost iteration and optimizes the weights of the optical features (dashed green box, Section 3.2.2). The second-level iteration contains the innermost iteration, which is for superpixel generation (dashed and dotted blue box, Section 3.1).

Fig. 1 Flow of optimization for parameters used in superpixel generation. The parameters to be optimized include the compactness factor and weights of the optical features. Three loops are nested in the optimization process: the solid red, dashed green, and dashed and dotted blue boxes indicate the loops for compactness factor optimization (Section 3.2.3), weight optimization, (Section 3.2.2), and superpixel generation (Section 3.1), respectively.

Download Full Size | PDF

4. Results

This section consists of two subsections. In the first subsection (Section 4.1), we present the validation results of the optimization method of superpixel parameters. In the second subsection (Section 4.2), we present applications of superpixels to in vivo posterior eye images.

4.1. Validation of the superpixel parameter optimization methods

4.1.1. Validation of feature weights optimization

We searched for the optimal weights as described in Section 3.2.2. Fig. 2 shows examples of the weights of each optical feature as functions of the index of the optimization iteration. The three graphs show three trials (trial-1 to -3) of optimization, which were performed for the same ONH image, but started from different random initial weights. The compactness factor was set to five, which was found to be the optimal compactness factor as described in the next subsection. The weights of each optical feature converged after approximately 10 iterations and the converged values of the three trials were almost the same.

Fig. 2 The alteration of the weights of the optical features during iterative optimization in three trials. (a) to (c) correspond to the first to third trials, respectively. Each plot color represents each optical feature as OCT intensity (purple, circle), OCTA (green, triangle), BR (red, square), and DOPU (blue, cross).

Download Full Size | PDF

Fig. 3(a) shows the cost function E(w) of trial-1 at each iteration. It is evident that the cost function successfully became very low and stable after a few iterations. The corresponding contribution metrics of each optical feature are shown in Fig. 3(b), where the colors of the plots of each feature are the same as those in Fig. 2. The contribution metrics of all features were successfully converged to a similar value after a few iterations. Thus, all the optical features equally contributed to superpixel formation with the optimized weights.

Fig. 3 (a) The alterations of the cost function E(w) and (b) the contribution metrics of each optical feature during the iterative optimization of the weights in the first trial. Each plot color in (b) represents each optical feature as OCT intensity (purple, circle), OCTA (green, triangle), BR (red, square), and DOPU (blue, cross).

Download Full Size | PDF

The optimal weight values varied among the measured tissue types, subjects, and measurement sessions. This issue will be discussed in Section 5.3.

4.1.2. Validation of the compactness factor optimization

Fig. 4 shows how the compactness factor (m) affected superpixel reshaping. The first column shows kernel-averaged OCT intensity images created using the superpixels as averaging kernels, and the second column shows superpixel images, where each superpixel is displayed in a randomly selected color. The weights of the optical features were optimized as described in Section 3.2.2 for each compactness factor. All superpixels generated with the compactness factor of 100 retained their initial hexagonal shape, even after reshaping [Fig. 4 (k)]. By contrast, the superpixels generated with the compactness factor 5 had a variety of shapes and sizes, and they adhered to the layered structure of the retina [Fig. 4 (q)]. Generally, we found that the reshaping flexibility of the superpixels increased as the compactness factor decreased. The mean normalized-ISPVs were 0.86, 0.80, 0.67, and 0.55 for m = 100, 70, 40, and 5, respectively.

Fig. 4 Superpixels generated with several compactness factor configurations. The first column ((a), (c), (e), and (g)) shows the kernel-averaged images with superpixel kernels. The second column ((b), (d), (f), and (h)) shows the superpixels where each superpixel is displayed with randomly assigned colors. The third column ((j)–(r)) shows magnified images of the images in the first and second columns, where the magnified regions are indicated by yellow boxes. Each row corresponds to the compactness factor of 100, 70, 40, and 5. The bottom row ((i) and (r)) shows the original (non-superpixelized) OCT intensity images.

Download Full Size | PDF

Fig. 5 shows an example of the mean normalized-ISPVs at each compactness factor, where a better compactness factor demonstrates a smaller mean normalized-ISPV. The mean normalized-ISPVs were computed from the same ONH image as that in Section 4.1.1, and the feature weights were optimized for each compactness factor independently. Generally, the mean normalized-ISPV decreased as the compactness factor decreased. This is because the lower compactness factor that reshapes superpixels is more highly dependent on optical feature similarity than spatial proximity. However, the mean normalized-ISPV increased if the compactness factor become too small. This can be explained as follows: A compactness factor that is too small results in a severe spatial fragmentation of superpixels at first, and then the small fragments are merged into the largest neighboring superpixels in the connectivity enforcement process. Because this merging is performed irrespective of the optical feature similarity, the mean normalized-ISPV becomes large.

Fig. 5 Mean normalized-ISPVs for several compactness factor configurations.

Download Full Size | PDF

We performed this evaluation for two retinal locations (macula and ONH) of two subjects, and found that a compactness factor of five always provided the minimal or nearly minimal normalized-ISPV. Thus, the compactness factor was set to five to obtain the results of Sections 4.1.1 and 4.2.

4.2. Superpixelization of in vivo posterior eye image

Figure 6 shows an example of an ONH. The superpixels were generated with optimal feature weights, which were 0.742 for intensity, 0.239 for OCTA, 0.393 for BR, and 0.487 for DOPU. The compactness factor was set to five. The final mean normalized-ISPV was 0.547.

Fig. 6 An example of an ONH. The first column shows the original images of (a) OCT intensity, (c) OCTA, (e) BR, and (g) DOPU. The second column ((b), (d), (f), and (h)) shows the kernel-averaged images with reshaped superpixels that correspond to the images in the first column. (i)–(k), (l)–(n), (o)–(q), and (r)–(t) show three types of images in the square windows in (b), (d), (f), and (h), respectively; (i), (l), (o), and (r) show the original images; (j), (m), (p), and (s) show the kernel-averaged images with initial (hexagon) superpixels; and (k), (n), (q), and (t) show the kernel-averaged images with reshaped superpixels.

Download Full Size | PDF

By observing the magnified images of OCT intensity [Fig. 6(i), (j), (k)], it is evident that the retinal layer structures that included retinal pigment epithelium (RPE) were well preserved by the reshaped superpixels. The shape of RPE can also be clearly recognized in the superpixelized DOPU image [Fig. 6 (t)]. Vessels in OCTA were clearly delineated by superpixels [Fig. 6(n)]. In the superpixelized BR image, high BR regions of sclera were depicted in the superpixelized BR image [Fig. 6(q)]. Thus, the superpixels were correctly reshaped to represent the tissue structures.

Figure 7 shows an example of a macula, where the alignment of subfigures are the same as those in Fig. 6. The optimized feature weights were 0.691 for OCT intensity, 0.313 for OCTA, 0.407 for BR and 0.508 for DOPU, while the compactness factor was set to five. Similar to the case of the ONH, it can be seen that the reshaped superpixels well adhered to the tissue boundaries. The contrast of the birefringence at deep choroid and sclera is lower in the superpixel-averaged image [Fig. 7(f)] than the original image [Fig. 7(e)]. It is because of the property of averaging of measured phase retardation. As described by Makita et al., the mean of measured phase retardation asymptotically approaches around 2/3 π radians as the effective signal-to-noise ratio (ESNR) of Jones matrix OCT signal decreases, and it is irrespective to the true phase retardation value [7]. Since the ESNR of the deep choroidal and scleral regions is low, the superpixel-kernel averaging erroneously up- or down-shifts the phase retardation and degrades the birefringence contrast. One possible solution is “re-estimation” of birefringence by a statistically relevant estimator in conjunction with the superpixel kernels. It is discussed in more details later in Section 5.6.

Fig. 7 An example of a macula. The first column shows the original images of (a) OCT intensity, (c) OCTA, (e) BR, and (g) DOPU. The second column ((b), (d), (f), and (h)) shows kernel-averaged images with reshaped superpixels. The square boxes in (b), (d), (f), and (h) show the locations of (i)–(k), (l)–(n), (o)–(q), and (r)–(t). (i), (l), (o), and (r) show the magnified original images; (j), (m), (p), and (s) show the kernel-averaged images with initial (hexagon) superpixels; and (k), (n), (q), and (t) show the kernel-averaged images with reshaped superpixels.

Download Full Size | PDF

The optimal feature weights varied between the ONH and the macula. The variation and universality of the optical weights are extensively discussed later in Section 5.3.

4.3. Conventional rectangle kernel and superpixel kernel

In order to demonstrate the advantage of superpixel kernel for local averaging, we compare averaged images with superpixel kernels and a standard rectangle kernel as shown in Fig. 8. Figures 8(a), 8(b), and 8(c) are an original image, moving averaging with a 6 × 6-pixel moving averaged image, and a locally averaged image with superpixels kernels. The mean number of pixels per superpixel is 36, which is equal to the number of pixels in the rectangle kernel. The insets are magnified images of the regions indicated by the yellow boxes.

Fig. 8 Comparison among intensity images of original (a), moving average with a rectangle kernel (b), and local average with superpixel kernels (c). The rectangle kernel has a size of 6 × 6 pixels and the average pixel number per superpixel kernel is 36. The insets are the magnified images of the regions indicated by the yellow box. (a) and (c) are the same with Figs. 6(a) and 6(b), respectively.

Download Full Size | PDF

It is evident that thin layered structures, such as RPE, are well preserved by the superpixel kernels, while it is blurred by the rectangle kernel. It is also noteworthy that the rectangle kernel image is generated by moving averaging, while the superpixel kernel averaging is non-moving averaging. So, the superpixel averaged image needs less information for its numerical representation than the rectangle kernel image, although it contains more clinically interpretable information.

4.4. Superpixels generated only form OCT intensity

One of the advantageous characteristics of superpixel kernel is that it follows the structures appeared in all contrast types. In order to highlight this characteristic, superpixel kernels are generated only form intensity OCT. Its kernel averaging images are compared with the original images and standard superpixel kernel averaged images as shown in Fig. 9. Here the standard superpixel kernels are generated by the method presented in Section 3.1. The first to third columns are the original images, standard superpixel kernel averaged images, and intensity-only superpixel kernel averaged images.

Fig. 9 The comparison of superpixel generated by the presented method (standard superpixel) and that generated only from intensity OCT. The first to third columns show original images, standard superpixel kernel averaged images, and intensity-only superpixel kernel averaged images, respectively.

Download Full Size | PDF

Although the structures appeared in the original intensity image [Fig. 9(a)] are more preserved in the intensity-only superpixel kernel averaged image [Fig. 9(c)] than the standard superpixel kernel averaged image [Fig. 9(b)], the structures appeared in the other contrasts are more clearly appeared in the standard superpixel kernel averaged images [Figs. 9(e), 9(h), and 9(k)] than intensity-only superpixel kernel averaged images [Figs. 9(f), 9(i), and 9(l)].

5. Discussions

5.1. Speckle effect on superpixels

OCT tomograms exhibit speckle. So, OCT intensity is broadly distributed in a negative exponential or a close relative of Rayleigh distributions as far as the speckle is fully developed [30]. Since the means of this distributions monotonically increase as the back scattering intensity increases, the superpixels can follow the tissue structures in the intensity OCT image despite of the broadness of the distributions. But it also can degrade the performance of superpixelization. Additional preprocessing for speckle reduction may improve the performance of superpixelization.

5.2. Optimal definition of spatial distance

The pixel size of the images in this study was 4.0 μm (axial) × 12.5 μm (lateral). Additionally, the spatial distance was expressed as the multiplication of the pixel size in the superpixel generation algorithm. Hence, the weight of the axial distance was approximately three times larger than that of the lateral distance. According to the final superpixel images, this imbalance between the weights of the spatial distances was acceptable. This could be partially because the retina has a layered structure, and hence it would be reasonable to apply a larger weight to the axial direction, which consists of a finer structure, than to the lateral direction.

5.3. Universality of the optimized weights of the optical features

In our algorithm, the weights of the four optical features were optimized for a specific single B-scan using the methods described in Section 3.2. We discuss how universally a specific set of weights optimized by a single B-scan is applicable to other B-scans. Universality is important because if the parameter set is universal, to some extent, the optimization is required only once for some range of B-scans. We discuss three types of universality. The first is intra-dataset universality, which is universality of the optimized weights among B-scans continuously obtained by a single measurement. The second is inter-dataset and intra-subject universality, where the weights are optimized and reapplied to datasets of the same subject but of different measurements. The third is inter-subject universality, which is the applicability of the optimized weights to a dataset from a different subject.

Optimization was performed to make the four optical features equally contribute to superpixel generation. Hence, the variance of the contributions among the optical features was used as a metric to evaluate the eligibility of a weight set, where the contribution is a quantity C_i defined by Eq. (6). This variance of the contributions becomes small if the weights are eligible for the B-scan. We consider that a set of weights is qualified to apply to a B-scan if the variance of the contributions is smaller than 0.01. As a reference, in Fig. 3 (b), the variance of the contributions became smaller than 0.01 after two iterations (0.592, 0.0126, and 0.00245 at the zeroth, first, and second iterations, respectively).

Additionally, the normalized-ISPV should be small for all optical features if the weights are eligible. Thus, the largest normalized-ISPV among the four optical features was used as a second metric to evaluate the weight set. A smaller maximum normalized-ISPV indicates that the weight set is more eligible for the B-scan. As a reference, we also generated superpixels without optimization because all weights were set to unity, and compared the normalized-ISPV to the optimized cases.

We computed these two metrics for B-scans taken from four ONHs and four maculae of four subjects and evaluated the three universalities. All subjects were 30 to 40 years old East Asians. Subjects 1 and 2 are the same subjects as those in the previous sections.

To evaluate the intra-dataset universality, we computed a set of optimal weights from a B-scan (training B-scan) and applied it to the other 3 B-scans (test B-scan) taken with the same acquisition sequence as the training B-scan. The results of the ONH are shown in Table 1. The variances of the contributions were smaller than 0.01 (qualified) for all test B-scans. All the largest normalized-ISPVs were also smaller than those without optimization; that is, all the largest normalized-ISPVs of the test B-scans with optimization were smaller than 0.578 and those without optimization were larger than 0.673. The normalized-ISPV was evidently improved by optimization. The same comparisons were performed in the macula, and the results are shown in Table 2. As with the ONH, the variances of the contributions were smaller than 0.01 for all test B-scans and the largest ISPVs were also smaller than those without optimization for all test B-scans. Hence, the optimal parameters can be considered to work well in other B-scans of the same measured volume.

Table 1. Validation for intra-dataset universality examined with an ONH of subject-1.

View Table | View all tables in this article

Table 2. Validation for intra-dataset universality examined with a macula of subject-2.

View Table | View all tables in this article

For a reference, the optimal weights for the Training in Table 1 were 0.797 (intensity), 0.207 (OCTA), 0.345 (BR), and 0.451 (DOPU). Those for Table 2 were 0.779 (intensity), 0.280 (OCTA), 0.366 (BR), and 0.425 (DOPU).

For the evaluation of the inter-dataset and intra-subject universality, we first computed optimal weights from a B-scan (training B-scan), and applied it to a test B-scan, which was taken from another measurement session, but the same subject. The results were also compared with those with unoptimized results, where all the weights were set to unity. The same tests were performed for three subjects. The results of the ONH are summarized in Table 3, where VOC denotes variance of the contributions and nISPV denotes normalized ISPV. The variances of the contributions were smaller than 0.01 in all test B-scans and the largest ISPVs were also smaller than those without optimization for all test B-scans. The same comparisons were performed in the macula, and the results are shown in Table 4. The variances of the contributions were smaller than 0.01 for all test B-scans. However, the largest ISPV of test-3 with optimization was 0.609, and was not substantially improved compared with that without optimization (0.627). We conclude that there was inter-dataset and intra-subject universality, to some extent; however, it was not always the case. Thus, it would be safe to optimize the weights for each individual measurement dataset.

Table 3. Validation of inter-dataset-and-intra-subject universality. The validation was performed with ONH datasets. VOC denotes variance of the contributions and nISPV denotes normalized-ISPV. “Optimal” indicates that the weights were optimized and w/o optimization indicates that the weights were not optimized but unities.

View Table | View all tables in this article

Table 4. Validation of inter-dataset-and-intra-subject universality. The validation was performed with macular datasets. VOC denotes variance of the contributions and nISPV denotes normalized-ISPV. “Optimal” indicates that the weights were optimized and w/o optimization indicates that the weights were not optimized but unities.

View Table | View all tables in this article

Training-1 in Table 3 and Training-1 in Table 4 were the same with the Trainings in Tables 1 and 2, respectively. So, these have the same optimal weights. The optimal weights of Training-2 of Table 3 were 0.799 (intensity), 0.219 (OCTA), 0.3412 (BR), and 0.444 (DOPU). Those of Training-3 were 0.839 (intensity), 0.209 (OCTA), 0.300 (BR), and 0.403 (DOPU). Those of Training-2 in Table 4 were 0.645 (intensity), 0.384 (OCTA), 0.410 (BR), and 0.517 (DOPU), while those of Training-3 were 0.623 (intensity), 0.332 (OCTA), 0.483 (BR), and 0.518 (DOPU).

For the inter-subject universality evaluation, we computed the optimal weights from a B-scan (training data) and applied them to the three B-scans (test data) of the other three subjects. We also computed the largest normalized-ISPV with unoptimized weights, all weights were set to unity for comparison. The results of the ONH are summarized in Table 5, where VOC denotes variance of the contributions and nISPV denotes normalized ISPV. The variances of the contributions were smaller than 0.01 in all test B-scans and the largest ISPVs were also smaller than those without optimization in all test B-scans. The same comparisons were performed in the macula, and the results are shown in Table 6. The variances of the contributions of Test-1 and Test-2 were larger than 0.01. The largest ISPV of Test-2 with optimization was 0.614, and was not substantially improved compared with that without optimization (0.641). Hence, we conclude that a set of optimized weights did not always work well in B-scans of other subjects. This lower universality among subjects could be partially explained by the inter-subject variation of the optical parameters of tissues. For example, the melanin concentration in choroid varies by age [31], which results in age-related variation of DOPU in the choroid.

Table 5. Validation of inter-subject universality. The validation was performed with ONH datasets. VOC denotes the variance of the contributions and nISPV denotes the normalized-ISPV. “Optimal” indicates that the weights were optimized and w/o optimization indicates that the weights were not optimized but unities.

View Table | View all tables in this article

Table 6. Validation of inter-subject universality. The validation was performed with macular datasets. VOC denotes the variance of contributions and nISPV denotes the normalized-ISPV. “Optimal” indicates that the weights were optimized and w/o optimization indicates that the weights were not optimized but unities.

View Table | View all tables in this article

Trainings in Tables 5 and 6 were the same with the Trainings in Tables 1 and 2, respectively. So, these have the same optimal weights.

To summarize, a specific set of optimal weights was applicable among the B-scans in the same measurement dataset. However, it was not always applicable to other datasets and other subjects.

5.4. Computational time

The computational time of superpixel generation was dominated by three aspects: 6-D distance calculation, centroid calculation, and connectivity enforcement process.

The distance calculation measures the distance between the centroids from each superpixel to each image pixels. The distance calculation was performed approximately 4N-times, where N is the number of image pixels. The factor of four was selected because the distance computation was performed only for the image pixels, which was in an area with a size of 2S × 2S centered at the centroid, where S is the interval among the initial superpixels. Thus, the distance was computed approximately four times on average for each image pixel. According to Eqs. (1)–(3), the computational time for a single distance calculation is approximately estimated as

\begin{array}{l} τ_{D} = [(τ_{a} + 2 τ_{m}) n_{o} + (n_{o} - 1) τ_{a} + τ_{s q r t}] + [(τ_{a} + τ_{m}) n_{s} + (n_{s} - 1) τ_{a} + τ_{s q r t}] \\ = 2 (n_{o} + n_{s} - 1) τ_{a} + (2 n_{o} + n_{s}) τ_{m} + 2 τ_{s q r t}, \end{array}

where τ_m, τ_a, and τ_sqrt are the computational times for single multiplication, single addition, and square root operations, respectively. n_o and n_s are the numbers of optical features and spatial dimensions, respectively. The first and second pairs of square brackets in the first line represent the computational times of Eqs. (1) and (2), respectively. The first terms of each part correspond to the computational time for each operation in the summation (∑), and the second terms represent the summation itself. The computational time for the total distance calculation is estimated to be 4Nτ_D.

The total computational time for the centroid τ_C is estimated to be

\begin{array}{l} τ_{C} = (n_{o} + n_{s}) [\sum_{k = 0}^{K - 1} (N_{k} - 1) τ_{a} + K τ_{d}] \\ = (n_{o} + n_{s}) [(N - K) τ_{a} + K τ_{m}], \end{array}

where N_k is the number of image pixels in the k-th superpixel and τ_d is the unit computational time for a single division. To derive the second line, we assumed that the unit computational times for single division and single multiplication were the same, that is, τ_d = τ_m.

The computational time for connectivity enforcement is highly dependent on its algorithm. Thus, we consider it as a black box, and denote it as τ_CE.

The total computational time for superpixel generation is approximately estimated as

τ_{t o t a l} ≃ [4 N τ_{D} + τ_{C}] ι + τ_{C E},

where ι is the number of iterations in the superpixel reshaping process. Although ι nonlinearly varies by the compactness factor, it can be regarded as a constant for the realistic values of the compactness factor. By substituting Eqs. (10) and (11) into Eq. (12), and assuming τ_a is negligibly small and τ_a ≪ τ_m, Eq. (12) becomes

τ_{t o t a l} \sim {[(8 n_{o} + 4 n_{s}) τ_{m} + 8 τ_{s q r t}] N + (n_{o} + n_{s}) τ_{m} K} ι + τ_{C E} .

By assuming τ_sqrt ≃ τ_m and substituting the parameters of our particular study, n_o = 4 and n_s = 2 into Eq. (13), the total computational time is approximately, but finally, estimated as

τ_{t o t a l} \sim ι (48 N + 6 K) τ_{m} + τ_{C E} .

This approximate estimation suggests that the number of image pixels (N) has approximately an eight times higher impact on the computational time than the number of superpixels (K) by omitting τ_CE.

As approximately, but theoretically, estimated above, the computational time of the superpixel generation algorithm varies by the number of superpixels, compactness factor, and particular algorithm implementation. In this study, the number of superpixels was 6,533 and the compactness factor was five. In our implementation, the core part of superpixel generation, which was indicated in the dashed and dotted blue square in the optimization chart (Fig. 1), took approximately 65 s for a B-scan. The weight optimizing process, including superpixel generation (dashed green box in Fig. 1) took approximately 22 min. Because weight optimization includes the superpixel generation process, the entire optimization (solid red box) took 22 min × n_i, where n_i is the number of iterations required to determine the optimal compactness factor. The computation was performed using an Intel Core i7 4219HQ CPU with a clock frequency of 2.3 GHz.

5.5. Image segmentation and superpixel

The majority of current retinal OCT segmentation methods uses layered structural property of retina. It is useful and accurate as far as the layered structure is preserved. However, it is hard to be applied for severe pathologic cases. For example, the RPE of a highly pathologic eye often does not preserve the layer structure, and hard to be segmented by the current segmentation methods.

This problem can be partially solved by polarization sensitive OCT. For example, Götzinger et al. selected RPE pixels by its DOPU value [14]. Such pixel-by-pixel classification method can be further improved by multifunctional JM-OCT. JM-OCT provides multiple optical features of each pixel. Due to these optical features, each pixel in a JM-OCT image can be projected into multi-dimensional feature space. It will enable application of machine learning based pixel-wise classifiers, such as support vector machine, to retinal tissue segmentation.

For general image processing, superpixel methods have been originally developed for accelerating the machine learning based segmentation by reducing the number of pixels [28]. For example, the learning time and classification time of standard classifiers, such as support vector machine, are proportional to the number of training and test data, i.e., pixels. By using the superpixel, we can reduce the number of training/test data from the numbers of pixels to that of superpixels. So, a few ten times acceleration is expected.

On the other hand, the specific motivation of the present study is to have an adaptive kernel to improve the local statistics. The improved statistics can also improve the accuracy of segmentation of JM-OCT [32]. In addition it also can accelerate machine learning based JM-OCT segmentation as it was intended by the original motivation of superpixel development. So, superpixel method would improve the segmentation accuracy of JM-OCT as well as its computation speed.

5.6. Re-estimation of optical properties by superpixel kernels

In this study, we have used the superpixel as an averaging kernel of birefringence, which is a constantly scaled version of phase retardation. However, Makita et al. found that the mean of measured phase retardation does not asymptotically approach the true phase retardation even the number of measurement becomes huge [7]. And hence, the averaging of the phase retardation, and also of birefringence, is not really rational.

One possible solution is combination of maximum a-posteriori (MAP) birefringence estimator [8, 13] and the superpixel kernel. Since our particular method has used the MAP estimator for the initial estimation of birefringence, this secondary estimation process can be considered as “re-estimation” of birefringence with an adaptive kernel. Our preliminary study showed that the re-estimation provides more rational estimation of retinal birefringence [33]. Not only the MAP estimator, other estimation and averaging methods, such as raw-OCT-intensity averaging [34], Jones matrix averaging [24, 35], Cloude-Pottier-decomposition based Jones matrix estimation [36], and other methods (see “advanced data processing” section of Ref. [5]), also can be used with the superpixel kernel.

The re-estimation is potentially applicable also to OCTA and DOPU. However, both of them are associated with the randomness of the OCT signal. So, the structural randomness within a kernel might affect these quantities. Superpixel kernels have variety of shapes and sizes, and hence each kernel has different degrees of structural randomness. So, its application to OCTA and DOPU re-estimation is not straight forward and it would be a future challenge.

6. Conclusion

We developed a new superpixel method specially tailored to multifunctional JM-OCT. A systemic optimization method for parameters in the superpixel algorithm was presented. The performance of the optimization method was evaluated in detail and it was found to work correctly. The superpixel method was applied to retinal OCT, and the generated superpixels were found to well preserve the tissue structures. Hence, the superpixel is expected to be a more suitable kernel for local statistics computation than the conventional uniform rectangular kernel. Local statistics with the superpixel kernel may enable more accurate quantitative analysis.

Funding

This research was supported in part by the Japan Society for the Promotion of Science (JSPS) (KAKENHI 15K13371 and 15K20887) and the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) through a contract of the Local Innovation Ecosystem Development Program.

Acknowledgments

The authors thank Masahiro Miura from Tokyo Medical University for his support in interpreting the retinal JM-OCT data. The research administrative work of Tomomi Nagasaka from the University of Tsukuba is also gratefully acknowledged.

Disclosures

AM, YJH, DK: Topcon (F), Tomey Corporation (F), Nidek (F). SM, YY Topcon (F), Tomey Corporation (F, P), Nidek (F, P). YJH is currently employed by Koh Young Technology.

References and links

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, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]

2. E. A. Swanson and J. G. Fujimoto, “The ecosystem that powered the translation of OCT from fundamental research to clinical and commercial impact [Invited],” Biomed. Opt. Express 8, 1638–1664 (2017). [CrossRef] [PubMed]

3. C. A. Puliafito, M. R. Hee, C. P. Lin, E. Reichel, J. S. Schuman, J. S. Duker, J. A. Izatt, E. A. Swanson, and J. G. Fujimoto, “Imaging of macular diseases with optical coherence tomography,” Ophthalmology 102, 217–229 (1995). [CrossRef] [PubMed]

4. M. Pircher, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Polarization sensitive optical coherence tomography in the human eye,” Prog. Retin. Eye Res. 30, 431–451 (2011). [CrossRef] [PubMed]

5. J. F. d. Boer, C. K. Hitzenberger, and Y. Yasuno, “Polarization sensitive optical coherence tomography – a review [Invited],” Biomed. Opt. Express 8, 1838–1873 (2017). [CrossRef] [PubMed]

6. S. Guo, J. Zhang, L. Wang, J. S. Nelson, and Z. Chen, “Depth-resolved birefringence and differential optical axis orientation measurements with fiber-based polarization-sensitive optical coherence tomography,” Opt. Lett. 29, 2025–2027 (2004). [CrossRef] [PubMed]

7. S. Makita, M. Yamanari, and Y. Yasuno, “Generalized Jones matrix optical coherence tomography: performance and local birefringence imaging,” Opt. Express 18, 854–876 (2010). [CrossRef] [PubMed]

8. D. Kasaragod, S. Makita, S. Fukuda, S. Beheregaray, T. Oshika, and Y. Yasuno, “Bayesian maximum likelihood estimator of phase retardation for quantitative polarization-sensitive optical coherence tomography,” Opt. Express 22, 16472–16492 (2014). [CrossRef] [PubMed]

9. M. Yamanari, S. Tsuda, T. Kokubun, Y. Shiga, K. Omodaka, Y. Yokoyama, N. Himori, M. Ryu, S. Kunimatsu-Sanuki, H. Takahashi, K. Maruyama, H. Kunikata, and T. Nakazawa, “Fiber-based polarization-sensitive OCT for birefringence imaging of the anterior eye segment,” Biomed. Opt. Express 6, 369–389 (2015). [CrossRef] [PubMed]

10. S. Sugiyama, Y.-J. Hong, D. Kasaragod, S. Makita, S. Uematsu, Y. Ikuno, M. Miura, and Y. Yasuno, “Birefringence imaging of posterior eye by multi-functional Jones matrix optical coherence tomography,” Biomed. Opt. Express 6, 4951–4974 (2015). [CrossRef] [PubMed]

11. W. C. Y. Lo, M. Villiger, A. Golberg, G. F. Broelsch, S. Khan, C. G. Lian, W. G. Austen, M. Yarmush, and B. E. Bouma, “Longitudinal, 3d Imaging of Collagen Remodeling in Murine Hypertrophic Scars In Vivo using Polarization-sensitive Optical Frequency Domain Imaging,” J. Invest. Dermatol. 136, 84–92 (2016). [CrossRef] [PubMed]

12. E. Li, S. Makita, Y.-J. Hong, D. Kasaragod, and Y. Yasuno, “Three-dimensional multi-contrast imaging of in vivo human skin by Jones matrix optical coherence tomography,” Biomed. Opt. Express 8, 1290–1305 (2017). [CrossRef] [PubMed]

13. D. Kasaragod, S. Makita, Y.-J. Hong, and Y. Yasuno, “Noise stochastic corrected maximum a posteriori estimator for birefringence imaging using polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 8, 653–669 (2017). [CrossRef] [PubMed]

14. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography,” Opt. Express 16, 16410–16422 (2008). [CrossRef] [PubMed]

15. C. Ahlers, E. Götzinger, M. Pircher, I. Golbaz, F. Prager, C. Schütze, B. Baumann, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Imaging of the retinal pigment epithelium in age-related macular degeneration using polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 51, 2149–2157 (2010). [CrossRef]

16. S. Makita, Y.-J. Hong, M. Miura, and Y. Yasuno, “Degree of polarization uniformity with high noise immunity using polarization-sensitive optical coherence tomography,” Opt. Lett. 39, 6783–6786 (2014). [CrossRef] [PubMed]

17. B. Baumann, S. Zotter, M. Pircher, E. Götzinger, S. Rauscher, M. Glösmann, J. Lammer, U. Schmidt-Erfurth, M. Gröger, and C. K. Hitzenberger, “Spectral degree of polarization uniformity for polarization-sensitive OCT,” J. Mod. Opt. 62, 1758–1763 (2015). [CrossRef] [PubMed]

18. B. Baumann, S. O. Baumann, T. Konegger, M. Pircher, E. Götzinger, F. Schlanitz, C. Schütze, H. Sattmann, M. Litschauer, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Polarization sensitive optical coherence tomography of melanin provides intrinsic contrast based on depolarization,” Biomed. Opt. Express 3, 1670–1683 (2012). [CrossRef] [PubMed]

19. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14, 7821–7840 (2006). [CrossRef] [PubMed]

20. L. An and R. K. Wang, “In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography,” Opt. Express 16, 11438–11452 (2008). [CrossRef] [PubMed]

21. Y. Jia, O. Tan, J. Tokayer, B. Potsaid, Y. Wang, J. J. Liu, M. F. Kraus, H. Subhash, J. G. Fujimoto, J. Hornegger, and D. Huang, “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express 20, 4710–4725 (2012). [CrossRef] [PubMed]

22. S. Makita, K. Kurokawa, Y.-J. Hong, M. Miura, and Y. Yasuno, “Noise-immune complex correlation for optical coherence angiography based on standard and Jones matrix optical coherence tomography,” Biomed. Opt. Express 7, 1525–1548 (2016). [CrossRef] [PubMed]

23. Y. Lim, Y.-J. Hong, L. Duan, M. Yamanari, and Y. Yasuno, “Passive component based multifunctional Jones matrix swept source optical coherence tomography for Doppler and polarization imaging,” Opt. Lett. 37, 1958–1960 (2012). [CrossRef] [PubMed]

24. M. J. Ju, Y.-J. Hong, S. Makita, Y. Lim, K. Kurokawa, L. Duan, M. Miura, S. Tang, and Y. Yasuno, “Advanced multi-contrast Jones matrix optical coherence tomography for Doppler and polarization sensitive imaging,” Opt. Express 21, 19412–19436 (2013). [CrossRef] [PubMed]

25. S. Kwon, Y. Yoon, B. Kim, W. H. Jang, B. Oh, K. Y. Chung, and K. H. Kim, “Dermoscopy guided dark-field multi-functional optical coherence tomography,” Biomed. Opt. Express 8, 1372–1381 (2017). [CrossRef] [PubMed]

26. X. Ren and J. Malik, “Learning a classification model for segmentation,” Proc. Ninth IEEE Int. Conf. Comp. Vis. pp. 10–17 vol.1 (2003).

27. R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, and S. Süsstrunk, “SLIC Superpixels,” EPFL Technical Report 149300 (2010).

28. R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, and S. Süsstrunk, “SLIC Superpixels Compared to State-of-the-Art Superpixel Methods,” IEEE Trans. Pattern Anal. Mach. Intell. 34, 2274–2282 (2012). [CrossRef] [PubMed]

29. American National Standard Institute, American National Standard for the Safe Use of Lasers ANSI Z 136.1-2014 (American National Standards Institute, New York, 2014).

30. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in Optical Coherence Tomography,” Journal of Biomedical Optics 4, 95–105 (1999). [CrossRef] [PubMed]

31. J. J. Weiter, F. C. Delori, G. L. Wing, and K. A. Fitch, “Retinal pigment epithelial lipofuscin and melanin and choroidal melanin in human eyes,” Investigative Ophthalmology & Visual Science 27, 145–152 (1986).

32. U. Bhaskar, Y.-J. Hong, M. Miura, and Y. Yasuno, “Five-dimensional analysis of multi-contrast Jones matrix tomography of posterior eye,” Proceedings of SPIE 8930, 893008 (2014). [CrossRef]

33. A. Miyazawa, Y.-J. Hong, S. Makita, D. Kasaragod, and Y. Yasuno, “Accurate measurement of retinal birefringence by jones matrix oct with maximum a posteriori estimator and superpixel kernel,” ARVO Imaging in the Eye Conference (2017).

34. E. Götzinger, M. Pircher, B. Baumann, T. Schmoll, H. Sattmann, R. A. Leitgeb, and C. K. Hitzenberger, “Speckle noise reduction in high speed polarization sensitive spectral domain optical coherence tomography,” Opt. Express 19, 14568–14584 (2011). [CrossRef] [PubMed]

35. Y. Lim, M. Yamanari, S. Fukuda, Y. Kaji, T. Kiuchi, M. Miura, T. Oshika, and Y. Yasuno, “Birefringence measurement of cornea and anterior segment by office-based polarization-sensitive optical coherence tomography,” Biomedical Optics Express 2, 2392–2402 (2011). [CrossRef] [PubMed]

36. M. Yamanari, S. Tsuda, T. Kokubun, Y. Shiga, K. Omodaka, N. Aizawa, Y. Yokoyama, N. Himori, S. Kunimatsu-Sanuki, K. Maruyama, H. Kunikata, and T. Nakazawa, “Estimation of Jones matrix, birefringence and entropy using Cloude-Pottier decomposition in polarization-sensitive optical coherence tomography,” Biomedical Optics Express 7, 3551 (2016). [CrossRef] [PubMed]

	Training B-scan 1	Test-1 B-scan 2	Test-2 B-scan 3	Test-3 B-scan 4
Variance of contributions	0.000111	0.00286	0.00377	0.00523
Largest normalized-ISPV with optimal weights	0.551	0.569	0.578	0.572
Largest normalized-ISPV without optimization	0.684	0.684	0.673	0.772

	Training B-scan 1	Test-1 B-scan 2	Test-2 B-scan 3	Test-3 B-scan 4
Variance of contributions	0.0000700	0.00206	0.00221	0.00404
Largest nISPV with optimal weights	0.553	0.589	0.581	0.593
Largest nISPV without optimization	0.648	0.662	0.655	0.668

	Training-1 Subject-1	Test-1 Subject-1	Training-2 Subject-2	Test-2 Subject-2	Training-3 Subject-3	Test-3 Subject-3
VOC	0.000111	0.00244	0.0000304	0.00609	0.0000787	0.00446
Largest nISPV, optimal	0.551	0.592	0.567	0.590	0.541	0.586
Largest nISPV, w/o optimization	0.684	0.750	0.731	0.777	0.795	0.781

	Training-1 Subject-2	Test-1 Subject-2	Training-2 Subject-4	Test-2 Subject-4	Training-3 Subject-1	Test-3 Subject-1
VOC	0.000111	0.00541	0.00000604	0.00307	0.00000785	0.00944
Largest nISPV, optimal	0.553	0.604	0.552	0.582	0.546	0.609
Largest nISPV, w/o optimization	0.648	0.676	0.641	0.656	0.630	0.627

	Training Subject-1	Test-1 Subject-3	Test-2 Subject-4	Test-3 Subject-2
VOC	0.000111	0.00368	0.00373	0.00516
Largest nISPV, optimal	0.551	0.603	0.571	0.587
Largest nISPV, w/o optimization	0.684	0.731	0.795	0.746

Generation and optimization of superpixels as image processing kernels for Jones matrix optical coherence tomography

Abstract

1. Introduction

2. JM-OCT system and measurement protocol

3. Generation and optimization of superpixels for JM-OCT

3.1. Superpixel generation

3.1.1. Generation of initial superpixels

3.1.2. Definition of distance in feature space

3.1.3. Superpixel reshaping

3.2. Optimization method of superpixel parameters

3.2.1. Intra-superpixel variance and contribution metric

3.2.2. Optimization of weights for the optical features

3.2.3. Optimization strategy of the compactness factor

3.2.4. Flow of optimizations

4. Results

4.1. Validation of the superpixel parameter optimization methods

4.1.1. Validation of feature weights optimization

4.1.2. Validation of the compactness factor optimization

4.2. Superpixelization of in vivo posterior eye image

4.3. Conventional rectangle kernel and superpixel kernel

4.4. Superpixels generated only form OCT intensity

5. Discussions

5.1. Speckle effect on superpixels

5.2. Optimal definition of spatial distance

5.3. Universality of the optimized weights of the optical features

5.4. Computational time

5.5. Image segmentation and superpixel

5.6. Re-estimation of optical properties by superpixel kernels

6. Conclusion

Funding

Acknowledgments

Disclosures

References and links

Cited By

Figures (9)

Tables (6)

Equations (14)

Biomedical Optics Express

	Training Subject-2	Test-1 Subject-3	Test-2 Subject-4	Test-3 Subject-1
VOC	0.0000700	0.0116	0.0165	0.00602
Largest nISPV, optimal	0.553	0.627	0.614	0.567
Largest nISPV, w/o optimization	0.648	0.731	0.641	0.630