Open Access
Add to BibSonomy
Abstract
Motion artifacts, from such sources as heartbeats, respiration, or peristalsis, often degrade microscopic images or videos of live subjects. We have developed a method using circular optical coherence tomography (OCT) scans to track the transverse and axial motion of biological samples at speeds ranging from several micrometers per second to several centimeters per second. We achieve fast and high-precision measurements of the magnitude and direction of the sample’s motion by adaptively controlling the circular scan pattern settings and applying interframe and intraframe analyses. These measurements are the basis of active motion compensation via feedback control for future in vivo microscopic and macroscopic imaging applications.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Motion artifacts, from such sources as heartbeats, respiration, or peristalsis, often degrade images or videos of live subjects. For example, during confocal or two-photon microscopy in neuroscience studies, sample motion artifacts hinder researchers in visualizing micrometer-level structures and recording long-period time-lapse data [1–3]. Moreover, since any biomedical imaging system has a confined field of view (FOV), severe sample motion can shift the target region out of the FOV during the image acquisition. Motion correction algorithms have been proposed to reduce and compensate the artifacts based on the acquired data [4–6]. However, these solutions may require multiple volumetric scans and the post-processing approach precludes real-time correction if the sample is moving out of the FOV during long-period acquisition. Specialized holders have been developed to limit the subject’s activity and keep the target within the FOV [7–9], but such artificial shackles can also change the behavior of living subjects [10,11].
One alternative is motion tracking, which can provide feedback to the imaging system’s control, adaptively following the subject within the FOV and decreasing or even eliminating motion artifacts. Several motion tracking algorithms [12–14] are in use for various applications, such as autonomous driving, but only a few [15–19] can be used at the microscopic level. An ideal motion tracking system for microscopy applications, without constraining the sample or the imaging system, should provide the high spatial and temporal resolution needed to achieve high measurement accuracy. Among many imaging modalities, optical coherence tomography (OCT), which is a non-invasive, label-free, real-time, three-dimensional (3D) imaging technique [20], is a promising candidate. Previous methods for tracking motion in OCT relied on speckle decorrelation analysis [21–25]. Liu et al. [23] explored the essential relationship between the displacement of two axial scans (A-scans) and their cross-correlation coefficient (XCC). Since the speckle size is on the same order of magnitude as the OCT resolution [21,26], speckle decorrelation analysis is more accurate than conventional frame-to-frame position tracking [12,13]. Speckle decorrelation analysis uses the principles of OCT signal formation and detection to connect two A-scans’ displacement and their XCC. We can derive the signal formation for each A-scan from the sample scattering distribution and the system’s point spread function [23]. Meanwhile, the XCC of two A-scans is calculated using their signal intensity profiles [26]. Therefore, the value of XCC can be used to estimate the position displacement between two A-scans.
Motion tracking based on speckle decorrelation analysis has been applied to track the nonconstant beam scanning speed of an OCT system [25]. It also has helped to correct motion artifacts caused by one-dimensional (1D) nonuniform probe scanning in hand-held OCT [22,23]. To track two-dimensional (2D) transverse motion, Liu et al. [24] scanned the sample circularly, capturing both the magnitude and the direction of the motion. However, their study reported the speed of the sample motion only in the range from 0.2 mm/s to 1.4 mm/s, a change of less than one order of magnitude. But actual sample motion often varies from the micro to the macro level, and tracking it requires a motion tracking method with a high dynamic range.
Here, we propose a method to track the transverse sample motion over an extended speed range, from several micrometers per second to several centimeters per second. We achieve this by adaptively controlling the circular scan pattern settings and applying interframe and intraframe analyses. In addition, we provide a detailed analysis of extracting the magnitude and direction of the sample motion. Further, we add axial motion tracking analysis to achieve 3D motion tracking. Several experiments validate our motion tracking method and demonstrate its ability to track programmed transverse motion patterns and the respiration of a mouse through the movements of the mouse’s skin.
2. Methods
When we repeatedly scan a light beam in a circle over a stationary sample (Fig. 1(A)), we will acquire the same image in each scan, and the distances between adjacent sampling points, or A-scans, on the circle are uniform. If the sample starts to move slowly in the transverse plane, or the XY plane (Fig. 1(B)), the low-speed movement will not significantly change the distances between the sampling points within the exposure time of a single A-scan. However, the transverse movement will accumulate over a full circular scan and cause the object to have a small lateral displacement, which will result in an offset between the second circular scan and the first one (Fig. 1(C)). If the sample moves fast, the relative speed between the sample motion and beam scanning speed will cause a large change, sufficient to change the spacing between A-scans within one circular scan (Fig. 1(D)). Hence, through precise measurements of the displacements between every circular scan or A-scan within each circle, we can accurately determine the speed and direction of the sample motion. Here, we propose two analysis models. One is called interframe analysis, and the other is called intraframe analysis. A detailed analysis of how to extract the magnitude and direction of sample motion from interframe and intraframe analyses is provided below.
Fig. 1. Circular scanning pattern over a sample. (A) Top view of the circular scan (red) on the sample (blue) when there is no sample motion. (B) Top view of the circular scan beam trace left on the sample when the sample is moving with a speed of ${v_m}$ at angle $\mathrm{\alpha }$. The green shape represents the position before the sample moved, and the blue represents the sample’s position after moving. (C) Beam traces of two successive circular scans when the sample moves along the positive $x$-axis. The red cycloid represents the first circular scan, and the blue cycloid represents the second. The blue and red beam spots in the black dashed rectangle are in the intersection between the two successive circular scans. (D) Representation of a single circular scan beam trace where ${v_m}$ has increased. The pink and green beam spots in the black dashed rectangles show the relative positions of adjacent A-scans in different regions along the trajectory.
Interframe analysis focuses on the intersection region between the adjacent circular scans and is aimed at tracking slow speed motion (micrometers per second to millimeters per second). Due to the slow motion, two successive circular scans will be offset, and intersect each other, where they scan over the same position. The black dashed rectangle in Fig. 1(C) outlines one of the intersection regions between adjacent circular scans. By measuring the intersection positions, we can derive the sample motion in the XY plane.
Intraframe analysis focuses on the displacement between adjacent A-scans within one circular scan, aiming at tracking fast motion (millimeters per second to centimeters per second). As the bottom black dashed rectangle in Fig. 1(D) shows, the displacement of two adjacent A-scans will become very close, because the beam scanning velocity and the sample motion have the same velocity direction (e.g., slow relative motion). Meanwhile, as the top black dashed rectangle in Fig. 1(D) shows, there will be little overlap between two adjacent A-scans when the beam scanning velocity and the sample motion have opposite velocity directions (e.g., fast relative motion). By analyzing the changes in the displacement of pairs of adjacent A-scans within one circular scan, we can learn how the sample moves in the XY plane.
To track axial motion, since each OCT A-scan contains depth information, we analyze the cross-correlation of images between successive circular scans to learn the displacement in the axial direction and further determine the sample’s axial motion.
With these analyses, we can track the speed and direction of the sample’s motion in 3D over a high dynamic range by repeatedly performing circular scans over the sample. Furthermore, because OCT has micron-scale resolution in both the axial and transverse dimensions and can acquire hundreds of thousands of A-scans per second, high spatial and temporal resolution measurements can be obtained.
2.1 Circular scan pattern settings
Three parameters control the circular scan pattern, the spacing between A-scans, and the OCT beam scanning speed: the radius of the circular scan (R, Fig. 1(A)), the number of sampling points (A-scans) acquired in one circular scan (NA), and the exposure time at each sampling point (${t_e}$). Thus, each circular scan will take a total time $T = {N_A}{t_e}$, which also determines the frame rate or temporal resolution of our motion tracking. The spacing between adjacent A-scans in the transverse plane is ${\delta _d} = \frac{{2\mathrm{\pi }R}}{{{N_A}}}$, and the beam scanning speed is ${v_s} = \frac{{2\mathrm{\pi }R}}{{{N_A}{t_e}}}$. Without loss of generality, we set our circular scan’s direction to be counterclockwise and its start phase to be zero.
2.2 Transverse motion tracking
2.2.1 Speckle decorrelation analysis
To calculate the displacement between adjacent A-scans from their intensity, we follow the speckle decorrelation analysis described in [23]. The XCC between two A-scans is calculated from their depth-dependent intensity profile, expressed by Eq. (1) [26], where N is the number of pixels per A-scan, $\zeta $ is the axial pixel index, ${I_i}(\zeta )$ and ${I_j}(\zeta )$ are the intensities of two different A-scans at axial pixel index $\zeta $, and $\bar{I}$ and ${\mathrm{\sigma }_I}$ are the mean and standard deviation of the A-scan intensity $I(\zeta )$.
The relationship between the square of the displacement (${d^2}$) between two A-scans and their XCC is expressed in Eq. (2) [23], where ${\mathrm{\omega }_0}$ is the OCT Gaussian beam waist as well as the transverse resolution of the OCT system.
Therefore, we can combine Eq. (1) and Eq. (2) to measure the ${d^2}$ between two A-scans by calculating their XCC from their intensity. Our motion tracking method will rely on this relationship.
2.2.2 Coordinate system setup for tracking transverse motion
We define the origin of the coordinate system to be the center of the circular scan. The x and y axes, shown in Fig. 1(A), define the transverse plane. As shown in Fig. 1(B), when the sample is moving at a speed ${v_m}$ in an angle $\mathrm{\alpha }$, a circular scan beam on the sample will trace a cycloid. Assuming the first A-scan started at $t = 0$, the coordinates $({x(t ),y(t )} )$ of the scanning beam at time t over the sample can be written as follows:
2.2.3 Interframe analysis
When a sample moves slowly, interframe analysis uses the intersection region between two successive circular scans to derive the magnitude and direction of the sample’s motion (Fig. 1(C)). For each A-scan on the first circular scan ${A_i}$ at ${t_i} = ({i - 1} ){t_e}$, we have the corresponding coordinates $({x({{t_i}} ),y({{t_i}} )} )$, and for each A-scan on the second circular scan ${A_j}$ at ${t_j} = ({j - 1} ){t_e} + T$, we have the corresponding coordinates $({x({{t_j}} ),y({{t_j}} )} )$. We use ${i^\ast }$ and ${j^\ast }$ to represent the corresponding indices of the intersecting A-scans between two circular scans. Hence, we have
Using $x({{t_{{i^\ast }}}} )\approx x({{t_{{j^\ast }}}} )$ from Eq. (5), we can derive ${v_{{m_x}}}$, which is the sample motion projected speed on the $x$-axis, where
Similarly, we can find ${v_{{m_y}}}$, which is the sample motion projected speed on the $y$-axis, from $y({{t_{{i^\ast }}}} )\approx y({{t_{{j^\ast }}}} )$, where
Thus, as $v_m^2 = v_{{m_x}}^2 + v_{{m_y}}^2$, we use ${{\textrm{Eq}.\,(7)} ^2} + {{\textrm{Eq}.\,(8 )}^2}$ to calculate the speed of sample’s motion:
Hence, the speed of the sample’s motion is written as:
For simplicity, we define the offset of the intersection indices as ${\mathrm{\varepsilon }^\ast } = {j^\ast } - {i^\ast }$. When ${N_A} \gg {\mathrm{\varepsilon }^\ast }$, we apply the first-order approximation and simplify Eq. (10):
Meanwhile, the angle can be calculated as $\tan (\mathrm{\alpha } )= \frac{{{v_{{m_y}}}}}{{{v_{{m_x}}}}} = \frac{{\textrm{Eq}.\; (8)}}{{\textrm{Eq}.\; (7)}}$, where
Therefore, the direction of the sample motion can be calculated from Eq. (12):
2.2.4 Intraframe analysis
When the sample speed increases, the distance between A-scans will no longer be uniform within one circular scan (Fig. 1(D)). Thus, to measure the sample’s motion, we switch to analyzing the spacing between adjacent A-scans within a single circular scan. We note the coordinates of $i$-th A-scan, $({x({{t_i}} ),y({{t_i}} )} )$, and the coordinates of $i + 1$-th A-scan, $({x({{t_{i + 1}}} ),y({{t_{i + 1}}} )} )$. The square of the distance ($d_{i,i + 1}^2$) between adjacent A-scans, ${A_i}$ and ${A_{i + 1}}$, on one circular scan is
From Eq. (3) and Eq. (4), we can further expand Eq. (14) to
Therefore, we can observe that $d_{i,i + 1}^2$ is a sine function of the A-scan index i, as
where $A = 4{v_m}{t_e}R\sin \left( {\frac{\mathrm{\pi }}{{{N_A}}}} \right)$, $C = 4{R^2}{\sin ^2}\left( {\frac{\pi }{{{N_A}}}} \right) + v_m^2t_e^2$, and $\mathrm{\psi } = \mathrm{\alpha } - \frac{{({2i - 1} )\mathrm{\pi }}}{{{N_A}}}$.Since ${N_A} \gg \mathrm{\pi }$, we have $\sin \left( {\frac{\mathrm{\pi }}{{{N_A}}}} \right) \approx \frac{\mathrm{\pi }}{{{N_A}}}$. We then use the amplitude of the sine function ($A$) to derive the speed of the sample’s motion,
Meanwhile, from the sine function, we can determine the index ${i^\ast }$ that reaches the minimum, where $\mathrm{\psi } = \frac{{3\mathrm{\pi }}}{2}.$ Thus, the direction of the sample’s motion is
2.3 Axial motion tracking
To track the axial motion, we analyze the cross-section image displacement between successive circular scans. Since the A-scan retrieves the depth information, if the sample motion has $z$-axis speed component, ${v_{{m_z}}}$, there will be a position offset $\Delta z$ from frame to frame.
Figure 2(A) is a flowchart of axial motion tracking with OCT. We first load two successive circular scans and designate the first one as our reference. On the second circular scan, we apply 2D XCC analysis [27] to calculate the similarity to the reference, searching from top to bottom. We record the displacement in the depth direction, $\mathrm{\Delta }z$, as where the maximum XCC position is. This displacement indicates that the sample moved by $\mathrm{\Delta }z$ in the axial direction after a circular scan cycle. Since each circular scan will take $T = {N_A}{t_e}$, the velocity of the sample’s axial motion can be expressed as
Fig. 2. Illustration of axial motion tracking. (A) Flowchart of axial motion tracking. (B) Sketch showing a cross-section view of a circular scan image, where the dashed line is the sample surface before movement and the solid line is sample surface after axial motion.
2.4 OCT system setup
A customized spectral domain OCT (SD-OCT) system, shown in Fig. 3, performs the circular scan for our motion tracking. We use an SLED (Exalos, EBD291023-02, ${\lambda _0} = 1300$ nm, $\mathrm{\Delta }\lambda = 175.6 \; \textrm{nm}$) as the broadband light source, which provides an axial resolution at about 5.3 $\mathrm{\mu}\textrm{m}$. We use a spectrometer (Cobra 1300, Wasatch Photonics) with a 2048-pixel InGaAs line-scan camera (Sensors Unlimited, GL2048) to measure the spectral interference pattern of the OCT signals. The maximum A-scan rate of the camera is 147 kHz, which provides a fastest A-scan exposure time of ${t_e} = 6.8\; \mathrm{\mu}\textrm{s}$. The slowest A-scan exposure is ${t_e} = 105\; \mathrm{\mu}\textrm{s}$. We use a 5X objective lens (Mitutoyo M Plan Apo NIR) to image the sample. The measured transverse resolution, which is the same as the beam waist ${\mathrm{\omega }_0}$, is 3.9 $\mathrm{\mu}\textrm{m}$.
The circular scan pattern projected over the sample is generated by applying two orthogonal sine wave voltages in the 2D galvanometer (GVS002, Thorlabs) on the sample arm, shown in Fig. 3. The voltage (${V_R}$) sent to the galvo system controls the radius of the circular scan, where $R = {k_V}{V_R}$. The coefficient ${k_V}$ represents the physical displacement of the probing beam at the focal plane corresponding to the galvo driving signal. It is calibrated by assigning different voltages to the 2D galvo system and measuring the corresponding physical scanning range on a US Air Force 1951 test target (DA052, Max Levy).
A two-axis precision motor stage (ALS130-100-NC-LTAS, Aerotech) and a Z-axis stepped motor stage (42BYG250Bk, Syntron) move the sample under the objective lens to control its speed and direction.
2.5 Sample preparations
To illustrate our motion tracking method, we used multilayer tapes and slices of defrosted chicken breast and onions as the imaging targets. Since mouse is a widely used experimental animal model in many applications, we also demonstrate our motion tracking by tracking a mouse’s breath under anesthesia by imaging the fluctuations of the mouse’s skin. A wild-type mouse was anesthetized using a mixture of Ketamine and Xylazine [28], had its hair removed exposing the skin, and was placed in a supine position over a heating pad to perform the experiments. Animal experiments were performed in accordance with protocols approved by the institutional animal care and use committee (IACUC) at Washington University in St. Louis.
3. Results
3.1 XCC sensitivity and noise floor
To verify how accurately and sensitively we can convert the XCC of A-scans to the square of their displacement, we first kept the sample stationary but varied ${N_A}$ from 100 to 8000 in one circular scan to change the density of the A-scans. From the geometry, we calculated the relationship between the square of the displacement of two adjacent A-scans in one circle ($d_{i,i + 1}^2$) and ${N_A}$, which is $d_{i,i + 1}^2 = \frac{{4{\mathrm{\pi }^2}{R^2}}}{{N_A^2}}$. Combining the result with Eq. (2), we get
We tested Eq. (20) with samples of tape, onion, and chicken breast. As seen in Fig. 4, for each circular scan, we calculated the average of the XCC for all pairs of adjacent A-scans. We also calculated the theoretical value of XCC based on Eq. (20) and plotted the result as the red line in Fig. 4. When ${N_A}$ is large enough, which implies ${d^2}$ between adjacent A-scans is small, the experimental XCC values of all three samples are close to the theoretical values. However, when ${N_A}$ is small, which implies ${d^2}$ between two A-scans is large, we observe mismatch. Due to the structure of the sample and the intrinsic noise of the OCT signal, an XCC noise floor (${\mathrm{\rho }_0}$) occurs where no matter how much we increase the displacement between two A-scans, and the value of XCC will never be lower than that floor. This noise floor affects the accuracy of the conversion from XCC to ${d^2}$. Hence, to avoid the noise floor effect, we must set circular scan pattern settings to give enough overlap between A-scans.
Fig. 4. The relationship between ${N_A}$ and XCC, showing the XCC noise floor for different samples.
3.2 Transverse motion tracking verification
3.2.1 Interframe analysis for slow motion
Using a slice of onion as a test example, we set the sample speed at 0.8 mm/s and the circular scan pattern settings as ${N_A} = 2000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = 0.9$ mm. We moved the sample along the positive $x$-axis to test the interframe analysis. Visualization 1 demonstrates the recording of several continuous circular scans for the interframe analysis. Figure 5(A) and (B) show two cross-section images of two adjacent circular scans. Highly similar structures are observed within the red and yellow rectangles. As indicated by the white arrows in Fig. 5(A), the structure in the red rectangle is shifted slightly to right in Fig. 5(B), while the structure in the yellow rectangle is shifted slightly to the left. These shifts indicate that the adjacent circular scans have scanned the same region over the sample, so the offset, or the intersection, has been clearly captured. To estimate the magnitude and direction of the sample motion from Eq. (11) and Eq. (13), we need to know the intersection indices for both the first circular scan and the second circular scan to derive the intersection index offset, ${\mathrm{\varepsilon }^\ast }$.
Fig. 5. Experimental results of interframe analysis in transverse motion tracking. (A) and (B) Cross-section images of the onion slice from the first and second circular scans. (C) Max XCC value for each $i$-th A-scan among all the $j$-th A-scans. (D) Detected intersection index offset ${\mathrm{\varepsilon }^\ast }$.
Based on Eq. (1), we computed the XCC between all pairs of A-scans in the first and the second circular scans. To find the intersection earlier described by Eq. (5), we recall that, according to Eq. (2), the higher the value of XCC, the closer the two A-scans will be. We first choose the $j$-th A-scan that most highly correlates $\textrm{with}$ each $i$-th A-scan, and we plot the results in Fig. 5(C), where the x-axis represents each $i$-th A-scan on the first circular scan, and the y-axis represents the maximum value of XCC among all the $j$-th A-scans. We applied XCC thresholding according to bimodal distribution to remove the background noise, where the A-scans have low correlation. Then, we can find the intersection index offset, ${\mathrm{\varepsilon }^\ast }$, plotted in Fig. 5(D), based on the corresponding ${i^\ast }$ and ${j^\ast }$.
We first validated Eq. (11) of the interframe analysis, which expresses the magnitude of slow sample motion. We fixed the sample speed at ${v_m} = 0.5$ mm/s, $R = 0.9$ mm, and ${t_e} = 20\,\mathrm{\mu}\textrm{s}$, but changed ${N_A}$ from 1000 to 8000 and plotted the relationship between $N_A^2$ and $|{{\mathrm{\varepsilon }^\mathrm{\ast }}} |$ in Fig. 6(A). We fixed ${N_A} = 6000$ and ${t_e} = 20\,\mathrm{\mu}\textrm{s}$, but changed R from 0.15 mm to 2.1 mm and plotted the relationship between $\frac{1}{R}$ and $|{{\mathrm{\varepsilon }^\mathrm{\ast }}} |$ in Fig. 6(B). Last, we fixed ${N_A} = 6000$ and $R = 0.9$ mm, but changed ${t_e}$ from 6.8 $\mathrm{\mu}\textrm{s}$ to 105 $\mathrm{\mu}\textrm{s}$, and plotted the relationship between ${t_e}$ and $|{{\mathrm{\varepsilon }^\mathrm{\ast }}} |$ in Fig. 6(C). To verify the use of different control variables, we fitted the experimental results and plotted them as red dashed line in each figure. The linearities in Figs. 6(A), 6(B), and 6(C) prove the correctness of Eq. (11).
Fig. 6. Verification of the interframe analysis. (A) Changing ${N_A}$. (B) Changing R. (C) Changing ${t_e}$ (D) Using different circular scan patterns to extract different detectable speed ranges of the sample motion. (E) Verification of detecting the direction of the sample motion.
Equation (11) also tells us that by adjusting different settings of the circular scan patterns, the interframe analysis can capture different ranges of sample motion. Therefore, we also tested three different circular scan pattern settings with three different speed ranges. We started with a circular scan pattern at ${N_A} = 6000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = 0.9$ mm. The blue dots with error bars in Fig. 6(D) show the experimental results, with sample motion speeds ranging from 10 µm/s to 0.9 mm/s. Our second circular scan pattern settings were ${N_A} = 4000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = 0.9$ mm. We varied the sample speed from 0.4 mm/s to 1 mm/s and plotted the experimental results as black dots with error bars in Fig. 6(D). The last circular scan pattern used settings of ${N_A} = 2000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = $ $0.9$ mm, tested over a sample speed range from 0.8 mm/s to 3 mm/s. The experimental results are plotted in green in Fig. 6(D). The red dashed fitted lines in Fig. 6(D) show the linearity between the sample speed and the absolute value of the intersection index offset, $|{{\mathrm{\varepsilon }^\ast }} |$, from Eq. (11). The three different circular scan pattern experiments show that by changing different settings of the circular scan patterns, with the interframe analysis, we can cover a detectable speed range from several micrometers per second to several millimeters per second. From Eq. (11) and our experimental results, we know that to detect slower speeds in this range we can increase ${N_A}$, decrease R, or increase ${t_e}$, while to detect faster speeds, we can decrease ${N_A}$, increase R, or decrease ${t_e}$.
To validate the direction detection estimated by the interframe analysis from Eq. (13), we fixed the sample motion speed at 0.5 mm/s and set our circular scan pattern at ${N_A} = 2000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = $ $1.5$ mm. We varied the direction of the sample motion over ${360^ \circ }$, and we recorded the corresponding ${i^\ast }$ and ${j^\ast }$ used in Eq. (13). We then calculated $\mathrm{\alpha }$ and plotted its changes corresponding to our experimental settings in Fig. 6(E). The red dashed line, $y = x$, shows the direction detection ability of the interframe analysis, with an average error of ${\sim} {1.3^ \circ }$.
3.2.2 Intraframe analysis for fast motion
Using a slice of defrosted chicken breast, we set the speed of sample motion at $7.5$ mm/s and set the circular scan pattern to ${N_A} = 10000$, ${t_e} = 60\; \mathrm{\mu}\textrm{s}$, and $R = 0.75$ mm. We moved the sample along the positive $x$-axis to demonstrate the intraframe analysis for fast motion. Figure 7(A) shows a single circular scan cross-section image of the chicken breast. The OCT image inside the red rectangle in Fig. 7(A), zoomed in view in Fig. 7(C), looks stretched because of the high concentration of A-scans within close positions, which are highly correlated. In contrast, the OCT image inside the green rectangle in Fig. 7(A), zoomed in view in Fig. 7(B), is well-delineated. The A-scans here are from distinct locations, which are lowly correlated. Based on Eq. (1), we calculated the XCC between all pairs of adjacent A-scans in one circular scan and plotted the results in Fig. 7(D). The maximum and minimum positions of the XCC value match our observation in Fig. 7(A). We further used Eq. (2) to convert the XCC into ${d^2}$ between the adjacent A-scans, as shown by blue dots in Fig. 7(E). We fitted the experimental results using a sine curve, as Eq. (16) describes and as shown by the red dashed line in Fig. 7(C). From the amplitude and the phase of the fitted sine curve, we then calculated the magnitude and direction of the sample motion.
Fig. 7. Experimental results of intraframe analysis. (A) Cross-section image of the chicken breast in one circular scan. (B) Zoomed in view of the green rectangle region where the speckles of the signal are well-delineated. (C) Zoomed in view of the red rectangle region where the speckles of the signal are stretched. (E) Curve of the value of XCC corresponds to (A). (E) Calculated ${d^2}$ value based on the value of XCC in (B). The scale bars in (A), (B) and (C) are 200 $\mathrm{\mu}\textrm{m}$.
To verify Eq. (17) of the intraframe analysis, where we derive the sample motion, we first set the sample speed $\textrm{to}\; {v_m} = 12$ mm/s and the circular scan pattern settings to $R = 1.2$ mm and ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, but changed ${N_A}$ from 6000 to 11000 and plotted the relationship between the amplitude of the sine curve of $d_{i,i + 1}^2$ (Eq. (16)) and $\frac{1}{{{N_A}}}$ in Fig. 8(A). We then set ${N_A} = 8000$ and ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, but changed R from 0.12 mm to 1.2 mm, and plotted the relationship between the amplitude of the sine curve and R in Fig. 8(B). Last, we set ${N_A} = 6000$ and $R = 0.9$ mm, but varied ${t_e}$ from 6.8 $\mathrm{\mu}\textrm{s}$ to 105 $\mathrm{\mu}\textrm{s}$, and plotted the relationship between ${t_e}$ and the amplitude of the sine curve in Fig. 8(C). In all cases, for the experimental results in Figs. 8(A), 8(B), and 8(C), we fitted them with the red dashed lines. The agreement of the linear relationship with the amplitude of the sine curve demonstrates the correctness of Eq. (17).
Fig. 8. Verification of the intraframe analysis. (A) Changing ${N_A}$. (B) Changing R. (C) Changing ${t_e}$. (D) Using different circular scan patterns to extract the different detectable speed ranges of the sample motion. (E) Verification of detecting the direction of the sample motion.
Similar to the interframe analysis, Eq. (17) in the intraframe analysis tells us that changing the circular scan patterns will lead to different ranges of detectable motions. As a test, we set three different circular scan patterns and varied the sample speed ranges. The first circular scan pattern settings were ${N_A} = 12000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = 1.2$ mm, and we varied the sample speed from 1 mm/s to 6 mm/s. The second pattern settings were ${N_A} = 6000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = 0.9$ mm and we varied the speed from 4 mm/s to 16 mm/s. The third pattern settings were ${N_A} = 5000$, ${t_e} = 20\; \mathrm{\mu}\textrm{s}$, and $R = $ $0.72$ mm, and we varied the speed from 14 mm/s to 22 mm/s. Figure 8(D) shows the amplitude of the sine curve (Eq. (16)) for each of these conditions in blue, black, and green dots with error bars. We also plotted the fitted lines in red dashed in Fig. 8(D) to show the linear correspondence between the sample speed and the amplitude of the sine curve from Eq. (17). The experimental results show that by changing different settings of the circular scan patterns, we can use intraframe analysis to cover the detectable speed range from several millimeters per second to several centimeters per second. From Eq. (17) and our results, we know that to detect the slower speeds in this range with intraframe analysis, we can decrease ${N_A}$, increase R, or $\textrm{increase}\; {t_e}$, while to detect faster speeds, we can increase ${N_A}$, decrease R, or decrease ${t_e}$.
Finally, to verify Eq. (18) for detecting the direction of the sample motion in the intraframe analysis, we fixed the speed of the sample motion at 7.5 mm/s and set our circular scan pattern at ${N_A} = 10000$, ${t_e} = 60\; \mathrm{\mu}\textrm{s}$, and $R = $ $0.75$ mm. We varied the direction of the sample motion over ${360^ \circ }$ and obtained the phase of the sine curve (Eq. (16)) from experimental data. We plotted the calculation results in Fig. 8(E) as blue dots with error bars, where $\textrm{the}\; x$-axis represents the change in the angle setting of the sample motion and the $y$-axis represents the angle change calculated from analysis. The red dashed line, a $y = x$, plotted in Fig. 8(E), shows the direction detection ability of the intraframe analysis, with an average error of ${\sim} {4.2^ \circ }$.
3.2.3 Transverse motion tracking extracts motion pattern
To show the capability of our method in the transverse plane, we designed a 2D transverse sample motion pattern, the joined letters “WU”, shown in Fig. 9(A) by the purple line and red arrows. This motion pattern demonstrates that our method can handle four basic motion conditions: 1) Both speed and direction change at the same time, 2) Speed changes, but direction does not, 3) Direction changes, but speed does not, and 4) Neither speed nor direction changes. We coded the motion pattern into the X-Y stage controller software to move a multilayer tape sample. Next, we started the circular scan acquisition and then started the stage movement. The circular pattern was set at ${N_A} = 4000$, $R = 0.9$ mm, and ${t_e} = 20\; \mathrm{\mu}\textrm{s}$. We calculated the magnitudes and directions of the sample motion, plotted in Fig. 9(B) and 9(C). We also plotted the designed values. The average accuracy of the tracked speed is 98.1%, where the error is ∼0.004 mm/s. The average accuracy of the tracked direction is ∼99.7%, where the error is ${\sim} {0.3^ \circ }$. With the calculated speed and angle, we drew the extracted sample motion pattern from our motion tracking in Fig. 9(A), with blue dots overlapping with the red designed pattern. Visualization 2 demonstrates this experiment in detail.
Fig. 9. Motion pattern extraction. (A) Designed and tracked motion pattern. (B) Designed and tracked magnitude of the sample motion. (C) Designed and tracked direction of the sample motion.
3.3 Axial motion tracking verification
To verify our axial motion tracking, we fixed the X-Y stage but set the Z-stage at different speeds to move the chicken breast sample up or down. To extract the $z$-axis speed component, ${v_{{m_z}}}$, we followed the flowchart shown in Fig. 2(A). Figure 10(A) and (B) are cross-section images of OCT circular scans before and after the Z-stage movement. We drew the red dashed line in Fig. 10(A) and 10(B) to represent the surface location before the axial motion. We also drew the yellow dashed line in Fig. 10(B) to represent the surface location after the axial motion. We can see an obvious axial shift, $\mathrm{\Delta }z$. Though we demonstrated the axial shift at the sample’s surface, we used the sample’s structural information near the surface for motion tracking. We changed the axial sample motion speed from $0.3$ mm/s to $2.7$ mm/s, then we plotted the axial pixel offset and the axial speed as blue dots with error bars in Fig. 10(C). We also plotted our fitted result as the red dashed line. The linearity again proves the correctness of Eq. (19), which shows the ability of our axial motion tracking.
Fig. 10. Verification of axial motion tracking. (A) and (B) Cross-section OCT circular scan images of the chicken breast before and after Z-stage moving. (C) Experimental results for different axial sample speeds.
3.4 Motion tracking of mouse breathing in vivo
Since a mouse usually has a very high breathing frequency, we needed a suitable frame rate for axial motion tracking. To track the mouse’s breath, we set the circular scan pattern at ${N_A} = 600$, $R = 0.15$ mm, and ${t_e} = 10\; \mathrm{\mu}\textrm{s}$. Thus, the frame rate is ∼166 Hz. We removed the mouse’s hair around its chest and performed circular scans over this region, as shown in Fig. 11(A). The axial motion tracking experimental results are shown in Fig. 11(B). The average breath rate calculated by the time differences from peak to peak is 195.22 breaths per minute, which matches a healthy mouse’s breath rate under anesthesia [28]. During breathing, the absolute axial motion speed of the chest reached ∼5.5 mm/s. Visualization 3 records the OCT circular scans during the mouse breathing. To test our motion tracking in 3D, we added a transverse motion to the mouse by moving the stage at 0.8 mm/s. Therefore, the XY plane motion was provided by the stage and the Z direction motion was provided by the mouse breathing. Considering the small size of the mouse’s chest, we applied interframe analysis and switched our circular scan pattern settings to ${N_A} = 1500$, $R = 0.3$ mm, and ${t_e} = 40\,\mathrm{\mu}\textrm{s}$. From Eq. (11), we then calculated the magnitude of the average transverse motion as 0.78 mm/s. These results demonstrated that our motion tracking could be applied to future in vivo imaging scenarios.
Fig. 11. Motion tracking of mouse skin. (A) Sketch of the mouse and where the circle scan was performed (red circle). The white position is where the hair on the chest has been removed. (B) Axial motion tracking results.
4. Discussion
Our results showed that by changing different settings of the circular scan pattern, both interframe and intraframe analysis could detect different speed ranges in transverse motion tracking. Furthermore, the detectable speed ranges overlapped between both analyses, which provided a bridge for switching different analysis approaches when the speed of the sample motion became faster or slower. However, our analyzes have minimum and maximum detectable speeds for any one circular scan pattern setting. When the magnitude of the sample motion exceeds the detectable speed range for that circular scan pattern, the scan pattern may be conveniently adjusted to accommodate the desired detection range.
The minimum detectable speed is limited by the interframe analysis, as determined by Eq. (11). Since the absolute value of the intersection index offset ${\mathrm{\varepsilon }^\ast }$ is an integer, its minimum value will be 1. In other words, if ${\mathrm{\varepsilon }^\ast }$ is lower than one, we cannot detect the motion between two adjacent circular scans, and the interframe analysis will assume that the sample does not move, i.e., two successive circular scans are almost identical. Therefore, we have the theoretical minimum detectable speed of interframe analysis, as expressed in Eq. (21):
If the speed that we want to detect is slower than what we calculated from Eq. (21), there are two options. The first is to analyze the third circular scan or even more subsequent circular scans and find the intersection region with the first circular scan. The other option is to consider changing the settings of the circular scan patterns. If we increase ${N_A}$ or ${t_e}$ to slow down the single circular scan period, the sample motion will traverse a noticeable distance between two adjacent circular scans. We can also decrease R or increase ${N_A}$ to get finer resolution for tracking the sample motion. However, since the distance of two adjacent A-scans on the circle is $\frac{{2\mathrm{\pi R}}}{{{N_A}}}$, we must be careful when changing the values of ${N_A}$ and R, where the pixel size must be greater than or equal to half of the beam waist ${\mathrm{\omega }_0}$ to match the sampling theorem and to avoid interference from the XCC noise floor. If we decrease the beam waist, in other words, increase the transverse resolution, we can obtain better resolution to detect lower speeds.
The maximum detectable speed is determined by the intraframe analysis and is limited by the XCC noise floor (${\mathrm{\rho }_0}$). When we convert the XCC from acquired data to the square of the distance (${d^2}$) between adjacent A-scans within one circular scan based on Eq. (2), once the relative speed between the sample motion and the beam scanning velocity is too large, the conversion will lose accuracy. To retrieve all the information of the fitted sine curve in Eq. (16) correctly, we need $\textrm{max}({d_{i,i + 1}^2} )\le d_0^2$, where ${d_0}$ is the maximum displacement that can be detected from the noise floor of the XCC, as calculated in Eq. (22):
Then, we have the maximum speed that can be detected derived in Eq. (23):
If the speed we want to detect is faster than Eq. (23), we also have two options for tracking the faster speed under two conditions. Intraframe analysis relies very heavily on a suitable relative speed between the beam scanning velocity and the sample motion. If the sample speed is slower than the scanning velocity, but it exceeds the value Eq. (23) defines, we must increase ${N_A}$ or decrease ${t_e}$ to reduce the beam scanning speed and get closer to the magnitude of the sample motion, or we can decrease R to increase the correlation of the two adjacent A-scans. In these ways, changing the circular scan pattern settings makes the sine curve of Eq. (16) better defined. If the sample speed is faster than the beam scanning speed, we can choose to focus on the minimum of the sine curve and its A-scan index, as Eq. (16) expressed. In this case, the maximum speed of detection can be extended to Eq. (24):
Two different analysis approaches and the adjustments of circular scan pattern settings are indispensable for our high dynamic range motion tracking. If we only use interframe analysis, when we decrease ${N_A}$ or increase R to detect high-speed motion, the displacement between A-scans will increase, and the intersection of two adjacent circles cannot provide enough overlap to calculate XCC accurately. Thus, the interframe analysis will lose its accuracy for high-speed motion. If we only use intraframe analysis, although the highest detection speed has been clarified in Eq. (23) and Eq. (24), when we decrease ${N_A}$ or increase R in order to test the low-speed motion, there will be insufficient overlap between adjacent A-scans for accurate XCC calculation. Both interframe and intraframe analysis are complementary in their speed detection range.
In the actual application of our 3D motion tracking, circular scanning will be continuously performed over the sample. We may not need to always start the analysis from the place where the start phase of scanning is equal to 0, and we also may not need to always scan counterclockwise. Introducing a new start phase variable or clockwise scanning direction will not change the essence of our analysis models. Prior knowledge or estimation of the measured speed range can be helpful in setting the initial parameters and choosing whether to use interframe or intraframe analysis. In practice, for the unknown speed range, scan parameter settings may be adaptively changed to achieve optimal performance. Scan parameters may be fine-tuned to clearly distinguish the index offset in interframe analysis or the amplitude of the sine curve in intraframe analysis. Meanwhile, the choice of interframe or intraframe analysis also depends on which method can more clearly extract the magnitude and direction of the sample motion.
The result of our motion tracking is the average value within the data acquisition time, which is also determined by our frame rate or temporal resolution. It equals the reciprocal of the multiplication of the adjustable parameters ${N_A}$ and ${t_e}$. When integrating our method with the feedback loop control in the future, to reduce the latency, we can consider integrating the analyses on the graphics processing unit (GPU) or field programmable gate array (FPGA).
5. Conclusions
This work demonstrates a 3D motion tracking method based on circular scans with OCT. We developed two analysis models, intraframe and interframe analysis, to cover a broad range of detectable speeds. Our experimental results demonstrated the accuracy of the motion tracking method and showed its large-scale detectable speed range, extending from several micrometers per second to several centimeters per second in transverse and axial directions. Further, we discussed the strategies for changing the circular scan parameters under two different analysis models to cover different speed ranges in tracking sample motion. The instructions for use will be very helpful when combining this auxiliary motion tracking system with the primary imaging system. In the future, automatic motion tracking will be implemented through feedback control, and both transverse tracking and axial tracking analysis can be carried out independently and compensated separately.
One limitation of this sample motion tracking method is that, within the circular scan region, we assume that the sample is rigid and has only three rectilinear velocity vectors. However, the sample may rotate, deform, or change its topology, so that not all A-scans will have the same axial shift. To overcome this limitation, parallel OCT imaging [29,30] may be performed on different regions of the sample to obtain a motion vector field in the future. Also, it is possible to combine the motion tracking methods with concentric-circle scanning OCT [31,32]. By changing the scan radius over time, after reconstructing the A-scans’ positions with motion tracking compensation, motion-artifact-free 3D OCT images can be obtained.
Funding
National Institutes of Health (R01EB025209, R01HL156265, R21EB03268401A1); Washington University in St. Louis (Start-up fund).
Acknowledgments
The authors wish to express their gratitude to Fei Wang for her valuable advice on video generation, Chao Ren for his critical suggestions regarding the validation experiments, and Mr. James Ballard for his constructive feedback on the manuscript.
Disclosures
SH and CZ are inventors of a patent application submitted by Washington University in St. Louis regarding methods described in this work.
Data availability
Data presented in this paper are available upon reasonable request to the authors.
References
1. D. Soulet, A. Paré, J. Coste, and S. Lacroix, “Automated filtering of intrinsic movement artifacts during two-photon intravital microscopy,” PLoS One 8(1), e53942 (2013). [CrossRef]
2. R. Weigert, N. Porat-Shliom, and P. Amornphimoltham, “Imaging cell biology in live animals: ready for prime time,” J. Cell Biol. 201(7), 969–979 (2013). [CrossRef]
3. M. S. Creamer, K. S. Chen, A. M. Leifer, and J. W. Pillow, “Correcting motion induced fluorescence artifacts in two-channel neural imaging,” PLoS Comput. Biol. 18(9), e1010421 (2022). [CrossRef]
4. M. F. Kraus, B. Potsaid, M. A. Mayer, R. Bock, B. Baumann, J. J. Liu, J. Hornegger, and J. G. Fujimoto, “Motion correction in optical coherence tomography volumes on a per A-scan basis using orthogonal scan patterns,” Biomed. Opt. Express 3(6), 1182–1199 (2012). [CrossRef]
5. M. F. Kraus, J. J. Liu, J. Schottenhamml, C.-L. Chen, A. Budai, L. Branchini, T. Ko, H. Ishikawa, G. Wollstein, J. Schuman, J. S. Duker, J. G. Fujimoto, and J. Hornegger, “Quantitative 3D-OCT motion correction with tilt and illumination correction, robust similarity measure and regularization,” Biomed. Opt. Express 5(8), 2591–2613 (2014). [CrossRef]
6. Y. Chen, Y.-J. Hong, S. Makita, and Y. Yasuno, “Three-dimensional eye motion correction by Lissajous scan optical coherence tomography,” Biomed. Opt. Express 8(3), 1783–1802 (2017). [CrossRef]
7. R. G. Presson, M. B. Brown, A. J. Fisher, R. M. Sandoval, K. W. Dunn, K. S. Lorenz, E. J. Delp, P. Salama, B. A. Molitoris, and I. Petrache, “Two-photon imaging within the murine thorax without respiratory and cardiac motion artifact,” Am. J. Pathol. 179(1), 75–82 (2011). [CrossRef]
8. W. Li, R. G. Nava, A. C. Bribriesco, B. H. Zinselmeyer, J. H. Spahn, A. E. Gelman, A. S. Krupnick, M. J. Miller, and D. Kreisel, “Intravital 2-photon imaging of leukocyte trafficking in beating heart,” J. Clin. Invest. 122(7), 2499–2508 (2012). [CrossRef]
9. W. Zong, R. Wu, S. Chen, J. Wu, H. Wang, Z. Zhao, G. Chen, R. Tu, D. Wu, Y. Hu, Y. Xu, Y. Wang, Z. Duan, H. Wu, Y. Zhang, J. Zhang, A. Wang, L. Chen, and H. Cheng, “Miniature two-photon microscopy for enlarged field-of-view, multi-plane and long-term brain imaging,” Nat. Methods 18(1), 46–49 (2021). [CrossRef]
10. Z. J. Hall and V. Tropepe, “Movement maintains forebrain neurogenesis via peripheral neural feedback in larval zebrafish,” eLife 7, e31045 (2018). [CrossRef]
11. M. Loos, B. Koopmans, E. Aarts, G. Maroteaux, S. van der Sluis, N.-B. M. P. Consortium, M. Verhage, and A. B. Smit, “Sheltering behavior and locomotor activity in 11 genetically diverse common inbred mouse strains using home-cage monitoring,” PLoS One 9(9), e108563 (2014). [CrossRef]
12. A. Cavallaro and E. Maggio, Video Tracking: Theory and Practice (John Wiley & Sons, 2011).
13. E. Trucco and K. Plakas, “Video tracking: a concise survey,” IEEE J. Oceanic Eng. 31(2), 520–529 (2006). [CrossRef]
14. A. P. Shukla and M. Saini, “Moving object tracking of vehicle detection: a concise review,” Int. J. Signal Process. Image Process. Pattern Recognit. 8, 169–176 (2015). [CrossRef]
15. T. T. Wu and J. Y. Qu, “Optical imaging for medical diagnosis based on active stereo vision and motion tracking,” Opt. Express 15(16), 10421–10426 (2007). [CrossRef]
16. K. V. Vienola, B. Braaf, C. K. Sheehy, Q. Yang, P. Tiruveedhula, D. W. Arathorn, J. F. de Boer, and A. Roorda, “Real-time eye motion compensation for OCT imaging with tracking SLO,” Biomed. Opt. Express 3(11), 2950–2963 (2012). [CrossRef]
17. R. D. Ferguson, D. X. Hammer, L. A. Paunescu, S. Beaton, and J. S. Schuman, “Tracking optical coherence tomography,” Opt. Lett. 29(18), 2139–2141 (2004). [CrossRef]
18. P. Mecê, J. Scholler, K. Groux, and C. Boccara, “High-resolution in-vivo human retinal imaging using full-field OCT with optical stabilization of axial motion,” Biomed. Opt. Express 11(1), 492–504 (2020). [CrossRef]
19. Y. Cai, K. Grieve, and P. Mecê, “Characterization and analysis of retinal axial motion at high spatiotemporal resolution and its implication for real-time correction in human retinal imaging,” Front. Med. 9, 1 (2022). [CrossRef]
20. 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(5035), 1178–1181 (1991). [CrossRef]
21. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]
22. A. Ahmad, S. G. Adie, E. J. Chaney, U. Sharma, and S. A. Boppart, “Cross-correlation-based image acquisition technique for manually-scanned optical coherence tomography,” Opt. Express 17(10), 8125–8136 (2009). [CrossRef]
23. X. Liu, Y. Huang, and J. U. Kang, “Distortion-free freehand-scanning OCT implemented with real-time scanning speed variance correction,” Opt. Express 20(15), 16567–16583 (2012). [CrossRef]
24. X. Liu and J. U. Kang, “Tracking both magnitude and direction of 2D transverse motion with optical coherence tomography,” in Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XVIII (SPIE, 2014), Vol. 8934, pp. 30–36.
25. Y. Wang, Y. Wang, A. Akansu, K. D. Belfield, B. Hubbi, and X. Liu, “Robust motion tracking based on adaptive speckle decorrelation analysis of OCT signal,” Biomed. Opt. Express 6(11), 4302–4316 (2015). [CrossRef]
26. J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199 (1998). [CrossRef]
27. D. Freedman, R. Pisani, and R. Purves, Statistics: Fourth International Student Edition (W.W. Norton & Company, 2007).
28. D. Wellington, I. Mikaelian, and L. Singer, “Comparison of ketamine–xylazine and ketamine–dexmedetomidine anesthesia and intraperitoneal tolerance in rats,” J. Am. Assoc. Lab. Anim. Sci. 52, 1 (2013).
29. C. Zhou, A. Alex, J. Rasakanthan, and Y. Ma, “Space-division multiplexing optical coherence tomography,” Opt. Express 21(16), 19219–19227 (2013). [CrossRef]
30. Y. Huang, M. Badar, A. Nitkowski, A. Weinroth, N. Tansu, and C. Zhou, “Wide-field high-speed space-division multiplexing optical coherence tomography using an integrated photonic device,” Biomed. Opt. Express 8(8), 3856 (2017). [CrossRef]
31. M. Niederleithner, M. Salas, R. A. Leitgeb, W. Drexler, and T. Schmoll, “Spiral scanning OCT angiography,” Invest. Ophthalmol. Vis. Sci. 60, PB070 (2019).
32. B. Tan, Z. Hosseinaee, and K. Bizheva, “Dense concentric circle scanning protocol for measuring pulsatile retinal blood flow in rats with Doppler optical coherence tomography,” J. Biomed. Opt. 22(11), 110501 (2017). [CrossRef]
