Theoretical model for en face optical coherence tomography imaging and its application to volumetric differential contrast imaging

A new formulation of the lateral imaging process of point-scanning optical coherence tomography (OCT) and a new differential contrast method designed by using this formulation are presented. The formulation is based on a mathematical sample model called the dispersed scatterer model (DSM), in which the sample is represented as a material with a spatially slowly varying refractive index and randomly distributed scatterers embedded in the material. It is shown that the formulation represents a meaningful OCT image and speckle as two independent mathematical quantities. The new differential contrast method is based on complex signal processing of OCT images, and the physical and numerical imaging processes of this method are jointly formulated using the same theoretical strategy as in the case of OCT. The formula shows that the method provides a spatially differential image of the sample structure. This differential imaging method is validated by measuring in vivo and in vitro samples.


Introduction
Several theoretical models of optical coherence tomography (OCT) [1] have been studied to understand its image formation. Because OCT is used in the biomedical field such as ophthalmology [2] and dermatology [3], theoretical modeling is important for interpreting OCT images and clarifying what is visualized in the sample tissue. Theoretical modeling is also useful for designing optimal OCT systems and developing image processing algorithms.
Several analytic and numerical models of OCT imaging have been developed. After an early trial of OCT modeling [4], the extended Huygens-Fresnel formulation was introduced [5,6]. In this model, propagation of the OCT probe beam in a tissue, i.e., an inhomogeneous medium, is modeled by extended Huygens-Fresnel formulation. Numerical simulations based on this model can recapitulate the forward multiple scatting in the tissue. Monte-Carlo methods were also widely used to numerically emulate the OCT imaging [7,8]. Although Monte-Carlo methods are nondeterministic and computationally expensive in general, they can be applied to wide varieties of imaging models, conditions, and samples. The interaction between the probe beam and the sample can be also accurately computed by using Maxwell's equation [9,10]. An extension of this method, the pseudo-spectral time domain method, enables the realistic simulation of three-dimensional OCT imaging [11]. Recently, the speckle was mathematically described as the interference of light backscattered by multiple scatterers in a resolution volume [12].
Aforementioned approaches aimed to recapture accurate and realistic OCT signals. By contrast, we aim to establish a simplified theoretical framework that enables the intuitive understanding of the point-scanning OCT imaging process and the interpretation of an en face OCT image. In this framework, we model a sample as a spatially slowly varying refractive index in which scatterers are randomly dispersed (Section 2). The density of the scatterers is also spatially slowly varying. This representation of the sample is denoted as a "dispersed scatterer model" (DSM). We formulate the OCT imaging process using the DSM (Section 3). This formulation identifies the relationship between the sample structure, that is, the slowly varying refractive index distribution and the scatterer density distribution, and the en face OCT signal. In this formulation, a meaningful OCT signal and speckle appear as different terms; hence, it provides a comprehensive understanding of the en face OCT image formation process.
By actively using this theoretical framework, we design and experimentally demonstrate an OCT-based volumetric differential contrast (VDC) imaging method (Section 4). We implement this method as the numerical signal processing flow of complex en face OCT images and provide differential contrast images similar to those of differential interference contrast microscopy (DIC) [13]. Similar to DIC, VDC highlights the micro-structure of a sample, but unlike DIC, VDC can be volumetric.
The DSM-based framework can also explain phenomena that occur in OCT-based imaging methods but have not been theoretically well explained. An example is artifacts found in A-scanwise Doppler OCT images caused by the structure of superior tissue. Another example relates to the specific imaging condition required for point-scanning computational multi-directional OCT [14]. In experiments, it was found that defocus is mandatory to achieve this method [15], but this necessity has not been explained in theory. We analyze these issues using the DSM-based framework in Section 5.

Dispersed scatterer model (DSM)
In our theoretical framework, we model a sample to be measured as a spatially slowly varying refractive index distribution in which scatterers are spatially randomly dispersed. We consider the density of the scatterers to be also spatially slowly varying. We refer to this theoretical model of the sample as a DSM. Figure 1 schematically illustrates the DSM, where (x, y, z) is the position in the sample as x and y are the lateral positions, and z is the depth position. The spatially slowly varying refractive index is represented by n(x, y, z), in which small scatterers are dispersed (i.e., embedded), where the scatterers are assumed to be significantly smaller than the spatial variation of n(x, y, z). The position of the i-th scatterer is (x i , y i , z i ) and all the scatterers have an identical refractive index of n s . The probe beam is incident from above of the sample, and the green region in the figure represents the coherence gate. In our formulations for OCT (Section 3) and VDC (Section 4), we consider en face images in which scattering only within this coherence gate contributes to imaging. The sample is represented as a spatially slowly varying refractive index n(x, y, z) in which infinitely small scatterers (black dots) are randomly dispersed. All the scatterers are assumed to have the same refractive index of n s . The green region represents the coherence gate of OCT and only the scatterers within the coherence gate contribute to OCT imaging. The probe beam is incident on the sample from the top and only lights backscattered toward the incident direction contribute to imaging.

Formulation of OCT imaging
In this section, we formulate the image formation process of standard scanning OCT imaging using DSM.

Complex OCT signal
We assume that a Gaussian beam is incident on the sample. Considering that N scatterers are involved in imaging, i.e., N scatterers exist in the coherence gate, a complex en face OCT signal becomes where x = (x, y) and x ′ = (x ′ , y ′ ) are the two-dimensional lateral positions on a sample and image, respectively. Here the contribution from each scatterer is that of single scattering. Namely, the first-order Born approximation has been employed. Throughout this paper, we use x and z to represent the position in the object (sample) space, and x ′ and z ′ to represent the position in the image space. z ′ 0 is the depth position in the image corresponding to the center depth of coherence gate z 0 in the sample. The magnification between the sample and the image is assumed to be 1.
is the lateral position of the i−th scatterer. B(x, z 0 ) is a quantity defined by the slowly varying refractive index distribution of the medium, n(x, y, z), and the refractive index of the embedded scatterers, n s . And hence, B is also spatially slowly varying as is n. It is noteworthy that the overall profile of B(x, z 0 ) is defined solely by the profile of the spatially slowly varying refractive index distribution n(x, y, z). The spatial inhomogeneity of the dispersed scatterers is accounted for later when we derive the continuous form of the imaging formula in Section 3.3. For the case of sufficiently small scatterers, B(x, z 0 ) can be considered as the scattering potential where λ 0 is the central wavelength of the probe beam [16].
︁ is a complex point spread function (PSF) at depth z ′ 0 . ϕ i is the phase offset caused by the depth position of the i-th scatterer with respect to the central depth of the coherence gate z 0 : where k 0 = 2π/λ 0 is the wavenumber as λ 0 is the center wavelength of the probe beam and n is the refractive index of the sample around the scatterers. Note that we assume that B is homogeneous along the depth within the coherence gate; hence, the variation of the scatterer's depth position does not affect B but only appears as phase offset.
where P a ≡ |P c | is the amplitude of the PSF and is typically Gaussian in the lateral direction, and φ is the phase of the complex PSF. We discuss a more practical shape for the complex PSF in Section 3.4. By substituting Eq.
As shown in Eq. (4), the contribution from each scatterer is weighted by the amplitude of the complex PSF, P a , and its phase is modulated by both the depth position of the scatterer and spatial phase pattern of the PSF.

Intensity-OCT signal in discrete form
The intensity of OCT is the absolute square of the complex OCT signal, that is, the product of the complex OCT signal and its conjugate. From Eq. (4), the en face OCT intensity I (x ′ ; z 0 ) becomes where x j = (x j , y j ) is the position of the j-th scatterer in the sample. On the last line of the equation, we did not explicitly write z 0 and z ′ 0 for simplicity. We also omitted them in the later part of the paper if it did not cause confusion. Eq. (5) is schematically depicted in Fig. 2(a). Each circle represents each term of the equation, that is, the contribution from each scatterer in the sample. The same color indicates the same term, that is, the contribution from the same scatterer, and the asterisk represents the complex conjugate. The red box consists of the terms in s   Physical depiction of the numerical interactions between two signals contributed by each pair of scatterers, which are mathematically described by Eq. (6). The interactions are classified into two groups: auto-interaction (AI) and cross-interaction (CI) groups. AI corresponds to the numerical interference of the electric field scattered from a scatterer and the conjugate field from the same scatterer; hence, it corresponds to the backscattered intensity from the scatterer. CI corresponds to the numerical interference of the electric field from a scatterer and the conjugate field from another scatterer; hence, is a phasor. AI and CI are also referred to as auto-interference and cross-interference, respectively.
The terms of Eq. (5) can be classified into two groups for i = j and i ≠ j as This equation is the representation of the OCT intensity signal in the discrete form of our theoretical framework. In this formulation, each term in the first group (the first line of the equation, for i = j, written in blue) is the numerical interaction between contributions from the same scatterer; hence, we refer to this group as an auto-interaction (AI) group. By contrast, each term in the second group (the second line, for i ≠ j, written in brown and red) is the numerical interaction between contributions from different scatterers. We refer to this group as a cross-interaction (CI) group. We classified CI into brown and red parts, and this classification is used in later sections. The terms of the AI group correspond to the diagonal terms in the orange box in Fig. 2(a), whereas those of the CI group are the off-diagonal terms.
It is noteworthy that these interactions can be considered as the numerical interferences of scattered electric fields from the scatterers as depicted in Fig. 2(b). The AI term corresponds to interference between an electric field from a scatterer and the complex conjugate of the same electric field. Hence, each term in the AI group represents the squared intensity of scattered light from each scatterer and the AI group is the incoherence summation of these scattered light intensities. By contrast, each term in the CI group is the numerical interference between the scattered complex field from a scatterer and the conjugate complex field from another scatterer, where two scatterers are sufficiently close to be included in the complex PSF. Note that each term in the CI group is a phasor with a practically random phase and each term in the CI group has its complex conjugate pair, as can be schematically seen from the diagonal symmetricity of the orange box in Fig. 2(a). Because of the random phase of each term, the summation of a term and its complex conjugate becomes a random real value. This suggests that the CI group causes a random signal, i.e., speckle, as we discuss in detail in Section 3.5.

Intensity-OCT signal in continuous form
In Eq. (6), the OCT intensity is expressed using the contributions from discrete scatterers; hence, the equation is also in discrete form. In this subsection, we derive a continuous-form representation of the OCT intensity, which is a convolution-based conventional image formation equation, and hence convenient to interpret.
Using the Dirac delta function, the intensity-OCT signal in discrete form [Eq. (6)] can be rewritten as where x p and x q are variables in real physical space in the sample. Each of the double integral is for x and y components of the vectorized integration variable x p , while the quadruple integral is for two vectorized integration variables x p and x q . The intervals of all the integrals are from −∞ to +∞. Hereafter, the integral intervals of all equations are from −∞ to +∞; we omit them for simplicity. ϕ pq is a phase that is originated from ϕ i − ϕ j in Eq. (6). Namely, ϕ pq is the phase difference between the contributions of i-th and j-th scatterer and equals to ϕ i − ϕ j . Similar to ϕ i − ϕ j , ϕ pq is unpredictable and random in practice, and ϕ pq = −ϕ qp . Note that, in the second line of Eq. (7), we can consider that x p and x q respectively denote the positions of i-th and j-th scatterers due to the Dirac delta function. The first line of Eq. (7) (blue) corresponds to the AI group [the first line of Eq. (6)] and the second to fourth lines of Eq. (7) (brown and red) correspond to the CI group [the second line of Eq. (6)]. Note that the summation in CI of Eq. (6) is for i ≠ j, whereas we removed this condition in Eq. (7). Alternatively, we subtract the fourth line, which is identical to the first line corresponding to the case of i = j. This removal of condition was necessary to convert the discrete equation [Eq. (6)] into the continuous equation [Eq. (7)] because the integral cannot properly account for such a condition. We introduce functions D and D ′ as which are two-and four-dimensional binary maps of the scatterers' positions, respectively. When the scatterer density is high enough to be considered as continuously distributed, that is, N is sufficiently large, D and D ′ can be considered as non-binary scatterer density maps in two-and four-dimensional spaces, respectively. Specifically, D is the scatterer density map of the sample. By substituting Eqs. (8) and (9) into Eq. (7), we obtain the continuous form of the equation: The first integral (i.e., the AI part) of Eq. (10) (blue) is "an imaging equation of the intensity-OCT image" in continuous form. According to this equation, the object to be imaged in the intensity-OCT image is the product of the squared scattering potential B 2 and the scatterer density D, and the PSF of intensity-OCT imaging is P 2 a . A more detailed interpretation of this intensity-OCT imaging formula is provided in the following subsections (Sections 3.4 and 3.5).

Intensity-OCT imaging with Gaussian PSF
The property of intensity-OCT imaging becomes clearer if we assume a complex Gaussian PSF with defocus in the continuous OCT imaging formula [Eq. (10)]. The complex Gaussian PSF is obtained from a Gaussian probe beam, and the Gaussian probe beam is typical for a single-mode-fiber-based point-scanning OCT. According to the formulation by Ralston et al. [17], the Gaussian probe beam in the paraxial approximation becomes where z 0 is the depth position of the focus, w(z; z 0 ) is the radius of the probe beam at depth z, w 0 = w(z = z 0 ; z 0 ) is the beam radius at the in-focus depth, that is, the beam waist radius, and n is the refractive index of the sample around the position of interest. The second term in the second exponential represents the quadratic phase induced by defocus and R(z; z 0 ) is the phase curvature.
Here after the squared vector represents the inner product of a vector with itself. In the third term of the second exponential, ψ(z; z 0 ) represents the Gouy phase. The beam waist radius w 0 , the radius of the probe beam w, and the phase curvature R are where f is the focal length of the objective and Φ is the 1/e 2 -diameter of the probe beam incident on the objective.
Because OCT is a reflective imaging modality, and the illumination and collection paths share the same optics, the complex PSF P c (x, z; z 0 ) is obtained by multiplying G(x, z; z 0 ) by itself as On the first line, each G(x, z; z 0 ) is for probe beam illumination and collection, respectively. By comparing this complex PSF and that of Eq. (3), we find that and By substituting this Gaussian complex PSF, or equally P a and φ, into Eq. (6), we obtain the discrete-form intensity-OCT imaging equation: We omit depth positions for simplicity. The first line of this equation (blue) corresponds the AI group, while second and third lines (brown and red) correspond to the CI group. Using the symmetricity of x i and x j , the second summation part can be rewritten in a real form from Eq. (18) as By applying the same manipulations used to convert Eq. (6) to Eq. (10), which we described in Section 3.3, we obtain the continuous form of the intensity-OCT equation with the Gaussian PSF with defocus from Eq. (18) as where R.H.S. means "the right hand side." The last line, which is identical to the AI group (blue), was subtracted from the second integral, so that the last three lines of the equation (brown and red) represent the CI group. This equation is the continuous-form imaging equation of an intensity-OCT image with a complex Gaussian PSF with defocus. The first line and other lines correspond to AI and CI, respectively. Using the symmetricity of x p and x q , and ϕ pq = −ϕ qp , the second integral of Eq. (20) is rewritten in a real form as
All the equations suggest that the OCT intensity consists of AI [i.e., the first summation or integral of each equation (blue) of Eqs. (6), (10), (19), and (21)], and CI (i.e., the residual terms, brown and red). AI forms a meaningful image, as we discuss in this paragraph, whereas CI corresponds to speckle, which we discuss in the next paragraph. As can be seen in the equations, for the intensity-OCT image formed by AI, the object to be measured is B 2 (x i ) for discrete forms and D(x)B 2 (x) for continuous forms. Additionally, it is noteworthy that this intensity-OCT imaging explicitly has a real-and positive-valued PSF: P a 2 . We found that the AIs of the continuous-form equations [Eqs. (10) and (21)] are the convolution of the real-positive PSF and the object to be measured. And hence, the AI is equivalent to the incoherent image of the object.
In contrast to AI, CI does not have an explicit PSF. To intuitively understand the image property of CI, we split the product parts inside the summation or integral into two parts: the object to be measured (B 2 (x i ) or D(x)B 2 (x)) and the residual part (written in red). We refer to the residual part as "pseudo-PSF" because it can be considered as the PSF for each scatterer pair. The pseudo-PSF varies for each combination of scatterers; hence, CI does not form a meaningful image. Namely, CI represents speckle. Specifically for the intensity-OCT imaging under a Gaussian complex PSF [Eqs. (19) and (21)], the pseudo-PSF is found to have the following characteristics. First, the pseudo-PSF is the product of the Gaussian envelope (the second exponential parts in CI) and a cosine wave. Second, the Gaussian envelope has the same width as the explicit PSF of AI (the exponential part of AI). It is consistent with a well-known fact; speckle size is a function of the OCT resolution [5,18]. Third, the center of the Gaussian is located at the center of the two scatterers of each scatterer pair; that is, the center of the Gaussian envelope walks around for each scatterer pair and is practically unpredictable and random. Fourth, the frequency and phase offset, that is, the shift of the cosine wave, depends on the positions of the scatterers of the scatterer pair; hence, they are practically unpredictable and random.
These formulations indicate that the intensity-OCT image is a summation of a speckle-free meaningful image, i.e., the incoherent image of the object, which corresponds to AI, and a random speckle image corresponding to CI.

Theory
In this section, we present a new OCT-based imaging method called volumetric differential contrast (VDC) imaging, which we design using the DSM and a similar theoretical framework presented in Section 3. This is an image processing method that manipulates the en face OCT signal and can provide a spatially differentiated image of a sample that highlights the sample structure as contrast, similar to DIC. In contrast to standard DIC, VDC is volumetric, that is, it provides a spatially differentiated image at each depth of a sample.
An en face VDC image is obtained by multiplying an en face complex OCT by its digitally shifted complex conjugate and then extracting the imaginary part. We describe the formulation of VDC in this section.
In the first step of an en face VDC computation, we compute the product of an en face complex OCT signal and the digitally laterally shifted complex conjugate of this signal. Using the same theoretical flow that we used to obtain Eq. (10), this product, which is an en face complex signal, is written as where ∆x is the digital shift amount between the two signals, typically one or two lateral pixels, and should be smaller than the spot size of the probe beam. The letter colors corresponds to those in the previous equations. The detailed derivation of this equation is in the Appendix.
By assuming the Gaussian complex PSF described in Section 3.4, that is, by substituting Eqs. (16) and (17) into Eq. (22), we obtain the en face complex signal: Similar to the conventional intensity-OCT image [Eq. (10)], this en face signal consists of the AI part [the first two lines of Eq. (23), blue] and the CI part (the next four lines, brown and red). The AI part forms an imaging formula in which the object to be measured is the same as the conventional intensity-OCT, that is, D(x)B 2 (x), the product of the scatterer density and the squared scattering potential. The PSF of this imaging formula is a complex PSF (the second line of the equation), which is the product of a Gaussian and a lateral linear phase slope. By contrast, the CI part is practically random, which is, again, similar to the conventional intensity-OCT. Now we define the VDC image as the imaginary part of Eq. (23) as In this equation, the phase parts of Eq. (23), i.e., the last exponential part of each integral, became sinusoidal function. The last four lines of Eq. (24) (brown and red) correspond to the CI part. Similar to the case of the intensity-OCT formula, ϕ pq and pseudo-PSF (red) are practically random; hence, CI becomes a random speckle pattern. By omitting this random CI part, we can define the VDC image [VDC(x ′ )] and write it as the first integral of Eq. (24): We interpret this imaging formula as the convolution of the object D(x)B 2 (x) with an imaging PSF that is a Gaussian-windowed sine function. This imaging PSF function is similar to the response function (i.e., the imaging PSF) of DIC [19] and also a differential operator used in digital image processing [20]. Hence, it works as a differential operator in VDC imaging and provides the differential contrast of the sample structure. The shape of the imaging PSF of VDC (i.e., the differential operator) is altered by the defocus amount z − z 0 and shift ∆x. Figure 3 exemplifies several PSFs with different shifts (a) and defocus amounts (b). We computed PSFs using OCT system parameters identical to those in the experimental setup described in Section 4.2: w 0 = 2.41 µm, n = 1.38, and k 0 = 7.48 × 10 6 rad/m. The shifting amounts were the integer multiples of the lateral pixel size (1.95 µm) and the defocus amounts were expressed as the distance from the focal depth, that is, z − z 0 in Eqs. (13) and (14). For Fig. 3(a), we set the defocus amount to 200 µm, whereas for Fig. 3(b), we set the shifting amount to 1.95 µm. The one-dimensional plots are the profiles at the center of the PSF, i.e., x ′ = x, and are in the digital shifting direction. The columns "Gaussian part" and "sine part" in the figure represent the Gaussian and sine parts within the integral of Eq. (25), respectively.
As shown in Fig. 3(a), the shift affected the frequency of the sine function in the imaging PSF and did not affect the Gaussian function. By contrast, the defocus amount affected both the Gaussian and sine functions of the imaging PSF [ Fig. 3(b)]. Please note that the absolute defocus amount was larger than the Rayleigh length (21.7 µm) in the examples in the figure.
By properly selecting the shift and defocus amounts, we can obtain a good differential operator with an arbitrary differential distance, which is defined as the distance between two main positive and negative peaks. The differential distances are shown at the top right of the 2D-PSF images. (By contrast, an example of the "inappropriate" selection of these amounts is a shift that is too large, which fluctuates the operator within the Gaussian and an appropriate differential contrast cannot be obtained.) When the defocus direction is reversed, that is, the plus/minus sign of the defocus is changed, the sine function and differential direction are flipped, as shown in Fig. 3(b). When the shift is zero, the sine function becomes zero; hence, VDC(x ′ ) = 0; that is, differential contrast is not obtained. Similarly, when the defocus amount is zero, that is, z = z 0 , R becomes infinity and the sine part of the operator becomes sin 0 = 0; hence, contrast is not obtained again. Thus, not only the shift but also the defocus is mandatory for VDC.
In our particular implementation, we computationally apply defocusing as detailed in Section 4.3.3; hence, both the shift and defocus amount can be selected after image acquisition. It is also noteworthy that the size and direction of the differential operator are analogous to the shear amount and the shear direction of DIC imaging. Because the shift and defocus amount can be set after image acquisition, VDC can generate multiple DIC-like images with arbitrary shear amounts and directions using post-acquisition signal processing.

System, samples, and image acquisition protocol
For the experimental validation of VDC imaging, we used a fiber-based spectral-domain OCT system [14]. The light source was a superluminescent diode (SLD, Superlum, Ireland) with a center wavelength and bandwidth of 840 nm and 100 nm, respectively. In the sample arm, the probe beam with a diameter of 4.0 mm passed through an objective (LSM02BB, Thorlabs, NJ) with a focal length of 18 mm and incident on the sample. The axial resolution was 3.8 µm in tissue and the in-focus lateral resolution was 4.9 µm. The spectral interference signal was acquired using a spectrometer (Cobra S800, Wasatch Photonics, NC). The scanning speed was 50,000 lines/s with 2,048 spectral pixels.
To reconstruct the OCT volume, spectral rescaling, numerical dispersion compensation, and fixed-pattern noise removal were applied. The depth-independent bulk phase error caused by system imperfection and air turbulence was numerically corrected after the reconstruction of  The larger the defocus amount, the wider the Gaussian width and the higher the frequency of the sine part, and vice versa, provided that the defocus amount is larger than the Rayleigh length. When the defocus amount is zero, the sine part becomes zero again; hence, the image is not formed. The appropriate selection of the shifting and defocus amounts leads to the PSF becoming a differential operator. The direction of differentiation can be flipped by flipping the shift and defocus directions. Note that the imaging PSF of VDC is analogous to the response function of DIC, and the shear amount and direction of DIC are analogous to the size and direction of the imaging PSF of VDC.
the OCT volume using the smart-integration-path method [21]. Note that, since our image processing exploits complex signals, high phase stability is required An 8-day-old zebrafish larva and a human breast cancer (MCF7) spheroid were measured. The zebrafish was anesthetized and embedded in agarose gel. The lateral scanning areas were 1 mm × 1 mm for the zebrafish and 0.8 mm × 0.8 mm for the spheroid, with 512 × 512 A-lines. Figure 4 shows

Processing of VDC images
The VDC processing flow is summarized in Fig. 5. (1) An en face plane was extracted from a complex OCT volume. (2) The en face signal was first refocused using computational refocusing based on a forward-propagation model [21,22] and then arbitrary numerical defocusing was applied (similar to that in Section 5.1.1 of Ref. [23]). (3) This complex en face OCT signal with arbitrary defocusing, s 1 (x ′ ), was copied and digitally shifted in the lateral direction to obtain s 2 (x ′ ). (4,5) Then the product of s 1 (x ′ ) and the complex conjugate of s 2 (x ′ ) was computed. (6) Finally, the en face VDC image was obtained as the imaginary part of the product.

Shifting direction dependency of VDC imaging
The En face VDC images of the (a) zebrafish and (b) spheroid are summarized in Fig. 6, which demonstrates the multi-directional differential imaging capability of VDC. The plots on the right show some representative line profiles of the VDC images. The positions of the line profiles are indicated by long yellow arrows on the VDC images. The shifting amounts ∆x = (∆x, ∆y) were the combination of -1, 0, and +1 pixels (i.e., -1.95, 0.0, and +1.95 µm for the zebrafish and -1.56, 0.0, and +1.56 µm for the spheroid) for the x and y directions. These shift amounts were selected to be the integer multiples of the pixel size to avoid sub-pixel interpolation of the complex OCT image. The defocus amount was +200 µm for all images. Note that this defocus was not a physical defocus but a numerically induced defocus by using complex signal processing as mentioned in Section 4.2.2. The combination of the shift and defocus amounts were selected to make the imaging PSF as a good differential operator. The small insets represent the imaging PSF.
As predicted by the imaging PSF shape of VDC, these images looked like the DIC image. In the VDC images, the edge structures normal to the differential direction (shift direction) were particularly emphasized. As expected, the VDC image was not obtained when there was no shift. Note that, because DIC cannot be performed on a thick sample, such as the zebrafish and spheroid, a direct comparison of VIC and DIC cannot be performed.

Defocus amount dependency of VDC imaging
The effect of the defocus amount is demonstrated in Fig. 8, where the defocus amounts were -200, -100, 0, +100, and +200 µm. The shifting amount was (+1, 0) pixel for the (x, y)-direction for all images. The first and second rows show the VDC images of the zebrafish and spheroid,     Fig. 7. En face VDC images of the (a)-(e) zebrafish and (f)-(j) spheroid with shifting amounts of -2, -1, 0, +1, and +2 pixels. The small insets show the corresponding imaging PSFs and the values above them show the differential distances. ±1-pix shift yielded a larger differential distance than ±2-pix shift, and hence emphasized larger structures than ±2-pixel-shift images.
respectively. The insets show the corresponding imaging PSFs and differential distances are shown above.
Increasing the absolute defocus amount increased the differential distance. When the defocusing direction was reversed, the differential direction was also reversed. The VDC contrast became very low for the in-focus (0-µm defocus) and [ Fig. 8(c) and (e)]. It should be noted that, in theory, the VDC is expected to vanish at in-focus, but low contrast was visible in the images. This residual contrast might be caused by the imperfection of numerical refocusing.   µm. The small insets show the corresponding imaging PSFs and the values above them show the differential distances. The differential distance increased as the amount of defocus increased. The differential direction reversed when the defocus direction reversed.

Volumetric differential imaging
As shown above, VDC yielded DIC-like differential images. However, unlike DIC, VDC achieved volumetric differential imaging by processing each en face plane in an OCT volume. Figure 9 shows examples of the zebrafish (a1-g1) and tumor spheroid (a2-g2). (a1)-(f1) and (a2)-(f2) are the en face VDC images extracted from the VDC volumes (from the upper side to the bottom side). The depth locations of the en face images are indicated in the OCT cross-sections (g1) and (g2). The defocus and shifting amounts were +200 µm and +1 pix along the x-axis for all images. The fly-through en face images for all depths are available as Visualization 1 (zebrafish) and Visualization 2 (tumor spheroid).

Advantage of VDC imaging
As shown in Section 4, VDC imaging provides images similar to DIC, but VDC imaging has two advantages over DIC. First, VDC is volumetric and can image thick tissue. Second, VDC can generate images with several differential conditions with a single acquisition; that is, the differential direction and differential distance, which correspond to the shear direction and amount of DIC, respectively, can be controlled by the parameters of post-acquisition signal processing. This means that VDC can be regarded as a volumetric DIC with post-hoc control of the shear direction and amount.

Difference between VDC and DIC
Despite their image similarity, the VDC and DIC have some fundamental differences in their contrast mechanisms. DIC's contrast is mainly from the optical path length difference between two adjacent lateral points. In addition, DIC is a transillumination microscope, and the sample is thin, low scattering, and typically, transparent. And hence, the DIC's contrast is mainly based on the refractive index distribution of the sample.
On the other hand, VDC is obtained by OCT, which is a reflection (back-scattering) measurement modality, and the sample is assumed to exhibit high scattering. And hence, VDC's contrast might be mainly based on the back-scattering.
One of the not-well-described effects in the current VDC theory is the effect of superior tissues. VDC may be sensitive to the refractive index of the superior tissues, because it alters the phase of the complex en face OCT signal at the depth of interest. However, the present theory only considers an image at a particular depth, and does not formulate the cumulative effect of the superior tissues. Extension of VDC theory to account for the cumulative effect is one of the future work.
In addition, three-dimensional imaging PSF of VDC has not yet been investigated. As demonstrated by Villiger et al. [24] and Ralston et al. [25], the three-dimensional imaging property of OCT can be formulated. The extension of VDC imaging theory into three-dimensional space is also one of the future works.
Well-characterized-phantom based experimental validation might be helpful to further understand the imaging property of VDC. However, for this validation, an inhomogeneous scattering phantom that has a slowly varying refractive index distribution and multiple domains with different scatterer densities are necessary. The establishment of such characterized-phantoms and further experimental validation of VDC might be one of the future studies.

Intensity based differential image generation
It is noteworthy that, as far as we consider only the AI part, the differential contrast image can be obtained by one of the standard intensity-based image processing methods. Figure 10 shows the differential contrast images obtained by the intensity-based image processing corresponding to Fig. 6. Here the images were obtained by convolving the linear-scale OCT intensity image with the imaging PSF of Fig. 6, and we denote these images as "intensity differential images." Two differences are noted between VDC (Fig. 6) and the intensity differential images (Fig. 10). First, the resolution of this intensity differential image might be slightly worse than that of VDC. It is because the original intensity-OCT image used for the intensity differential image has already been affected by the intensity PSF [P 2 a (x) of Eq. (10)]. However, the imaging PSF of VDC is much larger than the intensity PSF due to the numerical defocusing. The defocus amount was 200 µm in this particular case, while the Rayleigh length is 21.7 µm. And hence, this reduction of the resolution might be negligible. The second and more significant difference is that the intensity differential images exhibit more granular appearance than the VDC images. This difference can be explained by the effect of the CI part. In the intensity differential images, CI part, i.e., speckle is also differentiated by the same differentiation kernel with the meaningful AI part. And hence, "differential speckle image" is overlaid on the meaningful differential images. Since both of them have the same resolution, it is hard to distinguish the differential speckle from the meaningful differential image.
On the other hand, the CI part of the VDC [i.e., the last two lines of Eq. (24)] does not have an evident PSF. It has only a pseudo-PSF [see the paragraph just after Eq. (24)], and the width and the magnitude of pseudo-PSF alters for all combinations of the contributing scatterers. And hence, the speckle-originated signal (CI part) appears quite differently from the meaningful differential image (AI part) in VDC. Quantitative comparison of visibility of specific image features in specific applications between VDC intensity-based image-processing methods might be one of the future works in the transnational research.

Artifacts in Doppler OCT
Doppler OCT enables flow imaging in tissue [26]. Among several variations of Doppler OCT, conventional versions compute the phase difference between adjacent A-lines in a B-scan to obtain a Doppler OCT image [27,28]. This version of Doppler OCT is known to exhibit a phase noise pattern caused by the structure of the sample, that is, the structural phase artifact. This structural phase artifact can be mathematically described by extending the formulation of VDC as follows: This version of Doppler OCT can be written by modifying Eq. (23) to s(x ′ + ∆x/2) s * (x ′ − ∆x/2) e iϕ D , where the newly added phase φ D is the Doppler phase shift and ∆x is the separation of the adjacent A-lines. The AI part of this modified equation becomes and the Doppler OCT signal is defined as the phase of this equation. This formulation raises two points about the structural phase artifact. First, if the static structures of the adjacent A-lines are identical, i.e., s(x ′ + ∆x/2) = s(x ′ − ∆x/2), the Doppler OCT signal becomes the Doppler phase shift as Namely, no structural phase artifact appears if there is no structure. Second, if there is no defocus, i.e., R → ∞ in Eq. (26), the phase of Eq. (26) (the Doppler OCT signal) becomes the Doppler phase; otherwise, the Doppler OCT signal is not identical to the Doppler phase. It is noteworthy that, for the second point, we did not assume s(x ′ + ∆x/2) = s(x ′ − ∆x/2), i.e., the structural phase artifact does not appear, even if there is a structure as far as there is no defocus. To summarize, the structural phase artifact appears only when the structure and defocus coexist.
Another structure-related issue is known for Doppler OCT; that is, "structural decorrelation," which means that the measurement accuracy of the Doppler phase degrades as the distance between the adjacent A-lines (∆x) increases. This effect can be related to the first exponential part of Eq. (26); that is, the AI signal decreases as ∆x increases and leads to a low signal-to-noise ratio, and hence, low phase measurement accuracy.

Defocusing effect of computational multi-directional (CMD)-OCT
Oida et al. demonstrated point-scanning computational multi-directional (CMD)-OCT for qualitatively visualizing the microscopic phase structure of the sample [14]. This method exploits the phase angle of s(x ′ + ∆x/2)s * (x ′ − ∆x/2), and hence is similar to our VDC method.
Although it was empirically known that the contrast of CMD-OCT becomes very low if the image is in focus [15], its theoretical basis was not known. This low contrast has now been explained by Eq. (23) of the VDC formulation; that is, if there is no defocus (i.e., R → ∞), the phase of the AI part (meaningful part) of s(x ′ + ∆x/2)s * (x ′ − ∆x/2) becomes zero. This fact explains the contrast vanishing of CMD-OCT at in-focus.

Relation of our and previous OCT imaging theories
Zhou et al. has formulated the OCT signal as the collection of contributions from discrete scatterers (Section 7 of [12]), and gave a similar representation with our discrete equation of intensity-OCT [Eq. (6)]. Namely, they represented the OCT image and speckle as two independent additive components.
Our work added two points to their work. At first, we derived a continuous form of the intensity-OCT imaging equation [Eq. (10)]. This equation clarified that the meaningful OCT image is the convolution of the object [D(x)B 2 (x)] and a PSF of OCT intensity image P 2 a (x). Second, our work accounts for the phase of complex PSF. It has enabled the design of VDC (Section 4), interpretations of artifacts in Doppler OCT (Section 5.4), and understand the imaging property of computational CMD-OCT (Section 5.5).
Our OCT formulation does not represent the three-dimensional imaging property, similar to that discussed in Section 5.2 for VDC. It is because we have assumed that the lateral and axial imaging properties are independent of each other. Namely, low numerical aperture (NA) was tacitly assumed. Villiger et al. [24] and Ralston et al. [25] have formulated the three-dimensional imaging property of OCT, in which the lateral and axial imaging properties interact with each other. The extension of our OCT formulation to account for the axial imaging property and its interaction with the lateral imaging property might be a future work.
Another point to be discussed is the applicability of our formulation strategy to full-field (FF-) OCT. There are several works of formulation of FF-OCT imaging. For example, Tricoli et al. formulated imaging property of partially coherent FF-OCT [29], and Barolle et al. gave a theoretical model which recapitulates the manifestation of aberrations in FF-OCT [30]. Our formulation is potentially adoptable to FF-OCT. However, after Eq. (15), it is assumed that the illumination and collection pupil functions are identical. And this assumption is valid for point-scanning OCT but not true for FF-OCT with spatially incoherent light sources. We can make our theory applicable to FF-OCT by reformulating it from Eq. (10).

Conclusion
We presented a new mathematical model called DSM to represent a sample to be measured by OCT. Using this model, we reformulated lateral OCT imaging and also designed a new imaging method called VDC imaging.
Using the OCT reformulation, the meaningful signal and speckle were successfully formulated as two different mathematical entities. We also found that the meaningful signal and speckle originated from the auto-and cross-numerical interactions of the scatters' signal contribution, respectively.
VDC imaging is a combination of physical OCT imaging and subsequent complex signal processing, and we derived an imaging formula that expresses the physical and numerical process simultaneously. The imaging PSF that appeared in this formula suggested that VDC imaging yields an image similar to that yielded by DIC. We validated the VDC method by measuring biological samples.
Finally, the formulation strategy has the newly explained mathematical mechanism of the well-known but previously not well-explained artifacts of Doppler OCT, and the imaging property of CMD-OCT.
The presented formulation strategy may be applicable to other OCT-based imaging methods in the future. Although the present formulation focuses on the lateral imaging property, a similar approach based on DSM can be used to extend this formulation for axial and three-dimensional imaging properties.
OCT signal and its digitally laterally shifted complex conjugate as The first summation of the final expression (the third last line of the equation, blue) is for the case of i = j, i.e., for the AI group. The second summation (the last two lines of the equation, brown and red) is for the case of i ≠ j, that is, the CI group. Using the Dirac delta function, x p , and x q , the physical spatial coordinates in the sample, Eq. (28), are rewritten as The subtraction on the last line is because of the same reason as that in Eq. (7). By substituting the scatterer density maps [Eqs. (8) and (9)  Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.