# Data Stream Stabilization for Optical Coherence Tomography Volumetric
Scanning

Guiqiu Liao ID , Oscar Caravaca-Mora ID , Benoit Rosa, Philippe Zanne,
Alexandre Asch, Diego Dall’Alba, Paolo Fiorini, Michel de Mathelin, Florent
Nageotte and Michalina J. Gora Manuscript received January 31st, 2021; revised
May 12th, 2021, accepted June 15th, 2021. This work was supported in part by
the European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie grant agreement No 813782 - project ATLAS. Guiqiu Liao
is with ICube, UMR 7357 CNRS-University of Strabourg, Strasbourg, France, and
also with Department of Computer Science, University of Verona, Verona, Italy
(e-mail: liao.guiqiu@etu.unistra.fr). Oscar Caravaca-Mora, Benoit Rosa,
Philippe Zanne, Alexandre Asch, Michel de Mathelin, Florent Nageotte and
Michalina J. Gora are with ICube, UMR 7357 CNRS-University of Strabourg,
Strasbourg, France (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>gora@unistra.fr).Diego Dall’Alba and Paolo Fiorini are
with Department of Computer Science, University of Verona, Verona, Italy
(e-mail<EMAIL_ADDRESS>paolo.fiorini@univr.it).

###### Abstract

Optical Coherence Tomography (OCT) is an emerging medical imaging modality for
luminal organ diagnosis. The non-constant rotation speed of optical components
in the OCT catheter tip causes rotational distortion in OCT volumetric
scanning. By improving the scanning process, this instability can be partially
reduced. To further correct the rotational distortion in the OCT image, a
volumetric data stabilization algorithm is proposed. The algorithm first
estimates the Non-Uniform Rotational Distortion (NURD) for each B-scan by
using a Convolutional Neural Network (CNN). A correlation map between two
successive B-scans is computed and provided as input to the CNN. To solve the
problem of accumulative error in iterative frame stream processing, we deploy
an overall rotation estimation between reference orientation and actual OCT
image orientation. We train the network with synthetic OCT videos by
intentionally adding rotational distortion into real OCT images. As part of
this article we discuss the proposed method in two different scanning modes:
the first is a conventional pullback mode where the optical components move
along the protection sheath, and the second is a self-designed scanning mode
where the catheter is globally translated by using an external actuator. The
efficiency and robustness of the proposed method are evaluated with synthetic
scans as well as real scans under two scanning modes.

###### Index Terms:

Optical Coherence Tomography, Video Stabilization, Non-uniform Rotational
Distortion, Convolutional Neural Network.

##  I Introduction

Optical Coherence Tomography (OCT) is an emerging medical imaging modality for
diagnosing diseases in cardiovascular  system, airways and gastrointestinal
tract [1]. The main  challenges related to these  applications are delivery
and collection of light from the volumetric scanning of the area of interest.
Fourier domain OCT in a single measurement provides a one-dimensional in-depth
information and in order to collect three-dimensional information the optical
beam  has to be scanned over the tissue. Various designs of specialized micro-
optical catheters have been developed to achieve this goal [2].

A circumferential two-dimensional scan can be performed by rotation of an
optical beam reflected to the side of the probe using a micro-motor on the
distal tip, or by a proximal rotational actuation, which is remotely connected
to the distal optical components with a torque coil (see Fig. 1). Volumetric
scanning, in both cases, is then typically effectuated by pulling back the
rotating optical core to create a helical scan. This approach for volumetric
scanning was originally developed for cardiovascular applications where an OCT
catheter is inserted using a guidewire into the location of interest and a
pullback needs to be synchronized with a blood occlusion [3]. It was then
adopted in gastrointestinal imaging both in low-profile and balloon catheters,
which are inserted in the digestive system using  a working channel of an
endoscope [4]. Most recently, an internal pullback scanning system was also
developed for a  tethered capsule device [5]. A high precision short segment
pullback enabled high quality en-face imaging that could not be achieved with
a standard tethered capsule devices typically pulled back manually.

Figure 1: Schematics of (A) proximal-scanning with a fiber-optic rotary joint;
(B) distal-scanning endoscope with a micro-motor.

Even though motorized rotational and pullback actuation provides good control
of the OCT beam position, scanning distortions are still very common [4, 5].
Non-Uniform Rotational Distortion (NURD) is crucial for catheter-based imaging
systems such as endoscopic OCT but also intravascular ultrasound [6, 7].

In proximal rotation devices NURD can be observed due to mechanical friction
between rotating optical components and the protecting sheath, mostly caused
by bending of the catheter inserted in the tortuous path of the digestive,
cardiovascular or respiratory system. Miniaturization of motors enabled distal
scanning and more compact [1] and more easily miniaturized designs[11], in
which instabilities are only related to the quality of the motor itself.
However, these designs require electrical connection of the distally located
motor, which adds additional safety risks to the device and  causes partial
shadowing of the tissue by the wires. Non-linear transmission of displacement,
and imperfect synchronization between data acquisition and scanning speed also
causes distortions of the longitudinal scanning. This is typically observed as
drift in the 3D reconstruction.

We have previously investigated the NURD performance of an OCT catheter that
was developed by our team for real-time assistance during minimally invasive
robotized treatment of colorectal cancer [47]. This robotized OCT catheter can
be inserted in one of the two instrument channels of the robotized
interventional endoscope and teleoperated using a master controller. In this
article, we focus on studying the volumetric imaging instability of proximal
scanning endoscopic OCT catheters, and on stabilizing the OCT image by
proposing both a new robotic 3D scanning and a software approach. A new
scanning method was motivated by redundancy of two degrees of freedom (DoF) in
the OCT enhanced robotized interventional endoscope. The OCT catheter once
inserted in the robot has two DoF for volumetric imaging (rotation and
pullback of the inner optical core, provided by the OCT scanning system) and 3
DoF for maneuverability (rotation, translation and bending, provided by the
robotics system [48]). As can be seen both rotational and translation are
redundant, but OCT rotational scanning requires at least 3000 rpm, which
cannot be achieved by a motorized instrument driver. The required translation
speed of 2.5 mm/s, on the other hand, can be achieved by both inner core
pullback or outer tool pullback. The outer pullback can reduce the instability
of proximal scanning with a moving torque coil, but not perfectly, because the
residual friction with the protecting sheath cannot be completely eliminated,
and thus we propose a software correction to further decreased distortions.

When using the information within OCT images, two challenges need to be
addressed for the software  correction of the NURD. The first challenge is to
compensate the NURD on-line, which means that the algorithm should only take
historical OCT data to stabilize the scanning video. This can allow the
catheter to correctly access the cross-sectional structure of tissues in real-
time, which can dramatically reduce the time it takes for medical
doctors/physicians to make a diagnostic decision. Unlike the off-line
stabilization that uses the entire data after a scan to optimize the warping
parameters, on-line NURD correction typically follows an iterative processing
pipeline where the estimation error accumulates. In this case, common
iterative algorithms fail when a certain frame has a large stabilization error
or when the processing time is long. Using optical markers is an alternative
choice to stabilize rotational distortion, but it will block the OCT laser
beam and degrade the image quality.  Thus the second challenge  is that the
algorithm should not use any type of additional objects to track the
orientation shifts.

We formulate the problem as finding a warping vector for each frame, to shift
each individual scanning A-line to its correct position. This can be done
based on computing the cross-correlation or Euler distance between A-lines in
the new uncorrected frame and in previous stabilized frames [8, 9, 10, 11].
However, in a real application it is hard to ensure that every A-line of OCT
images can provide visible features for correlation because of the limited
Field of View (FoV) of OCT catheters (part of the OCT image is filled with
background noise when the tubular lumen is larger than the catheter FoV).
Given the prior knowledge [8, 9] that the warping path should be continuous
through the whole OCT frame stream, the problem of searching the path in a
computed correlation matrix is highly similar to a continuous boundary line
prediction problem. Motivated by this, we seek to leverage the state of art
machine learning techniques for the boundary line prediction method [22, 15]
as a solution to optimize the warping vector estimation.

Although the warping vector estimation accuracy can be improved, the
accumulative drift still exists in long data stream processing. We address
this issue by estimating an overall rotation value for each frame, and
compensating the drift error of the NURD warping vector using this overall
rotation value. To do this, the algorithm exploits the consistency of the OCT
protecting sheath as a cue for overall rotation estimation, and this overall
rotation is a direct and robust estimation (without time-integration process),
although it only provides a low accuracy because it relies on limited
information. We deploy a fusion algorithm [29] to integrate the overall
orientation and element-wise NURD estimation, which is inspired the fusion of
gyroscope and accelerometry for attitude angle estimation [30].

We demonstrate our method in two different 3D scanning modes: internal
pullback scanning and outer tool scanning (see Fig. 2 (A) and Fig. 2 (B)). For
the internal pullback scanning, we deploy an additional volumetric sheath
registration. We train the CNN based NURD estimation algorithm with semi-
synthetic data by intentionally distorting OCT images. We evaluate the
accuracy of the NURD estimation with semi-synthetic data, and test its
robustness on real data. An ablation study is demonstrated by disabling parts
of the full algorithm.

##  II Related Work

In this section we provide an overview of previous research on NURD correction
for catheter-based imaging systems, followed by an introduction to the state-
of-the-art video stabilization research and Convolutional Neural Network (CNN)
research for white light camera, which inspired the proposed method for
endoscopic OCT stabilization.

###  II-A NURD Correction

Several NURD correction methods have been proposed for catheter-based imaging
systems [6, 8, 9, 10, 11, 12, 14]. A frequency analysis of the texture of the
Intra-Vascular Ultrasound (IVUS) image was used in [6] to estimate the
rotational velocity. Matching based algorithms by minimizing the Euler
distance of warped neighboring frames [8, 9], or by maximizing the cross-
correlation [10, 11, 14] have been proposed to track the shift error of
A-lines. Specifically, in intravascular OCT imaging, intra-vascular stents may
serve as landmarks that help to capture rotational distortion in pullback
volumetric scanning [13]. However, matching based methods using cross-
correlation or Euler distance typically require highly correlated images. Even
with tractable features between successive frame pairs, it is difficult to
guarantee that the estimation of NURD is correct, which leads to iterative
drift in long data streams processing. Adding structural markers on the OCT
sheath provides an extrinsic reference for orientation change[9], but the
markers block the OCT light and thus remove information about the tissue.

###  II-B Data Fusion for Rotation Estimation

Similar to iterative NURD estimation, the estimation of attitude angle with
integral gyroscope data also has the problem of drift [27, 28]. The integral
drift of a gyroscope is usually compensated by another robust angle estimation
from sensors such as accelerometers and magnetometers [30, 31, 32, 33]. While
a gyroscope provides excellent information about rapid orientation changes, it
only provides relative orientation changes that gradually drift with their
lifetime and temperature. An accelerometer or magnetometer, on the other hand,
has a direct measurement of orientation but with lower accuracy. Various
classes of filters were demonstrated to fuse accurate rapid relative
(indirect) measurements with less accurate direct measurements such as with
the Extended Kalman Filter (EKF) [33], the complementary filter[29, 31, 32]
and a gradient based filter [30]. In a similar way to the role of an
accelerometer or a magnetometer, another additional overall rotation can be
estimated and fused with the NURD estimation, which compensates the
accumulative error.

###  II-C Video Stabilization

In the research field of image processing, CNN and deep learning related
approaches stand as state-of-the-art. For OCT image processing, CNN has been
applied to tissue layer segmentation [17, 18, 19], classification [20] and
cancer tissue identification [21], but there is no CNN application for OCT
stabilization. On the other hand, the literature in video stabilization is
richer for white light cameras than for medical imaging systems (including the
OCT).  We seek to fill in the gap between the common computer vision research
field and that of the OCT imaging, by relying on the CNN to enhance the
efficiency of OCT frame stream stabilization.

For white light camera video stabilization, there are two types of approaches
to model the problem. One seeks to directly estimate the camera path
(position), and the video stabilization can be considered as a camera path
smoothing problem [36]. This formulation aims to stabilize homographic
distortion cause by camera shaking, and recently a deep learning based method
has been developed to learn from data registered by a mechanical stabilizer
[37], which shows greater efficiency than traditional algorithms. The other
type of approach models the instability of the video (or frame stream) as an
appearance change [38]. This modeling methodology can be adapted to different
imaging systems beyond the white light camera. To formulate the appearance
change, features matching algorithms or optical flow [39, 40] can be used. A
recent study uses deep learning techniques to estimate an optical flow field
representing shift map of pixels in the video frames, and then applies another
CNN regression module to estimate a pixel-wise warping field from the optical
flow field [35]. Inspired by this, we deploy a CNN to estimate the element-
wise NURD warping vector from a correlation map which roughly interprets the
NURD distortion of a single frame.

###  II-D CNN for path searching

An essential step of the NURD estimation is to search a continuous optimal
path with large correlation value from the correlation map. Solutions applied
in previous OCT stabilization studies [10, 11, 14] are mainly based on graph
searching  (GS), and rely on local features of gradients, maxima, textures and
other prior information. This type of technique is also a traditional way of
contour tracing [23]. Recently deep learning based contour tracing techniques
[24, 25, 26, 22] has been demonstrated to be faster and more robust than
traditional methods.

The state-of-the-art deep learning models for pixel-wise segmentation are
based on adaptations of convolutional networks that had originally been
designed for image classification. To solve dense prediction problems such as
semantic segmentation , which are structurally different from image
classification, striding and dilated CNN [41] are proposed to systematically
aggregate multiscale contextual information without losing resolution. Path
detection can be achieved with a pixel-wise segmentation architecture (_e.g._
predict a binary map where the path position and background pixels have
different values). This is a high-cost approach , which usually adopts a
U-shape CNN [25, 26, 22] using up-convolution layers[42, 43]. We deploy a CNN
to predict a single vector representing path coordinates from the correlation
map. Compared with pixel-wise prediction, the proposed architecture can be
more efficient and no post processing is required. Moreover, it is easy to
integrate the continuity prior knowledge into the network training, by an
additional term in the loss function. Since we consider the situation when
part of the image has no feature, a simple cascade layer architecture [45] is
not sufficient for spacial exploitation. Therefore we use a parallel scheme
[44] which controls striding of convolution to extract hierarchical abstract
information from the correlation matrix.  Moreover, the proposed network has
shallower layers than VGG [45] or Resnet [46] and thus requires less inference
time. Our tests and results indicate that this choice is sufficient for the
NURD stabilization task.

##  III Methods

###  III-A OCT volumetric imaging

To compare the conventional pullback and robotic 3D scanning performance we
use an OCT instrument with an outer diameter of 3.5 mm [47]. The instrument is
terminated at the distal end with a transparent sheath on the tip, which
allows three-dimensional OCT imaging using an internal rotating side-focusing
optical probe with two proximal external scanning actuators. The instrument is
connected to an OCT imaging system built around the OCT Axsun engine, with a
1310 nm center wavelength swept source laser and 100 kHz A-line rate. The OCT
catheter is compatible with an instrument channel of a robotized flexible
interventional endoscope [48]. The distal end of the robotized endoscope can
be bent in 2 orthogonal planes, translated and rotated. The OCT catheter can
be translated, rotated and bent in one plane. In addition, the rotation and
translation of the inner OCT probe can also be controlled by a servo system.

We test two different types of volumetric scans with this system. One is an
internal pullback scanning where the torque coil of the catheter moves through
the protecting sheath, and another one is a robotic outer pullback of the
whole tool where the torque coil delivers only rotational scanning. Fig. 2
shows a schematic of the internal pullback scan and the outer tool scan.

Figure 2: (A) Schematic of internal pull back for proximal-scanning OCT; (B)
Schematic of robotic endoscopic scanning; (C) The rectangular test object; (D)
An example of OCT B-scan, the scale bar corresponds to 1 mm.

###  III-B Data Stream Stabilization

A  typical OCT scanning system makes the assumption that the focusing lens at
the distal tip rotates with a constant speed, so the OCT data acquisition
system arranges A-lines of each frame with equal spacing to cover a 360
degrees region in the polar domain. However, for various reasons such as
friction and non-constant motor speed the angular speed of the distal optical
components may not be uniform over time, resulting in a distortion of the
reconstructed B-scan in polar domain. The OCT rotational distortion correction
algorithm estimates a warping vector representing the angular position error
of each independent A-line. This vector is used to place A-lines at new
positions in order to limit the NURD artifact.

Figure 3: Algorithm pipeline of the proposed OCT stabilization method. The
original new frame is first pre-shifted based on the previous estimated
overall rotation.  Then the pre-shifted OCT image follows two paths in
parallel to estimate two levels of artifacts. In the first path (with blue
color) a correlation map between the pre-shifted frame and last stabilized
frame is computed, and a CNN estimates an A-line level NURD vector. In the
second path (with red color), the sheath image of the pre-shifted frame is
used to match the reference sheath stack to estimate an overall rotation
value. The final warping vector is a fusion of the A-line level NURD vector
(from the first path) and the single rotation value (from the second path),
which is used to shift each individual A-line of the new frame to obtain a new
corrected frame.

The architecture of the proposed stabilization algorithm is shown in Fig. 3,
and is composed  of two parts. The first part iteratively estimates the NURD
vector ${{\bar{P}}_{k}}$ with the last stabilized frame
${{\hat{\boldsymbol{F}}}_{k-1}}$ and the new frame
${{\tilde{\boldsymbol{F}}}_{k}}$, which is pre-shifted with a previous
estimated overall rotation ${{\bar{r}}_{k-1}}$. The second part estimates the
new overall rotation ${{\bar{r}}_{k}}$ with the pre-shifted
${{\tilde{\boldsymbol{F}}}_{k}}$ and a reference frame stack. The NURD warping
vector ${{\bar{P}}_{k}}$ is fused with ${{\bar{r}}_{k}}$, to compensate the
iterative drift of NURD estimation. By doing so, a final warping vector
${{\hat{P}}_{k}}$ is estimated. ${{\hat{P}}_{k}}$ is equal in length to the
polar image height $H$, which is the number of A-lines in a single OCT frame.
Details on these two parts and the fusion algorithm are presented in the
following two subsections.

####  III-B1 NURD estimation

The first part of the stabilization algorithm follows two steps to estimate
the NURD vector between the previous stabilized
${{\hat{\boldsymbol{F}}}_{k-1}}$ and pre-shifted frame
${{\tilde{\boldsymbol{F}}}_{k}}$: First a correlation map
${{\boldsymbol{M}_{k}}}$ is calculated using cross-correlation, and then a CNN
estimates the NURD vector taking the ${{\boldsymbol{M}_{k}}}$ as input.

Given two consecutive OCT frames with height $H$ (with $H$ A-lines in one
frame) and width $W$ in polar coordinates, we used the Pearson correlation
coefficient to reflect the similarity between two local image patches from
these two frames:

${{s}_{i,j}}=\frac{\sum\nolimits_{l=1}^{n}{{{\boldsymbol{f}}_{i,l}}{{\boldsymbol{f}}_{j,l}^{\prime}}-n{{{\bar{\boldsymbol{f}}}}_{i}}{\bar{{{\boldsymbol{f}}}_{j}^{\prime}}}}}{\sqrt{\sum\nolimits_{l=1}^{n}{\boldsymbol{f}_{i,l}^{2}-n\bar{\boldsymbol{f}}_{i}^{2}}}\sqrt{\sum\nolimits_{l=1}^{n}{{\boldsymbol{f}}_{j,l}^{\prime
2}-n{\bar{\boldsymbol{f}_{j}^{\prime}}}^{2}}}}$ (1)

where $\boldsymbol{f}_{i}\in{{\mathbb{R}}^{h\times W}}$ ($h\ll H$) is a local
image patch centered at index position $i,i\in[0,H)$ of the newest frame.
Correspondingly, ${{{\boldsymbol{f}}}_{j}^{\prime}}\in{{\mathbb{R}}^{h\times
W}}$ is a local image patch centered at index position $j$ of
${{\hat{\boldsymbol{F}}}_{k-1}}$ ( $j\in[i-w/2,i+w/2)$, $w$ is is the shifting
window width). ${s}_{i,j}$ is one element of the correlation map
${{\boldsymbol{M}_{k}}}\in{{\mathbb{R}}^{H\times w}}$. The pixel index $l$
will operate through the patch size $n=h\times W$. $\bar{\boldsymbol{f}_{i}}$
and $\bar{\boldsymbol{f}_{j}^{\prime}}$ are the mean values of patches
${\boldsymbol{f}_{i}}$ and ${\boldsymbol{f}_{j}^{\prime}}$ respectively.

Figure 4:  Schematics of the CNN for NURD estimation. This network has two
sub-branches that have different convolution kernels and striding.

The correlation map ${{\boldsymbol{M}_{k}}}$ contains information on the NURD,
and we design a CNN to further extract the NURD vector. The detailed diagram
of the NURD estimating nets is shown in fig. 4, where each block gives the
convolutional kernel size, the stride step length and the output feature map
depth of each convolution layer. Two convolution sub-branches are operated in
parallel to capture features and produce hierarchically coarse-to-fine
responses. The main difference between these two sub-branches is the strides
control. Both the left sub-branch and the right sub-branch of the NURD
estimation nets have 6 convolutional layers, and a LeakyReLU activation [49]
is used after each convolution layer.

The left sub-branch always keeps the vertical stride as 1. The objective of
this design is to emphasize the vertical spacial correspondence. By doing so,
the front 5 feature extraction layers can gradually reduce the feature map
width from $w$ to 1, while maintaining the feature map height $H$ as the
input’s height. The depth of each convolution operation’s output is twice as
deep as its input (here we set the output depth of the first layer as 8). The
$5th$ feature map $\boldsymbol{A}_{F}^{5}\in{{\mathbb{R}}^{H\times 1\times
128}}$ extracts 128 local features, which could include the minimal value
position, edge, and boundary position. A final layer with kernel size $1\times
1$ and channel depth 128, reorganizes the $5th$ feature map and decrease
channels to a sub-branch output of $\bar{P}$ with size $H\times 1$.

Considering that in certain situations ${{\boldsymbol{M}_{k}}}$ will miss
valid information for some row $m_{i}$ when there is no feature in a patch
(window) $\boldsymbol{f}_{i}$ of ${{\tilde{\boldsymbol{F}}}_{k}}$, the right
sub-branch is added. This sub-branch loosens the horizontal stride to 2 to
involve more spacial information. This is a more common architecture for many
deep learning applications [45]. Compared to the left sub-branch, this sub-
branch can extract higher abstract features which are less sensitive to local
missing information in the correlation map.  We put larger kernels in the
front layers to involve more common spatial information for neighboring
positions of convolution, which prevents the under-fitting of feature
extraction [16]. The kernel size is reduced in deeper layers to ensure a good
trade-off between kernel size and computational burden. After the six
convolutions, this subbranch outputs a matrix
$\bar{\boldsymbol{P}}_{m}^{\prime\prime}$ of size $(H/2^{6})\times w$. To get
an output of warping vector size $H\times 1$, a reshape operation is applied
to the output of this sub-branch (_i.e._ by sequentially connecting each
$1\times w$ row). A final convolution layer with kernel size $3\times 1$,
input depth 2 and output depth 1, combines the prediction of two sub-branches
to get an estimation of $\bar{P}$. The loss function for training the NURD
estimating nets uses the conventional $L_{2}$ loss and an additional
continuity loss. A standard $L_{2}$ loss is described by Eq. (2):

${{L}_{2}}=\frac{1}{{{n}_{p}}}{{\sum\nolimits_{i=1}^{{{n}_{p}}}{({{P}_{i}}-{{{\bar{P}}}_{i}})}}^{2}}$
(2)

where $P_{i}$ is the element of $P$, the true NURD vector (ground truth),
$n_{p}=H$ is the vector length. The $L_{2}$ loss function is commonly used for
value estimation, while for this estimation task, to take into account the
prior knowledge on the continuity of the distortion vector [8, 9, 10], we add
the continuity loss as follows:

$L_{c}=\frac{1}{{{n}_{p}}-1}{{\sum\nolimits_{i=1}^{{{n}_{p}-1}}{({{\bar{P}}_{k,i}}-{{{\bar{P}}}_{k,i+1}})}}^{2}}$
(3)

By calculating $L_{c}$, and combining it with $L_{2}$ in the network training,
the attraction towards local minima with discontinuous NURD estimation will be
suppressed. The final loss for branch (A) is:

$L_{A}=\alpha L_{c}+(1-\alpha)L_{2}$ (4)

####  III-B2 Drift compensation

It is possible to correct a short frame stack with just the NURD estimation
part presented in the previous subsection. To compensate the accumulative
error through iterations (especially for a long data stream), we estimate an
overall rotation value ${{\bar{r}}_{k}}$.

In conventional pullback scans only the sheath data is consistent no matter
what the environment outside the sheath is. According to this, the overall
rotation can be observed by matching real-time sheath images with pre-recorded
reference sheath images. However, when using the unstabilized pullback
scanning to record the reference images, it still suffers from rotational
distortion. To ensure  that the orientations of reference frames are correct,
we follow a calibration procedure. The setup of the sheath registration is
shown in Fig. 5 (A), which relies on a external calibration object (a straight
and flat ruler). As shown in Fig. 5 (B), the raw reference data still has the
rotational distortion, which leads to both the ruler and sheath images located
at wrong direction. We extract the contour surface of the ruler, and align the
raw reference frame stack by minimizing the surface distance of all frames. By
doing so, the rational error of the raw reference volumetric data is reduced
from 59.4∘ to 2.79∘ (see Fig. 5 (C)). This calibrated reference data composes
one of the inputs of the overall rotation estimation. Another input is
composed by real-time B-scans, where image outside the sheath is masked out
and only the sheath part is used.

With a recorded sheath image stack $\boldsymbol{S}\in{{\mathbb{R}}^{H\times
W\times N}}$ (N is the number of frames in the entire reference stack), we
compute Euler distances to estimate the rotation deviation $\Delta\bar{r}_{k}$
of pre-shifted image ${{\tilde{\boldsymbol{F}}}_{k}}$ (pre-shifted by
$\bar{r}_{k-1}$):

Figure 5: (A) Setup for sheath registration; (B) Original registration data in
3D view; (C) Calibrated registration data in 3D view.

${{d}_{r}}=\mathcal{L}({{\tilde{S}}_{k}},S_{k})$ (5)

where ${{\tilde{S}}_{k}}\in{{\mathbb{R}}^{H\times W\times l}}$ is buffered
pre-shifted original images, $S_{k}\in{{\mathbb{R}}^{H\times W\times l}}$ is
the corresponding selected buffer indexed at $k$ of $\boldsymbol{S}$ ($l\ll
N$, and assuming that in every pullback scanning the torque coil moves with
equal interval). The function $\mathcal{L}$ calculates Euler distance of two
input buffers. For the pullback scan, $S_{k}$ is the corresponding reference
image in the data stack, while to correct an outer robotic tool pullback 3D
scan the algorithm can be simplified and just  use the first frame
${{\boldsymbol{F}}}_{0}$ as the reference, because in this scanning condition
the relative longitude position between sheath and rotating ball-lens is
constant. For both ${{\tilde{S}}_{k}}$ and $S_{k}$ only the sheath part is
used for distance calculation. The index $r\in(-w,w)$ shifts within a given
window $w$ to calculate the distance $d_{r}$, and the position with the lowest
distance value is considered as the matched angular position of
$\Delta\bar{r}_{k}$. Finally, the overall rotation used to compensate the NURD
drift is computed by $\bar{r}_{k}=\bar{r}_{k-1}+\Delta\bar{r}_{k}$.

We use the concept of a _PI Complementary Filter_ [31] to fuse the
$\bar{P}_{k}$ vector with the $\bar{r}_{k}$ value, which is shown to be
efficient and able to solve the warping vector integral drift problem. A
discrete form of PI complementary filter is:

${{\hat{P}}_{k}}={{k}_{p}}{\bar{P}_{k}}+(1-{{k}_{p}}){{\bar{r}}_{k}}\boldsymbol{1}+{{k}_{i}}{{I}_{k}}$
(6) ${{I}_{k}}={{I}_{k-1}}+({{\bar{r}}_{k}}\boldsymbol{1}-{{\hat{P}}_{k}})$
(7)

where $k_{p}$ and $k_{i}$ are PI compensating gains. ${{I}_{k}}$ is the
integral component vector. $\boldsymbol{1}$ is a vector of ones. Each element
${{\hat{P}}_{k,i}}$ of $\hat{P_{k}}$ represents the angular distance between
the position of the $i$th A-line of ${{\tilde{\boldsymbol{F}}}_{k}}$ and its
correct position in polar domain. Applying $\hat{P_{k}}$ to this raw frame
${{\tilde{\boldsymbol{F}}}_{k}}$, a stabilized frame
${{\hat{\boldsymbol{F}}}_{k}}$ is obtained.

###  III-C Experiment

####  III-C1 Data

We test the proposed stabilization algorithm for both the internal
conventional pullback and  the outer scanning using a custom made OCT system
(see description in section III.A). We collect videos with the testing object
(fig. 2 (C)) to evaluate the performance for both types of scanning mode,
since it has a special geometry. 6 scans were collected for each scanning
mode, and each data stream comprises 500 frames. For the internal pullback,
one reference scan is collected following the procedure described in section
III.B.(2), and this reference scan is used in the correction of all internal
pullback scans. To perform quantitative analysis, we also use synthetic videos
generated from a variety of OCT images [50, 51, 52] to test the proposed
algorithm, since it is almost impossible to manually annotate the NURD in the
OCT data. We used the published videos from [50, 51, 52], and extracted all
the OCT images out and transformed them into the same resolution in the polar
domain. Approximately 5000 images which cover the catheter based OCT imaging
in intravascular, digestive tract and lung respiration airway are obtained
from these public videos.

####  III-C2 NURD estimation networks training

The network is only trained with synthetic videos by intentionally adding NURD
artifacts to the Optical Coherence Tomography (OCT) images, and tested on both
true OCT scans and synthetic scanning videos.  To ensure that the artificially
added distortion covers the distribution of artifacts present in real videos,
we first compute the correlation map $\boldsymbol{M}_{k}$ of every neighboring
raw image pairs in real videos, and then measure the maximum NURD error with a
graph path searching method [11]. The range of the GS-observed NURD was then
expanded by 1/3 on each side. For each synthetic image, the synthetic NURD
vector was then randomly generated within this expanded range. In addition, we
implemented data augmentation that included geometric transformations, noise
addition, brightness/contrast modification and OCT speckle/shadow simulation.
By doing so the domain of artifacts of synthetic videos should cover the
distribution of real artifacts. We randomly split all the collected images by
1:1:1 into train, validation and test data. For both the CNN training and
algorithm deployment, we use the same PC with an Nvidia Qt1000 graphic card
and Intel i5-9400H CPU. The code is implemented using the Pytorch
framework[53]. For training, Batch Normalization (BN) is used right after
convolution and before activation [54], dropout is not used [58] and weight
initialization is performed following the method described in [57]. The CNN is
trained with the Stochastic Gradient Descent (SGD) weights optimization method
(we used a weight decay of 0.0001 and a momentum of 0.9).

####  III-C3 Metrics

For synthetic videos we calculate the NURD estimation error to evaluate the
accuracy. For true videos, no guaranteed ground truth is available. We
calculate the normalized Standard Deviation (STD) $\sigma$ for a certain time
window of videos to access the stabilization performance [8]. The definition
of STD is:

$\sigma=\frac{1}{N_{sig}}\sum\nolimits_{i=1,j=1}^{N_{sig}}{\bar{\sigma}(f_{i,j})}$
(8)

where ${\bar{\sigma}(f_{i,j})}$ is the STD calculated with pixel $f_{i,j}$ in
one stacked frame stream, $i$ and $j$ are selected pixel indices on horizontal
and vertical axis respectively. $N_{sig}$ is the number of pixels used to
calculate $\bar{\sigma}$. Based on the STD calculation, we give a threshold to
differentiate unstable pixels from stable pixels, and count the number of
unstable pixels for each frame.

Figure 6: Result presented in RGB encoded image, 3D view and en-face
projection. (A) shows results of an original scan of rectangular phantom and
(B) shows the corresponding stabilized data. (C) shows results of an original
scan of layered colon tissue phantom and (D) shows the corresponding
stabilized data.  The yellow arrows point out rotation artifacts represented
by colorful pixels in B-scans, the yellow arrowheads point out colorful pixels
caused by the features changes of the object itself, and the red asterisks
mark out aliasing artifacts in en-face projections.

To evaluate the rotation error through a complete scan, we compute en-face
projections [10, 11] of the OCT videos where each A-line is accumulated to one
single value, so that the volumetric data in the polar domain are projected
into 2D images. In this case the vertical Y axis corresponds to a
circumferential scanning (B-Scan) and the horizontal X axis to a longitudinal
volumetric scanning (3D Scan). By counting the shift of max intensity pixels
in a whole stack, we can evaluate the precession angle. We also measure the
local fluctuation on the en-face image , which represents the local shaking
between neighboring frames.

##  IV Results

###  IV-A Pullback actuation comparison

To compare  distortions we effectuated two  proposed types of scans within a
21 mm long rectangular tube target 3D printed in polyjet VeroWhite material
(Stratasys, Rehovot, Israel), as shown in Fig. 2 (C). Table I summarizes
statistical analysis of distortion of these 2 types of scans. We have also
compared them to  a stationary scanning without moving the optical components
on the catheter tip, which is used as a baseline for hardware instability. It
can be seen that the stationary scan data has lower precession, local rotation
angle and STD (we also separate the sheath image STD from the full image STD)
compared to any of the volumetric scanning.

TABLE I:  Statistics of different scanning methods Method |  Preces- sion(∘)$\downarrow$ |  local (∘)$\downarrow$ |  STD$\downarrow$ |  Sheath STD$\downarrow$
---|---|---|---|---
Stationary |  11.6 |  12.9$\pm$2.5 |  16.6$\pm$2.2 |  13.3$\pm$2.5
Conventional |  79.6 |  13.6$\pm$3.1 |  18.6$\pm$2.4 |  19.8$\pm$2.6
Robotic 3D |  24.7 |  12.3$\pm$1.9 |  17.2$\pm$2.9 |  16.9$\pm$2.6

These two volumetric scanning methods introduce instability into the data, and
in particular the conventional internal pullback scanning presents larger
image rotation between successive 2D radial frames, with a maximum precession
angle of 79.6°. When performing volumetric robotic scanning of the OCT tool
the image rotation (precession) is reduced to within the range of 24.6°. This
constitutes an improvement and reduces the difficulty of software
stabilization (see section III.B). Other metrics such as the local rotation
and STD reflect the short term instability. Although the 3 scanning methods
share similar short term instability, the stationary scan is sightly better
than the internal pullback scan and robotic endoscopic scan.

###  IV-B Software Improvement

####  IV-B1 Quantitative results

The proposed stabilization algorithm is applied to both scanning modes. A
pixel-wise stack smoothing is deployed to be the baseline method (average
filtering of each pixel using 3 consecutive frames).

Table II shows the precession of the whole data stream, local fluctuation
angle and mean pixel counts (MPC) of unstable pixels. For both scanning modes,
the precession angle is reduced by 75%-80%, and for the robotic scan the
precession is stabilised to only 6∘. The baseline method can reduce the
intensity derivation by blurring the 3D data, but does not improve the
rotations and angular fluctuations. The proposed method not only reduces the
overall rotations, but also reduces the pixel-wise instability without
blurring the image or changing other information within every OCT frame (see
the images in the subsection for qualitative results).

TABLE II:  Comparison with baseline method  Method | Prec(∘) $\downarrow$ | local (∘)$\downarrow$ | MPC($\times 10^{3}$)$\downarrow$
---|---|---|---
Original (conventional) | 79.6 | 13.6$\pm$3.1 | 8.72$\pm$4.46
Baseline (conventional) | 79.1 | 13.1$\pm$3.2 | 2.94$\pm$2.31
Proposed (conventional) | 16.0 | 2.61$\pm$0.4 | 1.61$\pm$0.72
Original (robotic) | 24.7 | 12.3$\pm$1.9 | 10.9$\pm$5.90
Baseline (robotic) | 24.6 | 12.1$\pm$2.3 | 3.79$\pm$2.83
Proposed (robotic) | 6.05 | 3.02$\pm$0.8 | 2.76$\pm$1.03

####  IV-B2 Ablation Study

TABLE III:  ablation study for two scanning modes and synthetic scans Method | Prec(∘) $\downarrow$ | MPC($\times 10^{3}$) $\downarrow$
---|---|---
Proposed full (conventional) | 16.0 | 1.61$\pm$0.72
w/o overall (conventional) | 21.6 | 2.28$\pm$0.95
w/o NURD (conventional) | 23.4 | 6.86$\pm$4.32
Proposed full (robotic) | 6.05 | 2.67$\pm$1.03
w/o overall (robotic) | 61.4 | 6.02$\pm$3.22
w/o NURD (robotic) | 13.4 | 8.91$\pm$4.83
original (synthetic) | 132 | 19.8$\pm$5.30
Proposed full (synthetic) | 2.61 | 2.60$\pm$1.43
w/o overall (synthetic) | 136 | 3.01$\pm$1.72
w/o NURD (synthetic) | 12.9 | 5.45$\pm$2.95

Figure 7: Ablation comparison results for scans of rectangular phantom. RGB
encoded images and en-face projections of (A) original data, (B) stabilized
data without (w/o) overall rotation estimation, (C) stabilized data with only
overall rotation estimation and (D) stabilized data of the full proposed
algorithm.  The yellow arrows point out rotation artifacts represented by
colorful pixels in B-scans, the yellow arrowheads point out colorful pixels
caused by the feature changing of the object itself, and the red asterisks
mark out aliasing artifacts in en-face projections.

Figure 8: Ablation comparison results for scans of sponge surface. RGB encoded
images and en-face projections of (A) original data, (B) stabilized data
without overall rotation estimation, (C) stabilized data with only overall
rotation estimation and (D) stabilized data of the full proposed algorithm.
The yellow arrows point out rotation artifacts represented by colorful pixels
in B-scans, the yellow arrowheads point out colorful pixels caused by the
feature changing of the object itself, and the red asterisks mark out aliasing
artifacts in en-face projections.

Figure 9: Comparison with GS method. The left figure shows the max precession
angles of difference data output, the right figure shows the in-stable pixel
MPC comparison.

We use both the conventional internal pullback and the robotic scan and the
synthetic video to study the functionality of the two parts of the algorithm
(the element-wise NURD estimation part, and the overall rotation estimation
part). We did a thorough ablation comparison by disabling any of those two
parts, and the results are shown in Table III. The most detrimental
modification is to disable the overall estimation part. In this case the
proposed algorithm can no longer provide a valid angular value to compensate
the drift error, and the precession error is no longer corrected. Sometimes it
will even introduce a precession larger than the original precession, for
example in the case of  the real robotic scanning and synthetic videos.
Disabling the element-wise NURD estimation also affects the reduction of the
precession angle, but not as much as disabling the overall estimation.
Compared to our full stabilization algorithm, disabling the NURD estimation
obviously increases the unstable pixel count, because the overall rotation
estimation can not align individual A-lines. This is worse for real scans
where the overall rotation may have a larger estimation error due to the
sparse features in sheath matching.

####  IV-B3 Qualitative results

Fig. 6 shows qualitative results of pullback scans. To illustrate the NURD
instability, we encode each 3 consecutive frames in independent channels of
RGB color images[8]. In this RGB encoding, a colorful part of the image
indicates the NURD artifact, and a well stabilized frame sequence should have
minimal colorful pixels. We transform the data stream into the en-face
projection (see description in section IV.A). All these formats of results are
presented in Fig. 6. It can be seen that the “colorful” pixels of the original
sequence (RGB encoded) is reduced by the proposed algorithm, and without
changing the image quality. The 3D results are consistent with the en-face
projections, where the angular rotation is reflected by the shift of high
intensity pixels of en-face images. It can be seen that the stabilization
algorithm straightens the high intensity lines of the original en-face
projection, which maintains the intensity distribution of the whole data
stream.  More qualitative visual comparisons can be seen in Supplementary
Movies.

Qualitative results on the ablation study are presented in fig. 7 and 8.
Keeping only the NURD estimation will introduce iterative distortion to the
image, and also limit the precession reduction (fig. 7 (B) and fig. 8 (B)).
Fig. 7 (C) and fig. 8 (C) show the results with only the overall rotation
estimation from the reference data stream. In this case the shaking of the
frame sequence can be enlarged due to a large angular estimation error. The
fusion stage from the full algorithm can effectively limit these phenomena.

TABLE IV: NURD estimation evaluation  Method |  MSE(Deg2)$\downarrow$ |  Prec. Reduction (%) $\uparrow$ |  MPC Reduction (%) $\uparrow$
---|---|---|---
GS [10] |  6.47 $\pm$ 1.62 |  12.85 |  70.26 $\pm$ 14.3
Proposed |  0.38 $\pm$ 0.09 |  78.86 |  77.78 $\pm$ 8.90

####  IV-B4 Comparison with state of the art

In addition, we compare the NURD estimation part of the proposed CNN based
method and a  state-of-the-art approach based on the graph path searching (GS)
[11],  which is not a learning based method. Since this online correction
method has no drift compensation, we disable the overall rotation estimation
to compare the NURD estimation performance, using the indirect metrics mean
pixel account (MPC) for in-stable pixels of real scans. We also compare them
for the synthetic scans, using mean square error (MSE) of NURD estimation as a
direct metric. Results from original scans and different algorithms are shown
in fig. 9. Table IV shows statistics of these comparisons, the NURD estimation
error of the proposed method is significantly lower than graph searching based
method and can better reduce the number of unstable pixels (MPC).

##  V Conclusion

We have developed a new solution to tackle the rotational distortion problem
using deep CNN, which can be generalized for scanning situations with
different targets and catheters. We proposed a new warping vector estimation
net to estimate NURD between adjacent frames, which has a higher accuracy and
robustness compared with the conventional approach in situations where the
images have few features. Moreover, we solved the problem of drift error
accumulation in iterative video process, with a group rotation estimation net.
We were able to apply the CNN based algorithm trained on synthetic data to
real videos acquired in various scanning conditions. Compared to the A-line
level NURD vector estimated with tissue information, the orientation
estimation using the sheath is less accurate because the feature information
provided by the sheath is even more sparse. Nevertheless, this overall
rotation is robust and will not be affected by integral drift. The overall
rotation is a redundant complementary estimation besides the NURD estimation.
Fusing the NURD warping vector with the overall rotation value can form a
robust A-line level warping vector which suppresses the integral drift.

An ablation study has been performed on different scanning data by calculating
a quantitative stabilization metric and qualitative analysis, which shows that
both the NURD estimation and overall rotation are essential for the
stabilization of the OCT volumetric scanning. The proposed algorithm
outperforms other state-of-the-art methods[11] that iteratively estimate the
warping path using a graphic searching approach without a drift compensation
module.

In the literature of OCT stabilization, very little research study can be
found, and existing hardware solutions are mainly based on optical marker for
calibration [10]. On the other hand, software methods are for stationary
scanning without  longitudinal translation[8] or did not consider the
precession [11]. These image based methods will certainly be affected by
accumulative error in the volumetric scanning. The proposed stabilization
algorithm is designed for both conventional pullback scan and free 3D scan.
For the conventional pullback scan we design a sheath registration and
calibration procedure to record reference data. The new proposed 3D scan with
robotic tool pullback in a endoscope channel provides two benefits for OCT
volumetric imaging. The first is that it reduces the unstable friction and
reduces the precession from hardware level. The second benefit is that the
sheath in the first frame can be reference for overall rotation estimation,
which reduces the procedures of the further stabilization algorithm.

The proposed method assumes that the appearance change between consecutive
images is more affected by the rotational artifacts than by real changes of
tissue. The performance of the algorithm can be impaired in some cases where
the appearance of tissue in the field of view changes rapidly, for example
when the probe is brought to the tissue very fast or in the presence of
peristalsis. However, in the majority of clinical situations, to allow
thorough inspection of the tissue and to avoid its damaging, the probe is
maneuvered slowly, and thus the algorithm will achieve good stabilization.

## Acknowledgment

This work was supported by the ATIP-Avenir grant, the ARC Foundation for
Cancer research, the University of Strasbourg IdEx,  Plan Investissement
d’Avenir and by the ANR (ANR-10-IAHU-02 and ANR-11-LABX-0004-01) and funded by
ATLAS project from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 813782.

## References

  * [1] M. Gora, M. Suter, G. Tearney and X. Li, “Endoscopic optical coherence tomography: technologies and clinical applications [Invited]”, _Biomedical Optics Express_ , vol. 8, no. 5, pp. 2405-2444, 2017.
  * [2] M. Atif, H. Ullah, M. Y. Hamza, and M. Ikram, “Catheters for optical coherence tomography,” Laser Physics Letters, 2011.
  * [3] T. Okamura et al., “In vivo evaluation of stent strut distribution patterns in the bioabsorbable everolimus-eluting device: an OCT ad hoc analysis of the revision 1.0 and revision 1.1 stent design in the ABSORB clinical trial”, _EuroIntervention_ , vol. 5, no. 8, pp. 932-938, 2010.
  * [4] H.-C. Lee, O. O. Ahsen, K. Liang, Z. Wang, C. Cleveland, L. Booth, B. Potsaid, V. Jayaraman, A. E. Cable, H. Mashimo, R. Langer, G. Traverso, and J. G. Fujimoto, “Circumferential optical coherence tomography angiography imaging of the swine esophagus using a micromotor balloon catheter,” _Biomed. Opt. Express_ , 7(8), pp. 2927–2942, 2016.
  * [5] K. Liang, et. al, “Ultrahigh speed en face OCT capsule for endoscopic imaging”, _Biomedical optics express_ , 6(4), pp.1146-1163, 2015.
  * [6] Y. Kawase, Y. Suzuki, F. Ikeno, R. Yoneyama, K. Hoshino, H. Q. Ly, G. T. Lau, M. Hayase, A. C. Yeung, R. J. Hajjar, and I.-K. Jang, “Comparison of nonuniform rotational distortion between mechanical IVUS and OCT using a phantom model,” _Ultrasound in Medicine & Biology_, vol. 33, no. 1, pp. 67–73, 2007.
  * [7] W. Kang, H. Wang, Z. Wang, M. W. Jenkins, G. A. Isenberg, A. Chak, and A. M. Rollins, “Motion artifacts associated with in vivo endoscopic OCT images of the esophagus,” _Optics Express_ , vol. 19, no. 21, p. 20722, 2011.
  * [8] G. V. Soest, J. Bosch, and A. V. D. Steen, “Azimuthal Registration of Image Sequences Affected by Nonuniform Rotation Distortion,” _IEEE Transactions on Information Technology in Biomedicine_ , vol. 12, no. 3, pp. 348–355, 2008.
  * [9] O. O. Ahsen, H.-C. Lee, M. G. Giacomelli, Z. Wang, K. Liang, T.-H. Tsai, B. Potsaid, H. Mashimo, and J. G. Fujimoto, “Correction of rotational distortion for catheter-based en face OCT and OCT angiography,” _Optics Letters_ , vol. 39, no. 20, p. 5973-5976, 2014.
  * [10] N. Uribe-Patarroyo and B. E. Bouma, “Rotational distortion correction in endoscopic optical coherence tomography based on speckle de-correlation,” _Optics Letters_ , vol. 40, no. 23, p. 5518-5521, 2015.
  * [11] E. Abouei, A. M. D. Lee, H. Pahlevaninezhad, G. Hohert, M. Cua, P. Lane, S. Lam, and C. Macaulay, “Correction of motion artifacts in endoscopic optical coherence tomography and auto fluorescence images based on azimuthal en face image registration,” _Journal of Biomedical Optics_ , vol. 23, no. 01, p. 1-14, 2018.
  * [12] S. Sathyanarayana, “Nonuniform rotational distortion (NURD) reduction,” U.S. Patent 7 024 025 B2, April 4, 2006.
  * [13] G. J. Ughi, “Automatic three-dimensional registration of intravascular optical coherence tomography images,” _Journal of Biomedical Optics_ , vol. 17, no. 2, p. 026005, 2012.
  * [14] C. Gatta, O. Pujol, O. Leor, J. Ferre, and P. Radeva, “Fast Rigid Registration of Vascular Structures in IVUS Sequences,” _IEEE Transactions on Information Technology in Biomedicine_ , vol. 13, no. 6, pp. 1006–1011, 2009.
  * [15] H. Zhang, W. Yang, H. Yu, H. Zhang, and G.-S. Xia, “Detecting Power Lines in UAV Images with Convolutional Features and Structured Constraints,” _Remote Sensing_ , vol. 11, no. 11, p. 1342, 2019.
  * [16] I. Goodfellow, Y. Bengio, and A. Courville, “Part II: Modern Practical Deep Networks,” in _Deep learning_ , Cambridge, MA: MIT Press, 2017.
  * [17] J. V. D. Putten, F. V. D. Sommen, M. Struyvenberg, J. D. Groof, W. Curvers, E. Schoon, J. J. Bergman, and P. H. N. D. With, “Tissue segmentation in volumetric laser endomicroscopy data using FusionNet and a domain-specific loss function,” _Medical Imaging 2019: Image Processing_ , 2019.
  * [18] D. Li, J. Wu, Y. He, X. Yao, W. Yuan, D. Chen, H.-C. Park, S. Yu, J. L. Prince, and X. Li, “Parallel deep neural networks for endoscopic OCT image segmentation,” _Biomedical Optics Express_ , vol. 10, no. 3, p. 1126, Jul. 2019.
  * [19] Y. L. Yong, L. K. Tan, R. A. Mclaughlin, K. H. Chee, and Y. M. Liew, “Linear-regression convolutional neural network for fully automated coronary lumen segmentation in intravascular optical coherence tomography,” _Journal of Biomedical Optics_ , vol. 22, no. 12, p. 1, 2017.
  * [20] J. V. D. Putten, M. Struyvenberg, J. D. Groof, T. Scheeve, W. Curvers, E. Schoon, J. J. Bergman, P. H. D. With, and F. V. D. Sommen, “Deep principal dimension encoding for the classification of early neoplasia in Barretts Esophagus with volumetric laser endomicroscopy,” _Computerized Medical Imaging and Graphics_ , vol. 80, p. 101701, 2020.
  * [21] Y. Zeng, S. Xu, W. C. Chapman, S. Li, Z. Alipour, H. Abdelal, D. Chatterjee, M. Mutch, and Q. Zhu, “Real-time colorectal cancer diagnosis using PR-OCT with deep learning,” _Theranostics_ , vol. 10, no. 6, pp. 2587–2596, 2020.
  * [22] K.-K. Maninis, J. Pont-Tuset, P. Arbelaez, and L. V. Gool, “Convolutional Oriented Boundaries: From Image Segmentation to High-Level Tasks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 40, no. 4, pp. 819–833, 2018.
  * [23] M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision, _2nd ed. Pacific Grove_ , CA: PWS, 1999.
  * [24] W. Shen, X. Wang, Y. Wang, X. Bai, and Z. Zhang, “DeepContour: A deep convolutional feature learned by positive-sharing loss for contour detection,” _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015.
  * [25] G. Bertasius, J. Shi, and L. Torresani, “DeepEdge: A multi-scale bifurcated deep network for top-down contour detection,” _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015.
  * [26] J. Yang, B. Price, S. Cohen, H. Lee, and M.-H. Yang, “Object Contour Detection with a Fully Convolutional Encoder-Decoder Network,” _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * [27] J. L. Crassidis, F. L. Markley, and Y. Cheng, “Survey of Nonlinear Attitude Estimation Methods,” _Journal of Guidance, Control, and Dynamics_ , vol. 30, no. 1, pp. 12–28, 2007.
  * [28] A. Philipp and S. Behnke, “Robust sensor fusion for robot attitude estimation.” _2014 IEEE-RAS International Conference on Humanoid Robots_ , pp. 218-224., 2014.
  * [29] R. Mahony, T. Hamel and J.M. Pflimlin, “Complementary filter design on the special orthogonal group SO (3)”. _In Proceedings of the 44th IEEE Conference on Decision and Control_ , pp. 1477-1484, 2005.
  * [30] J. Justa, V. Šmídl, and A. Hamáček, “Fast AHRS Filter for Accelerometer, Magnetometer, and Gyroscope Combination with Separated Sensor Corrections,” _Sensors_ , vol. 20, no. 14, p. 3824, 2020.
  * [31] J. Wu, Z. Zhou, H. Fourati, and Y. Cheng, “A Super Fast Attitude Determination Algorithm for Consumer-Level Accelerometer and Magnetometer,” _IEEE Transactions on Consumer Electronics_ , vol. 64, no. 3, pp. 375–381, 2018.
  * [32] D. Gebre-Egziabher, R. Hayward, and J. Powell, “Design Of Multi-sensor Attitude Determination Systems,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 40, no. 2, pp. 627–649, 2004.
  * [33] Y. S. Suh, “Simple-Structured Quaternion Estimator Separating Inertial and Magnetic Sensor Effects,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 55, no. 6, pp. 2698–2706, 2019.
  * [34] S. O. H. Madgwick, A. J. L. Harrison, and R. Vaidyanathan, “Estimation of IMU and MARG orientation using a gradient descent algorithm,” z _2011 IEEE International Conference on Rehabilitation Robotics_ , 2011.
  * [35] J. Yu and R. Ramamoorthi, “Learning Video Stabilization Using Optical Flow,” _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [36] M. Grundmann, V. Kwatra, and I. Essa, “Auto-directed video stabilization with robust L1 optimal camera paths,” _2011 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2011.
  * [37] M. Wang, G.-Y. Yang, J.-K. Lin, S.-H. Zhang, A. Shamir, S.-P. Lu, and S.-M. Hu, “Deep Online Video Stabilization With Multi-Grid Warping Transformation Learning,” _IEEE Transactions on Image Processing_ , vol. 28, no. 5, pp. 2283–2292, 2019.
  * [38] S. Liu, L. Yuan, P. Tan, and J. Sun, “SteadyFlow: Spatially Smooth Optical Flow for Video Stabilization,” _2014 IEEE Conference on Computer Vision and Pattern Recognition_ , 2014.
  * [39] D. Sun, S. Roth, and M. J. Black, “Secrets of optical flow estimation and their principles,” _2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , 2010.
  * [40] E. Ilg, N. Mayer, T. Saikia, M. Keuper, A. Dosovitskiy, and T. Brox, “FlowNet 2.0: Evolution of Optical Flow Estimation with Deep Networks,” _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 2462-2470, 2017.
  * [41] F. Yu and V. Koltun, “Multi-scale context aggregation by dilated convolutions”. _Preprint at arXiv_ , https://arxiv.org/abs/1511.07122, 2015.
  * [42] M. D. Zeiler, G. W. Taylor, and R. Fergus, “Adaptive deconvolutional networks for mid and high level feature learning”. _IEEE International Conference on Computer Vision (ICCV)_ , pp. 2018–2025, 2011.
  * [43] J. Long, E. Shelhamer, and T. Darrell, “Fully convolutional networks for semantic segmentation”. _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015.
  * [44] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015.
  * [45] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _Int. Conf. Learn. Represent._ , 2015.
  * [46] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition”. _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770-778. 2016.
  * [47] O. C. Mora, P. Zanne, L. Zorn, F. Nageotte, N. Zulina, S. Gravelyn, P. Montgomery, M. de Mathelin, B. Dallemagne, and M. J. Gora, “Steerable OCT catheter for real-time assistance during teleoperated endoscopic treatment of colorectal cancer,” _Biomedical Optics Express_ , vol. 11, no. 3, p. 1231, 2020.
  * [48] L. Zorn, F. Nageotte, P. Zanne, A. Legner , B. Dallemagne, J. Marescaux, and M. de Mathelin “A Novel Telemanipulated Robotic Assistant for Surgical Endoscopy: Preclinical Application to ESD,” _IEEE Transactions on Biomedical Engineering_ , p. 797-808, Vol. 65, 4, 2018.
  * [49] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” _ICMLWorkshop on Deep Learning for Audio, Speech, and Language Processing (WDLASL 2013)_ , vol. 28., 2013.
  * [50] T. Wang, et al., “Heartbeat OCT: in vivo intravascular megahertz-optical coherence tomography”, _Biomedical optics express_ , 6(12), pp.5021-5032. 2015.
  * [51] M.J. Gora, et al., “Tethered capsule endomicroscopy enables less invasive imaging of gastrointestinal tract microstructure”. _Nature medicine_ , 19(2), pp.238-240. 2013.
  * [52] S.W. Lee, et al., “Quantification of airway thickness changes in smoke-inhalation injury using in-vivo 3-D endoscopic frequency-domain optical coherence tomography”. _Biomedical optics express_ , 2(2), pp.243-254, 2011.
  * [53] A. Paszke, et al., “Automatic differentiation in pytorch”. In NIPS-W, 2017.
  * [54] S. Ioffe and C. Szegedy, “ Batch normalization: Accelerating deep network training by reducing internal covariate shift”. _In International Conference on Machine Learning (ICML)_ , 2015.
  * [55] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization.” _2015 International Conference on Learning Representations_ , 2015.
  * [56] Y. Bengio, “Global Optimization Strategies,” in _Learning deep architectures for AI_ , pp.99 - 104., 2009.
  * [57] K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.” _In IEEE International Conference on Computer Vision (ICCV)_ , 2015.
  * [58] G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov, “Improving neural networks by preventing coadaptation of feature detectors.” arXiv:1207.0580, 2012.


# Standard Model prediction for paramagnetic EDMs

Yohei Ema<EMAIL_ADDRESS>William I. Fine Theoretical Physics Institute,
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA  Ting Gao<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  Maxim Pospelov
<EMAIL_ADDRESS>William I. Fine Theoretical Physics Institute, School of
Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA

###### Abstract

Standard Model $CP$ violation associated with the phase of the Cabibbo-
Kobayashi-Maskawa quark mixing matrix is known to give small answers for the
EDM observables. Moreover, predictions for the EDMs of neutrons and
diamagnetic atoms suffer from considerable uncertainties. We point out that
the $CP$-violating observables associated with the electron spin (paramagnetic
EDMs) are dominated by the combination of the electroweak penguin diagrams and
$\Delta I=1/2$ weak transitions in the baryon sector, and are calculable
within chiral perturbation theory. The predicted size of the semileptonic
operator $C_{S}$ is $-7\times 10^{-16}$ which corresponds to the equivalent
electron EDM $d_{e}^{\rm eq}=-1.0\times 10^{-35}e\,{\rm cm}$. While still far
from the current observational limits, this result is three orders of
magnitude larger than previously believed.

††preprint: UMN-TH-4115/22, FTPI-MINN-22-06

Introduction — The searches for EDMs of elementary particles Graner:2016ses ;
Cairncross:2017fip ; Andreev:2018ayy ; nEDM:2020crw  represent an important
way of probing the TeV scale new physics Ginges:2003qt ; Pospelov:2005pr ;
Engel:2013lsa . Recent breakthrough sensitivity to $CP$ violation connected to
electron spin (that we will refer to as “paramagnetic EDMs”) Andreev:2018ayy
established a new limit on the linear combination of the electron EDM $d_{e}$
and semileptonic nucleon-electron $\bar{N}N\bar{e}i\gamma_{5}e$ operators,
commonly parametrized by a $C_{S}$ coefficient. Given rapid progress of the
last decade, as well as some additional hopes for increased accuracy (see e.g.
Vutha:2017pej ; Vutha:2018tsz ) makes one to revisit the Standard Model (SM)
sources of $CP$ violation, and the expected size of the paramagnetic EDMs in
the SM.

SM has two sources of $CP$-violation. First source, undetected thus far,
corresponds to the non-perturbative effects parametrized by the QCD vacuum
angle $\theta$. Recently it has been shown Flambaum:2019ejc  that paramagnetic
EDMs are dominated by the two-photon exchange mechanism, and the leading
chiral behavior of the hadronic part of the diagram is given by the
$t$-channel exchange by $\pi^{0},\eta$. $CP$ violation due to $\theta$ comes
through the $\pi^{0}(\eta)\bar{N}N$ coupling. The result, in combination with
the experimental bound Andreev:2018ayy , sets the independent limit on
$|\theta|<3\times 10^{-8}$, which is still subdominant to the limit provided
by $d_{n}(\theta)$.

The second source of the SM $CP$-violation is the celebrated Kobayashi-Maskawa
(KM) phase $\delta_{\rm KM}$ Kobayashi:1973fv , which is now observed to
rather good accuracy in a plethora of flavor transitions in $B$ and $K$
mesons. Observations are often matched by rather precise theoretical
predictions, starting from Bigi:1981qs . The predictions of EDM-like
observables induced by $\delta_{\rm KM}$ thus far can be summarized by two
adjectives: small and uncertain. The suppression comes from the necessity to
involve at least two $W$-bosons and multiple loops Shabalin:1978rs ;
Khriplovich:1985jr ; Pospelov:1991zt  involving all three generations of
quarks. As a result, short distance contributions to quark EDMs do not exceed
$10^{-33}\,e\,{\rm cm}$ level Czarnecki:1997bu . At the same time, it is clear
that long-distance nonperturbative contributions, typically described as a
combination of two transitions changing strangeness by one unit, $\Delta S=\pm
1$, dominate $d_{n}$ and nucleon-nucleon forces Khriplovich:1981ca ;
Gavela:1981sk ; Flambaum:1984fb ; Donoghue:1987dd ; McKellar:1987tf . More
recent estimate Seng:2014lea  places $d_{n}$ in the ballpark of ${\rm
few}\times 10^{-32}\,e\,{\rm cm}$ with a wide order-of-magnitude expected
range. It is fair to say that magnitudes of $d_{n}$ and nucleon-nucleon forces
(that feeds into the nuclear-spin-dependent atomic EDMs) cannot be accurately
predicted at this point.

What is the size of paramagnetic EDMs induced by $\delta_{\rm KM}$? Recent
estimates of $d_{e}$ Yamaguchi:2020eub  (dominated again by long-distance
effects) converge at the tiniest value of $\sim 6\times 10^{-40}\,e\,{\rm
cm}$, presumably with considerable uncertainties corresponding to hadronic
modelling of quark loops. This result is subdominant to the $C_{S}$ estimate
due to the two-photon exchange mechanism in combination with $\Delta S=\pm 1$
transitions Pospelov:2013sca , that corresponds to equivalent $d_{e}$ of $\sim
10^{-38}\,e\,{\rm cm}$. To introduce useful notations, this is ${\rm
EW}^{2}{\rm EM}^{2}$ order effect, where EW/EM stands for
electroweak/electromagnetic.

In this Letter we demonstrate that the dominant contribution to paramagnetic
EDMs associated with the KM $CP$-violation is given by the semileptonic
$C_{S}$ induced in ${\rm EW}^{3}$ order. It has an unambiguous answer in the
flavor-$SU(3)$ chiral limit, and is calculable to $\sim 30\%$ accuracy that
can be further improved. Remarkably, the result reaches the level of $\sim
10^{-35}\,e\,{\rm cm}$ in terms of the $d_{e}$ equivalent, which is three
orders of magnitude larger than previously believed Pospelov:2013sca .

Our starting point is the expression for the equivalent $d_{e}$ that follows
from atomic/molecular theory, and defines the linear combination of two Wilson
coefficients constrained by the most precise paramagnetic EDM measurements
performed with ThO molecule:

$d_{e}^{\rm equiv}=d_{e}+C_{S}\times 1.5\times 10^{-20}\,e\,{\rm cm},$ (1)

where $e$ is the positron charge. Current experimental limit Andreev:2018ayy
stands as $|d_{e}^{\rm equiv}|<1.1\times 10^{-29}\,e\,{\rm cm}$. As per
convention, $C_{S}$ is defined with the Fermi constant factored out, and
$\gamma_{5}$ corresponds to the $\frac{1}{2}\gamma_{\mu}(1-\gamma_{5})$
definition of the left-handed current:

${\cal
L}_{eN}=C_{S}\frac{G_{F}}{\sqrt{2}}(\bar{e}i\gamma_{5}e)(\bar{p}p+\bar{n}n).$
(2)

Our goal is to calculate $C_{S}(\delta_{\rm KM})$.

Leading chiral order $C_{S}$ calculation — Because of the conservation of the
electron chirality in the SM, it is clear that $C_{S}\propto m_{e}$. This in
turn rules out single photon exchange (EM penguin) as origin of
$m_{e}\bar{e}i\gamma_{5}e$, and one would need either a two-photon mechanism
Pospelov:2013sca ; Flambaum:2019ejc  or the EW penguin $Z$-boson
exchange/$W$-box diagram. The most crucial property of EW penguins is that
although they are formally of the second order in weak interactions, their
size is enhanced by the heavy top, so that the result scales as
$G_{F}^{2}m_{t}^{2}$. EW penguins111As is well known, EW penguins must also
include $W$-box diagrams, and we include both. induce
$B_{s,d}\to\mu^{+}\mu^{-}$ decays, and dominate the dispersive part of
$K_{L}\to\mu^{+}\mu^{-}$ amplitude. Dropping the vector part of the lepton
current (as not leading to $m_{e}\bar{e}i\gamma_{5}e$), and integrating out
heavy $W,Z,t$ particles, one can concisely write down the semileptonic
operator as

$\displaystyle{\cal L}_{\rm EWP}={\cal P}_{\rm
EW}\times\bar{e}\gamma_{\mu}\gamma_{5}e\times\bar{s}\gamma^{\mu}(1-\gamma_{5})d+(h.c.),$
(3)

where

$\displaystyle{\cal P}_{\rm EW}=\frac{G_{F}}{\sqrt{2}}\times
V_{ts}^{*}V_{td}\times\frac{\alpha_{\rm
EM}(m_{Z})}{4\pi\sin^{2}\theta_{W}}I(x_{t}),$ (4)

and the loop function is given by Inami:1980fz

$I(x_{t})=\frac{3}{4}\left(\frac{x_{t}}{x_{t}-1}\right)^{2}\log
x_{t}+\frac{1}{4}x_{t}-\frac{3}{4}\frac{x_{t}}{x_{t}-1},~{}x_{t}=\frac{m_{t}^{2}}{m_{W}^{2}}.$
(5)

These results are well established, and unlike the case of four-quark
operators, the subsequent QCD evolution of (3) introduces only small
corrections (see e.g. Buras:2012ru ).

Figure 1: EW3 order diagram that dominates in the chiral limit. The top vertex
is the $CP$-odd, $P$-even $K_{S}\bar{e}i\gamma_{5}e$ generated in EW2 order,
and the bottom vertex is $CP$-even, $P$-odd $K_{S}\bar{N}N$ coupling generated
at EW1 order.

The most convenient representation of the CKM matrix is when $\delta_{\rm KM}$
enters mostly in $V_{td}$. It enters the imaginary part of ${\cal P}_{\rm EW}$
and couples the axial vector current of leptons to the
$\bar{s}\gamma_{\mu}(1-\gamma_{5})d-\bar{d}\gamma_{\mu}(1-\gamma_{5})s$ quark
current. This current can create/annihilate $CP$-even combination of the
neutral kaons that (in neglection of small $\epsilon_{K}$) can be identified
with $K_{S}$ field. Same operator in the muon channel induces
$K_{S}\to\mu^{+}\mu^{-}$ meson decay Isidori:2003ts ; DAmbrosio:2017klp .
Therefore, at the hadronic level, retaining only the lightest open-strangeness
fields, we can re-write (3) as
$\partial_{\mu}K\times\bar{e}\gamma^{\mu}\gamma_{5}e$, and upon application of
the equation of motion for electrons arrive to

${\cal L}_{Kee}=2\sqrt{2}f_{0}m_{e}\bar{e}i\gamma_{5}e\left(K_{S}\times{\rm
Im}{\cal P}_{\rm EW}+K_{L}\times{\rm Re}{\cal P}_{\rm EW}\right).$ (6)

In this expression, $f_{0}$ is the meson coupling constant, that in the
$SU(3)$ symmetric limit is equal to $\simeq 134$ MeV, and we follow Ref.
Bijnens:1985kj  conventions. Subsequent $m_{s}$-dependent corrections
renormalize this coupling to $f_{0}\to f_{K}\simeq 160$ MeV. While other
$s$-quark containing resonances may also contribute, the neutral kaon
exchange, Fig. 1, will give the only $m_{s}^{-1}$-enhanced contribution in the
chiral limit.

We now need to find out how the neutral kaons couple to the nucleon scalar
densities, $\bar{p}p$ and $\bar{n}n$ that occur due to $\Delta S=\pm 1$
transitions in the EW1 order. Instead of attempting such calculation from
first principles (see e.g. Shifman:1975tn ) we will use flavor $SU(3)$
relations and connect this coupling to the $s$-wave amplitudes of hadronic
decays of strange hyperons, following Bijnens:1985kj . It is well known that
empirical $\Delta I=1/2$ rule holds for hyperon decays, and the leading order
$SU(3)$ relations fit $s$-wave amplitudes with $O(10\%)$ accuracy. It is
strongly suspected that these amplitudes are indeed induced by strong penguins
(SP), although this assumption is not crucial for us. With that, one can write
down the two types of couplings consistent with $(8_{L},1_{R})$ transformation
properties:

${\cal L}_{\rm SP}=-a{\rm Tr}(\bar{B}\\{\xi^{\dagger}h\xi,B\\})-b{\rm
Tr}(\bar{B}[\xi^{\dagger}h\xi,B])+(h.c.),$ (7)

where $h_{ij}=\delta_{i2}\delta_{j3}$. In this expression, $B$ is the baryon
octet matrix, and $\xi$ is the exponential of the meson octet $M$,
$\xi=\exp[iMf_{0}^{-1}]$. Assuming $a$ and $b$ to be real, and taking
$f_{0}=f_{\pi}$, they are fit to be

$a=0.56G_{F}f_{\pi}\times[m_{\pi^{+}}]^{2};~{}~{}b=-1.42G_{F}f_{\pi}\times[m_{\pi^{+}}]^{2}.$
(8)

Brackets over $m_{\pi^{+}}$ indicate that these are numerical values taken,
$139.5$ MeV, rather than $m_{u}+m_{d}$-proportional theoretical quantity
$m_{\pi}$. In the assumption of $a$ and $b$ being real, only the $K_{S}$ meson
couples to nucleons, $2^{1/2}f_{0}^{-1}((b-a)\bar{p}p+2b\bar{n}n)K_{S}$, which
will provide the dominant contribution. This type of coupling breaks $P$ but
respects $CP$ symmetry. Restoring the CKM factors, one can also include much
subdominant coupling to $K_{L}$ so that we have:

$\displaystyle{\cal
L}_{KNN}\simeq-\frac{\sqrt{2}\,G_{F}\times[m_{\pi^{+}}]^{2}f_{\pi}}{|V_{ud}V_{us}|f_{0}}\times
2.84(0.7\bar{p}p+\bar{n}n)$ (9) $\displaystyle\times\left({\rm
Re}(V_{ud}^{*}V_{us})K_{S}+{\rm Im}(V_{ud}^{*}V_{us})K_{L}\right).$

At the last step, we integrate out the $K$ mesons as shown in Fig. 1. Adopting
it for a nucleus containing $A=Z+N$ nucleons, one arrives to a straightforward
prediction for the $\delta_{\rm KM}$-induced size of the electron-nucleon
interaction:

$C_{S}\simeq-{\cal
J}\times\frac{N+0.7Z}{A}\times\frac{13[m_{\pi^{+}}]^{2}f_{\pi}m_{e}G_{F}}{m_{K}^{2}}\times\frac{\alpha_{\rm
EM}I(x_{t})}{\pi\sin\theta_{W}^{2}},$ (10)

where ${\cal J}$ is the rephasing invariant combination of the CKM angles,

${\cal J}={\rm Im}(V_{ts}^{*}V_{td}V_{ud}^{*}V_{us})\simeq 3.1\times 10^{-5},$
(11)

that carries about $\sim 6\%$ uncertainty. Notice that the $f_{0}$ factor in
the numerator of (6) cancels against $f_{0}$ in the denominator of (9), and
this cancellation would persist even one changes $f_{0}$ for $f_{K}$.

The overall scaling of this formula in the chiral limit and at large $x_{t}$
is

$G_{F}C_{S}\propto{\cal J}G_{F}^{3}m_{t}^{2}m_{e}m_{s}^{-1}\Lambda_{\rm
hadr}^{2}.$ (12)

where $\Lambda_{\rm hadr}$ is a typical hadronic energy/momentum scale. Notice
that this is far more singular behavior with $m_{q}$ of a light quark than
that arising in the chiral-loop-induced expressions for $d_{n}$. Also notice
that the $K_{S}$ exchange dominates for any conventional parametrization of
the CKM matrix, and the role of $K_{L}$ exchange is to add small pieces of the
amplitude that take ${\rm Re}(V_{ud}V_{us}^{*}){\rm Im}(V_{ts}V_{td}^{*})$,
arising from $K_{S}$ exchange, to full ${\cal J}$. Substituting all SM
parameters, we obtain the following leading order result:

$C_{S}({\rm LO})\simeq-5\times 10^{-16}.$ (13)

In order to estimate accuracy of the LO $\sim O(m_{s}^{-1})$ result, one could
try to evaluate the NLO corrections in the expansion over small $m_{s}$. These
corrections can be divided into two groups: A. corrections to the $K\bar{N}N$
vertex at $m_{s}\log m_{s}$ order, B. diagrams that do not reduce to the
$t$-channel $K$-meson exchange. Type A corrections involve essentially same
diagrams as those appearing in the corresponding corrections to the $s$-wave
hyperon decays Bijnens:1985kj ; Jenkins:1991bt . The analysis of Ref.
Jenkins:1991bt  showed that when the loop corrections are included with the
tree-level $a$ and $b$ parameters and the total theoretical result is fit to
experimental data, one notices that the tree level values for $a$ and $b$ come
out smaller than in (8), while total result is rather close to the tree-level
fit for $a,b$. This comes mostly from the renormalization of the meson and
baryon wave functions. The lesson from this is that the corrections of type A
for $KNN$ weak coupling are expected to mirror results of Ref. Jenkins:1991bt
for $s$-wave amplitudes, and therefore would not deviate substantially from
Eq. (9).

Figure 2: The baryon pole diagrams that contribute to $C_{S}$ at the NLO level
in the chiral limit. The left vertex is the nucleon-hyperon mixing induced by
Eq. (7), while the top vertex is induced by Eq. (3). The vertices without
black dots are the strong interaction with the coupling constants $D$ and $F$.
The diagrams with the nulceon-hyperon mixing on the right side give the same
amount of contribution.

We then estimate type B corrections. It turns out that they parametrically
dominate over other types of corrections, as the baryon pole diagrams, Fig. 2,
contribute. The $m_{s}$ scaling of these corrections is set by the ratio of
the loop integral, proportional to $m_{K}$ (at $m_{K}^{2}\gg m_{\pi}^{2}$
limit), divided by mass splitting $\Delta m_{B}$ in the baryon octet, e.g.
$m_{\Lambda}-m_{n}$. This quantity scales as $m_{s}^{-1/2}$ and therefore
these baryon pole diagrams dominate the NLO contributions in the chiral limit.
They are fully calculable (i.e. do not depend on unknown counterterms), and
the result for these corrections are:

$\displaystyle\frac{C_{S,NLO}(p)}{C_{S,LO}(p)}=\frac{m_{K}^{3}(0.77D^{2}+2.7DF-2.3F^{2})}{24\pi
f_{0}^{2}(m_{\Sigma^{+}}-m_{p})}~{}~{}~{}$ (14)
$\displaystyle\frac{C_{S,NLO}(n)}{C_{S,LO}(n)}=\frac{m_{K}^{3}}{24\pi
f_{0}^{2}}\left(\frac{(a/b+3)}{2\sqrt{6}(m_{\Lambda}-m_{n})}~{}~{}~{}\right.$
(15) $\displaystyle\times(-0.44D^{2}+3.2DF+1.3F^{2})$
$\displaystyle+\frac{a/b-1}{2\sqrt{2}(m_{\Sigma^{0}}-m_{n})}\left.(-0.53D^{2}-1.9DF+1.6F^{2})\right).$

It has been obtained using heavy baryon chiral perturbation theory, and $D,F$
are the coupling constants characterizing the strength of the
$SU(3)$-invariant baryon-meson strong interaction, with $F=0.46,\,D=0.8$
typically used Bijnens:1985kj . Since the dominant contribution comes from
loops with $K-\pi$ transition, it is appropriate to take $f_{0}^{2}\simeq
f_{\pi}f_{K}$. Using these numbers, we discover that NLO corrections interfere
constructively with LO, and give 30% correction for the proton, and 40% for
the neutron, correspondingly. Combining LO and NLO, we arrive at our final
result,

$\displaystyle C_{S}({\rm LO+NLO})\simeq-6.9\times 10^{-16}$
$\displaystyle~{}\Longrightarrow~{}d_{e}^{\rm equiv}\simeq-1.0\times
10^{-35}\,e\,{\rm cm}.$ (16)

The size of the NLO corrections also allows us to estimate the accuracy of
this computation as $O(30\%)$.

As stated in the introduction, this result is much larger than previously
believed, and exceeds any contributions of $d_{e}$ into $d_{e}^{\rm equiv}$ by
at least four orders of magnitude. The enhancement of $C_{S}$ at EW3 order
compared to EW2EM2 can be roughly ascribed to $\alpha_{W}/\alpha_{\rm
EM}^{2}\sim O(10^{3})$. We note that, although translating $C_{S}$ to
$d_{e}^{\mathrm{equiv}}$ depends on atoms/molecules that one considers (ThO
above), this dependence is mild and $d_{e}^{\mathrm{equiv}}$ is within the
same ballpark if we instead consider, _e.g._ , Tl or YbF Pospelov:2013sca .

Stepping away from chiral expansion, one can formulate the necessary hadronic
matrix element that will be required to generate $C_{S}$ in combination with
the dominant Im${\cal P}_{\rm EW}$ channel of Eq. (3). The corresponding
$d$-to-$s$ transitions need to be taken in the first order, EW1, that break
$P$ and $C$ separately but conserve $CP$.

$\displaystyle\langle
N|i(\bar{s}\gamma_{\mu}(1-\gamma_{5})d-\bar{d}\gamma_{\mu}(1-\gamma_{5})s)|N\rangle_{\rm
EW^{1}}$ (17)
$\displaystyle=\frac{f_{S}}{m_{N}}iq_{\mu}\bar{N}N+\frac{f_{T}}{m_{N}}q_{\nu}\bar{N}\sigma_{\mu\nu}\gamma_{5}N.$

In this formula, $q_{\mu}$ stands for the momentum transfer. It turns out that
there are only two form factors on the r.h.s. of this expression that have the
same $CP$ properties as the left-hand side. Moreover, $f_{T}$ in combination
with Im${\cal P}_{\rm EW}$ leads to $CP$-odd $P$-even interactions that do not
induce EDMs. Therefore only $f_{S}$ form factor (that sometimes is called
induced scalar) at $q^{2}\to 0$ is relevant. We have provided first two terms
in the chiral expansion of $f_{S}$ for neutrons and protons, so that
effectively $f_{S}\propto a(b)\times m_{N}m_{K}^{-2}+...$. While we use chiral
perturbation theory, in principle, calculation of $f_{S}$ can be attempted
using lattice QCD methods.

Finally, we note that other semi-leptonic operators such as
$\bar{e}e\bar{N}i\gamma_{5}N$ that lead to nuclear-spind-dependent effects are
not generated the same way at EW3 order and therefore will be suppressed
compared to (10).

Conclusions — We have shown that $\delta_{\rm KM}$ induces the $CP$-odd
electron nucleon interaction at the level much larger than previous estimates
Pospelov:2013sca . The main mechanism is not a two-photon exchange, EW2EM2,
between electron and the nucleus, but the combination of weak non-leptonic EW1
transition with the semileptonic EW2 electroweak penguin. Although the result
is still small, it is not unthinkable that the progress in sensitivty to
paramagnetic EDMs may reach the level of $d_{e}^{\rm equiv}$ in the future.
Indeed, some novel proposals Vutha:2017pej  envision that statistical
sensitivity to paramagnetic EDMs can be brought down to $d_{e}\sim
O(10^{-35}\mathchar 45\relax\mathchar 45\relax 10^{-37})\,e\,\mathrm{cm}$.

It is not surprising that the $C_{S}$ operator can be predicted, at least in
the $SU(3)$ chiral expansion, rather precisely. This puts in clear distinction
with $d_{n}(\delta_{\rm KM})$ estimate that carries an order of magnitude
uncertainty with unclear prospects for improvement. In contrast, the only
significant source of uncertainty in $C_{S}$ is in the induced scalar form
factor (17) that can be improved in the future with the use of lattice QCD
methods.

Even if one takes chiral $SU(3)$ expansion sceptically, it is clear that
unique $m_{s}^{-1}$ (LO) and $m_{s}^{-1/2}$ (NLO) contributions to $C_{S}$
identified in our work would not be cancelled - unless completely accidentally
- by other contributions, mirroring a similar argument of Crewther:1979pi
made for $d_{n}(\theta)$. Therefore, $10^{-35}\,e\,{\rm cm}$ should be adopted
as the robust $\delta_{\rm KM}$-induced SM benchmark value for all experiments
attempting the search of $d_{e}$ using electron spins in heavy atoms and
molecules. It also allows for establishing the maximum sensitivity to CP-
violating New Physics via $d_{e}$. Taking a one-loop perturbative scaling,
$d_{e}\propto(\alpha/\pi)m_{e}\Lambda_{\rm NP}^{-2}$, and equating it to
$d_{e}^{\rm equiv}(\delta_{\rm KM})$ one arrives at the maximum scale that is
possible to probe with paramagnetic EDMs: $\Lambda_{\rm NP}^{\rm max}\sim
5\times 10^{7}$ GeV.

###### Acknowledgements.

Acknowledgments — Y.E. and M.P. are supported in part by U.S. Department of
Energy Grant No. desc0011842. The Feynman diagrams in this paper are drawn
with TikZ-Feynman Ellis:2016jkw .

## References

  * (1) B. Graner, Y. Chen, E. G. Lindahl, and B. R. Heckel, “Reduced Limit on the Permanent Electric Dipole Moment of Hg199,” Phys. Rev. Lett. 116 no. 16, (2016) 161601, arXiv:1601.04339 [physics.atom-ph]. [Erratum: Phys. Rev. Lett.119,no.11,119901(2017)].
  * (2) W. B. Cairncross, D. N. Gresh, M. Grau, K. C. Cossel, T. S. Roussy, Y. Ni, Y. Zhou, J. Ye, and E. A. Cornell, “Precision Measurement of the Electron?s Electric Dipole Moment Using Trapped Molecular Ions,” Phys. Rev. Lett. 119 no. 15, (2017) 153001, arXiv:1704.07928 [physics.atom-ph].
  * (3) ACME Collaboration, V. Andreev et al., “Improved limit on the electric dipole moment of the electron,” Nature 562 no. 7727, (2018) 355–360.
  * (4) nEDM Collaboration, C. Abel et al., “Measurement of the permanent electric dipole moment of the neutron,” Phys. Rev. Lett. 124 no. 8, (2020) 081803, arXiv:2001.11966 [hep-ex].
  * (5) J. S. M. Ginges and V. V. Flambaum, “Violations of fundamental symmetries in atoms and tests of unification theories of elementary particles,” Phys. Rept. 397 (2004) 63–154, arXiv:physics/0309054 [physics].
  * (6) M. Pospelov and A. Ritz, “Electric dipole moments as probes of new physics,” Annals Phys. 318 (2005) 119–169, arXiv:hep-ph/0504231 [hep-ph].
  * (7) J. Engel, M. J. Ramsey-Musolf, and U. van Kolck, “Electric Dipole Moments of Nucleons, Nuclei, and Atoms: The Standard Model and Beyond,” Prog. Part. Nucl. Phys. 71 (2013) 21–74, arXiv:1303.2371 [nucl-th].
  * (8) A. C. Vutha, M. Horbatsch, and E. A. Hessels, “Oriented polar molecules in a solid inert-gas matrix: a proposed method for measuring the electric dipole moment of the electron,” Atoms 6 no. 1, (2018) 3, arXiv:1710.08785 [physics.atom-ph].
  * (9) A. C. Vutha, M. Horbatsch, and E. A. Hessels, “Orientation-dependent hyperfine structure of polar molecules in a rare-gas matrix: A scheme for measuring the electron electric dipole moment,” Phys. Rev. A 98 no. 3, (2018) 032513, arXiv:1806.06774 [physics.atom-ph].
  * (10) V. V. Flambaum, M. Pospelov, A. Ritz, and Y. V. Stadnik, “Sensitivity of EDM experiments in paramagnetic atoms and molecules to hadronic CP violation,” Phys. Rev. D 102 no. 3, (2020) 035001, arXiv:1912.13129 [hep-ph].
  * (11) M. Kobayashi and T. Maskawa, “CP Violation in the Renormalizable Theory of Weak Interaction,” Prog. Theor. Phys. 49 (1973) 652–657.
  * (12) I. I. Y. Bigi and A. I. Sanda, “Notes on the Observability of CP Violations in B Decays,” Nucl. Phys. B 193 (1981) 85–108.
  * (13) E. P. Shabalin, “Electric Dipole Moment of Quark in a Gauge Theory with Left-Handed Currents,” Sov. J. Nucl. Phys. 28 (1978) 75.
  * (14) I. B. Khriplovich, “Quark Electric Dipole Moment and Induced $\theta$ Term in the Kobayashi-Maskawa Model,” Phys. Lett. B 173 (1986) 193–196.
  * (15) M. E. Pospelov and I. B. Khriplovich, “Electric dipole moment of the W boson and the electron in the Kobayashi-Maskawa model,” Sov. J. Nucl. Phys. 53 (1991) 638–640.
  * (16) A. Czarnecki and B. Krause, “Neutron electric dipole moment in the standard model: Valence quark contributions,” Phys. Rev. Lett. 78 (1997) 4339–4342, arXiv:hep-ph/9704355.
  * (17) I. B. Khriplovich and A. R. Zhitnitsky, “What Is the Value of the Neutron Electric Dipole Moment in the Kobayashi-Maskawa Model?,” Phys. Lett. B 109 (1982) 490–492.
  * (18) M. B. Gavela, A. Le Yaouanc, L. Oliver, O. Pene, J. C. Raynal, and T. N. Pham, “CP Violation Induced by Penguin Diagrams and the Neutron Electric Dipole Moment,” Phys. Lett. B 109 (1982) 215–220.
  * (19) V. V. Flambaum, I. B. Khriplovich, and O. P. Sushkov, “On the Possibility to Study $P$ Odd and $T$ Odd Nuclear Forces in Atomic and Molecular Experiments,” Sov. Phys. JETP 60 (1984) 873.
  * (20) J. F. Donoghue, B. R. Holstein, and M. J. Musolf, “Electric Dipole Moments of Nuclei,” Phys. Lett. B 196 (1987) 196–202.
  * (21) B. H. J. McKellar, S. R. Choudhury, X.-G. He, and S. Pakvasa, “The Neutron Electric Dipole Moment in the Standard $K^{-}$m Model,” Phys. Lett. B 197 (1987) 556–560.
  * (22) C.-Y. Seng, “Reexamination of The Standard Model Nucleon Electric Dipole Moment,” Phys. Rev. C 91 no. 2, (2015) 025502, arXiv:1411.1476 [hep-ph].
  * (23) Y. Yamaguchi and N. Yamanaka, “Large long-distance contributions to the electric dipole moments of charged leptons in the standard model,” Phys. Rev. Lett. 125 (2020) 241802, arXiv:2003.08195 [hep-ph].
  * (24) M. Pospelov and A. Ritz, “CKM benchmarks for electron electric dipole moment experiments,” Phys. Rev. D89 no. 5, (2014) 056006, arXiv:1311.5537 [hep-ph].
  * (25) T. Inami and C. S. Lim, “Effects of Superheavy Quarks and Leptons in Low-Energy Weak Processes k(L) —$>$ mu anti-mu, K+ —$>$ pi+ Neutrino anti-neutrino and K0 $<$—$>$ anti-K0,” Prog. Theor. Phys. 65 (1981) 297. [Erratum: Prog.Theor.Phys. 65, 1772 (1981)].
  * (26) A. J. Buras, J. Girrbach, D. Guadagnoli, and G. Isidori, “On the Standard Model prediction for BR(Bs,d to mu+ mu-),” Eur. Phys. J. C 72 (2012) 2172, arXiv:1208.0934 [hep-ph].
  * (27) G. Isidori and R. Unterdorfer, “On the short distance constraints from K(L,S) —$>$ mu+ mu-,” JHEP 01 (2004) 009, arXiv:hep-ph/0311084.
  * (28) G. D’Ambrosio and T. Kitahara, “Direct $CP$ Violation in $K\to\mu^{+}\mu^{-}$,” Phys. Rev. Lett. 119 no. 20, (2017) 201802, arXiv:1707.06999 [hep-ph].
  * (29) J. Bijnens, H. Sonoda, and M. B. Wise, “On the Validity of Chiral Perturbation Theory for Weak Hyperon Decays,” Nucl. Phys. B 261 (1985) 185–198.
  * (30) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Light Quarks and the Origin of the Delta I = 1/2 Rule in the Nonleptonic Decays of Strange Particles,” Nucl. Phys. B 120 (1977) 316–324.
  * (31) E. E. Jenkins, “Hyperon nonleptonic decays in chiral perturbation theory,” Nucl. Phys. B 375 (1992) 561–581.
  * (32) R. J. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten, “Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics,” Phys. Lett. B 88 (1979) 123. [Erratum: Phys.Lett.B 91, 487 (1980)].
  * (33) J. Ellis, “TikZ-Feynman: Feynman diagrams with TikZ,” Comput. Phys. Commun. 210 (2017) 103–123, arXiv:1601.05437 [hep-ph].