In this section, we study the performance of the proposed model on two real-
world datasets.

##  Experimental Settings

###  Datasets

We evaluate all methods using two real-world city-wide datasets, Beijing (BJ)
and Shanghai (SH), from Meituan111https://www.meituan.com which is one of the
largest location-based service providers in China. We construct relationships
between POIs from user logs extracted from these two datasets. Following the
previous work [KDD15_itemrelation, WSDM18_itemrelation, CIKM20_DecGCN], we
generate the ground truth of competitive relationship for a POI pair if the
two POIs in the pair are ”viewed/clicked together” by users within a query
session, and we generate the ground truth of complementary relationships for a
POI pair if the two POIs in the pair are ”also viewed/clicked” by same users
across different query sessions. Note that user logs are often unavailable or
insufficient, and thus the way of generating ground-truth cannot be
generalized to cities where they are unavailable, or are insufficient to cover
all the POIs. Our proposed methods do not use user logs, i.e., the information
of generating the ground-truth relationships for evaluation, and are
applicable when the user logs are unavailable, or can be used to enrich the
POI relationships when the initial data sources are incomplete. Details of the
extracted graphs are shown in Table .

Table : Statistics of the datasets

0.9!  2*Dataset Taxonomy 2*#POIs 2*#Relational Edges #Non-leaf nodes
#Categories _Beijing_ 2*95 2*805 13,334 122,462 _Shanghai_ 10,090 112,848

###  Baselines and Evaluation Metrics

We compare our model with five types of baseline methods, including rule-based
methods (CAT, CAT-D), random walk based graph embedding methods (Deepwalk,
node2vec), vanilla GNN (GCN, GAT), heterogeneous GNN (HAN, HGT, R-GCN,
CompGCN), and the state-of-the-art relationship inference methods (DecGCN,
DeepR).

  * •

CAT and CAT-D. CAT determines the relation type by a threshold of POI category
distance on the category taxonomy, and CAT-D considers both POI geographical
distance and POI category distance thresholds. We search the thresholds that
achieve the best results on two datasets.

  * •

Deepwalk [Deepwalk] and node2vec [node2vec] learn node representations from
sequences generated by random walks on a graph. node2vec adopts biased random
walks, which is more generalized compared to the unbiased ones adopted in
Deepwalk.

  * •

GCN [GCN] and GAT [GAT] are two vanilla GNN models which produce
representations by aggregating from neighbor nodes. GCN treats neighbor nodes
equally while GAT considers different weights for them via attention
mechanism. These two methods do not model different relation types between
POIs.

  * •

HAN [WWW19_HAN] is a meta-path based heterogeneous GNN model, which performs
node-level and semantic-level attention to model multiple relation types.

  * •

HGT [WWW20_HGT] handles graph heterogeneity with relation-specific mutual
attention mechanism which models neighbor aggregation in different latent
spaces for each relation type.

  * •

R-GCN [RGCN] is an extension of GCN to better model relational graphs by
assigning different weight matrices for different relation types in neighbor
aggregation.

  * •

CompGCN [ICLR20_CompGCN] jointly learns the embeddings of nodes and relations
for heterogeneous graph and leverages multiple composition functions to update
the node representations.

  * •

DecGCN [CIKM20_DecGCN] is a state-of-the-art product relationship inference
method that decomposes a heterogeneous graph into sub-graphs, one for each
relation type. GNN is applied in each sub-graph followed by a co-attention
mechanism to extract supplementary information from other sub-graphs.

  * •

DeepR [KDD20_POI] is a state-of-the-art method that aims to infer the
competitiveness between POIs. It considers spatial features by splitting the
neighbors of a POI into different sectors based on the coordinates and
performs aggregation from each sector. To make it applicable to multiple
relationships, we extract a sub-graph for each relation type, and apply the
method to each sub-graph.

To evaluate the performance of different methods, we adopt the widely-used
_Micro-F1_ and _Macro-F1_ in previous studies [WSDM18_itemrelation, KDD20_POI]
to jointly consider precision and recall. A higher score indicates a better
performance.

###  Parameter Settings

In our experiments, we randomly sample 10% of the edges as validation data,
20% of the edges as test data, and different fractions of the remaining edges
are used as training data. Apart from the relational edges, 16,000 POI pairs
which do not have relationships are randomly sampled for the non-relation type
for testing. In our model, we set the threshold $d$ for determining spatial
neighbors to be $1.15km$, the scaling factor $\theta$ in the radial basis
function (RBF) kernel to be 2, the number of heads to be 4, the category
embedding size to be 128, and the number of negative samples $\omega$ to be 5.
For random walk based models, we set window size to 5, walk length to 30, and
walks per node to 20. For all the GNN based methods, we set the number of
layers to be 3. To achieve fair comparison, POI embedding size is set to be
128 for all the methods. We use Adam optimizer with learning rate of 0.001 and
batch size of 512 to train our model.

Table : Results on the two datasets in terms of Macro-F1 and Micro-F1 (with
best in bold and second-best underlined)

0.9!  Dataset Metric Train% CAT CAT-D Deepwalk node2vec GCN GAT HAN HGT R-GCN
CompGCN DecGCN DeepR PRIM 8*BJ 4*Macro-F1 40% 0.464 0.519 0.638 0.640 0.707
0.724 0.782 0.779 0.789 0.794 0.757 0.783 0.845 50% 0.464 0.519 0.691 0.692
0.737 0.748 0.811 0.814 0.814 0.832 0.801 0.815 0.870 60% 0.464 0.519 0.731
0.734 0.755 0.776 0.839 0.842 0.820 0.860 0.811 0.842 0.882 70% 0.464 0.519
0.757 0.761 0.770 0.795 0.857 0.857 0.828 0.870 0.823 0.861 0.895 4*Micro-F1
40% 0.559 0.579 0.707 0.710 0.729 0.753 0.817 0.813 0.808 0.827 0.805 0.820
0.879 50% 0.559 0.579 0.758 0.762 0.766 0.776 0.842 0.845 0.827 0.859 0.823
0.847 0.895 60% 0.559 0.579 0.783 0.784 0.780 0.804 0.867 0.869 0.832 0.882
0.826 0.871 0.907 70% 0.559 0.579 0.816 0.817 0.796 0.821 0.882 0.883 0.839
0.892 0.843 0.887 0.913 8*SH 4*Macro-F1 40% 0.443 0.509 0.652 0.655 0.673
0.692 0.788 0.790 0.777 0.776 0.765 0.786 0.822 50% 0.443 0.509 0.682 0.684
0.704 0.711 0.807 0.811 0.802 0.808 0.797 0.808 0.844 60% 0.443 0.509 0.698
0.702 0.713 0.729 0.834 0.838 0.828 0.837 0.817 0.839 0.861 70% 0.443 0.509
0.707 0.711 0.731 0.735 0.853 0.858 0.839 0.857 0.832 0.852 0.875 4*Micro-F1
40% 0.551 0.573 0.724 0.727 0.744 0.769 0.843 0.841 0.822 0.836 0.824 0.845
0.886 50% 0.551 0.573 0.766 0.768 0.775 0.788 0.861 0.864 0.837 0.860 0.846
0.863 0.896 60% 0.551 0.573 0.783 0.786 0.793 0.806 0.886 0.888 0.864 0.891
0.856 0.890 0.909 70% 0.551 0.573 0.796 0.798 0.804 0.811 0.897 0.901 0.872
0.903 0.869 0.898 0.920

Table : Results on multiple relationships in terms of Macro-F1 and Micro-F1
(with best in bold and second-best underlined)

0.75!  Dataset Metric Train% Deepwalk node2vec GCN GAT HAN HGT R-GCN CompGCN
DeepR PRIM 8*BJ 4*Macro-F1 40% 0.475 0.477 0.579 0.577 0.605 0.620 0.591 0.603
0.626 0.664 50% 0.518 0.522 0.611 0.603 0.641 0.654 0.632 0.647 0.656 0.678
60% 0.543 0.544 0.625 0.621 0.660 0.675 0.657 0.676 0.673 0.694 70% 0.576
0.579 0.632 0.626 0.688 0.701 0.676 0.687 0.697 0.721 4*Micro-F1 40% 0.638
0.637 0.673 0.675 0.678 0.716 0.680 0.699 0.714 0.759 50% 0.667 0.669 0.696
0.691 0.717 0.744 0.709 0.716 0.742 0.789 60% 0.689 0.693 0.718 0.714 0.724
0.761 0.721 0.751 0.752 0.799 70% 0.709 0.711 0.723 0.722 0.756 0.776 0.755
0.764 0.778 0.804 8*SH 4*Macro-F1 40% 0.472 0.471 0.474 0.473 0.524 0.547
0.532 0.531 0.543 0.582 50% 0.503 0.505 0.504 0.506 0.556 0.579 0.555 0.559
0.576 0.604 60% 0.520 0.523 0.526 0.523 0.601 0.622 0.597 0.575 0.623 0.642
70% 0.533 0.537 0.542 0.529 0.613 0.631 0.614 0.611 0.634 0.659 4*Micro-F1 40%
0.635 0.638 0.627 0.622 0.676 0.718 0.688 0.678 0.713 0.753 50% 0.657 0.659
0.674 0.664 0.714 0.742 0.694 0.694 0.744 0.778 60% 0.677 0.681 0.679 0.669
0.744 0.761 0.742 0.746 0.756 0.790 70% 0.702 0.704 0.694 0.695 0.752 0.778
0.768 0.771 0.774 0.806

##  Performance Comparison

The results of different methods are reported in Table , where we use
different fractions of datasets from 40% to 70% as training data (denoted by
“Train%”). We have several observations: (1) Rule-based methods, random walk
based methods and vanilla GNN methods perform poorly on both datasets,
demonstrating that simple rules and methods which do not consider relational
features are not sufficient to provide accurate predictions. (2) Methods that
can model the relation heterogeneity achieve relatively good results as
compared to other baselines. Moreover, the methods that decompose a graph into
sub-graphs for different relation types (DecGCN, DeepR) perform worse than
those that model multiple relation types within a graph in a unified way (HGT,
CompGCN). (3) HGT and CompGCN outperform other baselines in most cases, which
affirms the effectiveness of attention mechanism in aggregation from neighbor
nodes and learning relation representations. (4) Our proposed PRIM outperforms
all the baselines in all metrics and in all cases. Compared to HGT and
CompGCN, PRIM captures unique spatial features of POIs with novel spatial-
aware attention mechanism and spatial context extractor, and employs distance-
specific scoring function to better utilize relation representations.

To further validate the effectiveness of handling multiple relation types,
based on the number of times appearing in user logs, we further create finer-
grained relationships to differentiate the degree of
competitiveness/complementarity between POI pairs. Finally, we end up with 6
relationships and the experimental results for GNN based models are shown in
Table . It can be observed that the proposed PRIM still performs the best,
which further demonstrates the superiority of our model.

##  Model Scalability

Given that the above datasets are relatively small, we study the model
scalability on a large POI dataset in Singapore that consists of 251,219 POIs.
Due to the lack of ground-truth relationships, for each POI we randomly assign
8 relationships to others. We evaluate the scalability for POIs ranging from
50K to 250K.

First, we evaluate the training efficiency for GNN based methods. The results
are shown in Figure . It can be observed that the models that handle
homogeneous graphs achieve the best efficiency. Apart from R-GCN, all the
models that can handle multiple relationships have comparable training
efficiency. Furthermore, as discussed in our time complexity analysis earlier,
our proposed method grows linearly against the input size in terms of the
training time. The results also show that our model achieves a great trade-off
between effectiveness and efficiency in training phrase.

Next, we evaluate the prediction efficiency. Given POI pairs as queries, we
first obtain the embeddings for these POIs according to Section 4.2-4.4. After
that, we only need to index the generated embeddings and apply scoring
functions in Equation (11) & (12) to get predictions. In this case, the
prediction time does not change with the increase of POI size for all the GNN
based methods. Based on such a procedure, we first generate the POI embeddings
and then conduct 10K queries to test the prediction efficiency. The
experiments show that the average prediction time per query is 1.57ms for our
proposed method. In addition, when distance-specific hyperplane projection is
not performed (Equation (11)), which is the case for the other GNN based
methods, the average query time reduces to 0.61ms. In practice, both cases can
satisfy the requirement in daily company services.

In summary, the training and prediction efficiency results demonstrate that
our proposed method can scale to large-scale datasets and can be applied in
practical use.

[width=0.7]figures/scalability_1.pdf

Figure : Training scalability

[width=]figures/ablation_BJ1.pdf

Macro-F1 on BJ

Micro-F1 on BJ

[width=]figures/ablation_SH1.pdf

Macro-F1 on SH

Micro-F1 on SH

Figure : Results of ablation study

##  Ablation Study

We conduct an ablation study by removing different components of PRIM to
demonstrate their contributions to the performance. Specifically, we compare
with the following model variants: 1)-T: the category taxonomy constraint is
removed and we learn each category embeddings independently (Section
LABEL:subsec:taxonomy); 2)-S: the spatial context is removed (Section
LABEL:subsec:spatial_context); 3)-D: the distance-specific hyperplane
projection (Equation LABEL:hyper) is removed (Section
LABEL:subsec:score_func). We also combine these three variants to test the
model with two or more components being removed, i.e., -DS, -DT, -ST and -DST.
In addition, we present the results of the best baseline model for comparison,
denoted by Base.

Figure  presents the results of all variants of our model, and we have the
following observations: (1) The removal of different components leads to
degradation of model performance, which demonstrates the effectiveness of each
component. Moreover, with more components removed, we can see a larger
decrease in model performance. This shows that different components capture
different aspects of characteristics in POIs which are complementary to
enhance the performance. (2) The performance gap between model variants
becomes larger when dealing with less training data. It reveals that context
information extracted from different components would be important to
alleviate the issue when the data is relatively small. (3) The model with all
the three components removed (equivalent to WRGNN) still outperforms the best
baseline in most cases, which validates the superiority of WRGNN over baseline
models.

##  Model Analysis

###  Analysis on sparse cases

To achieve practical usage in applications, PRIM is expected to be robust to
all cases. That is to say, the model should deal with POIs effectively even if
they have very few relationships (i.e., sparsity). For this purpose, we
construct the test set for POIs with fewer than 3 relationships in the
training data to check their performance. We only report the results for 4
best performing baselines, since others show consistently worse results.

[width=]figures/sparse_BJ.pdf

Macro-F1 on BJ

Micro-F1 on BJ

[width=]figures/sparse_SH.pdf

Macro-F1 on SH

Micro-F1 on SH

Figure : Results on sparse cases

The results of sparse POI cases are reported in Figure . We can observe from
the results that our model outperforms all the baselines on both metrics. It
demonstrates that our model is more capable of handling sparse cases.
Furthermore, the performance of the baselines degrades more significantly than
our model. For example, in the Shanghai dataset, the Macro-F1 decreases by
5.1% on average for our model, and by 8.0%, 8.4%, 6.1% and 7.3% respectively
for HAN, HGT, CompGCN and DeepR. These results show the superiority of PRIM
for exploiting the context information specific to POIs to compensate for the
lack of relationships.

###  Analysis on unseen cases

Apart from model robustness on sparse cases, it is also important that the
model should be effective in inductive settings. In other words, the model
should generalize well on POIs which are not seen during training. Such a
capability is desired in real scenarios since the model can be directly
applied to newly arrived POIs without being frequently updated. To compare the
model performance in inductive settings, we randomly hide 20% of the POIs and
remove the relation edges connecting to them in each dataset to serve as
unseen cases. Then we train all the models with the remaining POIs and the
same parameter settings. Same as sparse cases, we only report the results for
4 best performing baselines.

Table : Results on unseen cases

0.75!  Dataset BJ SH Macro-F1 Micro-F1 Macro-F1 Micro-F1 HAN 0.844 0.875 0.794
0.860 HGT 0.837 0.864 0.793 0.852 CompGCN 0.841 0.872 0.790 0.848 DeepR 0.815
0.852 0.764 0.836 PRIM 0.880 0.905 0.814 0.885

Table : Model performance on different areas

0.9!  2*Metric 2*Train% Area BJ core area BJ suburb BJ overall SH 4*Macro-F1
40% 0.846 0.844 0.845 0.722/0.822 50% 0.871 0.869 0.870 0.718/0.844 60% 0.881
0.883 0.882 0.729/0.861 70% 0.896 0.894 0.895 0.741/0.875 4*Micro-F1 40% 0.874
0.890 0.879 0.797/0.886 50% 0.890 0.906 0.895 0.813/0.896 60% 0.903 0.915
0.907 0.801/0.909 70% 0.911 0.918 0.913 0.823/0.920

The results of unseen POI cases are listed in Table . We observe that all the
compared models achieve relatively good results. This is because they are all
based on GNN models, which have been shown to handle inductive settings well
[GraphSAGE]. Specifically, DeepR performs the worst among all the baselines,
while HAN, HGT and CompGCN have comparable performance. In addition, PRIM
outperforms all the compared methods, which verifies the effectiveness of our
model on unseen POIs.

###  Analysis on different regions

To gain a better understanding of our model’s performance on different
regions, we split Beijing into core area and suburb. The core area is less
than 15% of the area in Beijing with over 53% POIs, which is much denser than
the suburb. We report the model’s performance for these two regions in Table .
The results show that the performance gap is very small between the core area
and the suburb, which demonstrates the robustness of our model on regions with
different sparsity.

We further investigate the performance on regions with more complex
differences by directly applying the model learned from Beijing on Shanghai.
The results are presented in Table . In the column of Shanghai dataset, the
first and second values indicate the performance of model trained on Beijing
and model trained on Shanghai itself respectively. We can observe larger
performance drop compared to regions within a city. However, the overall
performance is still good where most Micro-F1 metrics are over 0.8.

[width=]figures/different_scale.pdf

Macro-F1

Micro-F1

Figure : Results on datasets with different characteristics

###  Analysis on datasets with different characteristics

We conduct experiments to evaluate the model performance on datasets with
different characteristics. To this end, we propose to generate datasets by
sampling subsets of the Beijing dataset. Specifically, we randomly select
40%/60%/80% POIs in Beijing and only keep the edges among the selected POIs to
serve as datasets with different scale/density/spatial distances. As a result,
sparser datasets have lower density and larger spatial distances. Then we
split the relational edges into 60%/20%/20% as train/validation/test data. The
experimental results are shown in Figure . We can observe that PRIM
outperforms all the baselines, thus validating the superiority of our proposed
model in handling datasets with different characteristics.


††thanks: The authors gratefully acknowledge University of São Paulo and
support from SHELL Brazil (subsidiary company of Royal Dutch Shell) and
FAPESP, through the Research Centre for Gas Innovation (RCGI) hosted by the
University of São Paulo (FAPESP Grant Proc. 2014/50279-4). We would also like
to thank FAPESP and CNPq for their support, by grants number FAPESP
2016/02175-0 and CNPq 303992/2017-4.

# Estimating the radii of air bubbles in water using passive acoustic
monitoring

P. Hubert  Technology and Data Science Department, EAESP - FGV
Acoustics and Environment Lab - EP - USP  L. Padovese  Acoustics and
Environment Lab - EP - USP

###### Abstract

The study of the acoustic emission of underwater gas bubbles is a subject of
both theoretical and applied interest, since it finds an important application
in the development of acoustic monitoring tools for detection and
quantification of underwater gas leakages. An underlying physical model is
essential in the study of such emissions, but is not enough: also some
statistical procedure must be applied in order to deal with all uncertainties
(including those caused by background noise). In this paper we take a
probabilistic (Bayesian) methodology which is well known in the statistical
signal analysis communitiy, and apply it to the problem of estimating the
radii of air bubbles in water. We introduce the _bubblegram_ , a feature
extraction technique graphically similar to the traditional spectrogram but
tailored to respond only to pulse structures that correspond to a given
physical model. We investigate the performance of the bubblegram and our model
in general using laboratory generated data.

##  I Introduction

Underwater gas leakages are a matter of concern for several industries.
Detecting and quantifying such leakages is thus a problem with theoretical
interest but also with important engineering implications.

In this paper we are particularly interested in the use of Passive Acoustic
Monitoring (PAM) technologies. This approach is cheaper compared to other,
more complex options involving active monitoring, but presents the challenge
of analyzing the acquired signal which will be contaminated by background
noise and other acoustic events at a wide frequency range.

In this work we consider the single sensor case and propose a probabilistic
(Bayesian) approach to the problem of modelling the acoustic behavior of
underwater air bubbles and bubble plumes. Our approach is independent of
physical models; we adopt the model of (Strasberg, 1956) as our first choice
in this paper, but generalization for other models is straightforward.

Based on this approach we propose a graphical evaluation tool, the
_bubblegram_ , which bears resemblance to the usual spectrogram graphs widely
used in signal processing. However, instead of representing the signal’s
energy at different frequency ranges as the spectrogram does, the bubblegram
represents the posterior density for the presence of bubbles of any radius $R$
at a given time $t$.

The paper is organized as follows: section II describes the Bayesian
methodology for the analysis of exponentially decaying sinusoids in noisy
signals. Secion III describe the physical model of (Strasberg, 1956) and shows
how to plug this model into the probabilistic framework. Section IV presents
the bubblegram and section V presents the results of our model in laboratory
data. Section VI concludes the paper.

##  II Probabilistic modelling

Consider a real-valued signal along a single time dimension, $y(t)$, sampled
(not necessarily uniformly) at points $\\{t_{i}\\}$. We adopt the notation
$y_{t}$ for the discretely sampled signal.

Define a pulse as an exponentially decaying sinusoidal function

$x_{t}=A\cdot e^{-\lambda t}\cdot cos(\omega t+\phi)$ (1)

where $A$ is the (assumed constant) amplitude, $\lambda$ is the decay factor,
$\omega$ is the frequency and $\phi$ the pulse phase.

This equation assumes that $t=0$ is the starting time of the pulse.

Instead of writing the pulse equation depending non-linearly on the phase
$\phi$, we adopt the parameterization

$x_{t}=e^{-\lambda(t-t_{0})}\left(A\cdot cos(\omega(t-t_{0}))+B\cdot
sin(\omega(t-t_{0}))\right)\cdot\mathbf{1}_{t>t_{0}}$ (2)

which is equivalent to the former but replaces the nonlinear parameter $\phi$
by a new, linear parameter (the amplitude $B$). Also we included the parameter
$t_{0}$ to generalize the model for pulses starting at times different from
$t=0$, and the function $\mathbf{1}_{U}$ is the indicator function, taking the
value $1$ when $U$ is true, and $0$ otherwise.

The acquired signal $y_{t}$ is assumed to be formed by a single pulse starting
at time $t_{0}$, corrupted with Gaussian white noise $r_{t}$, with
$E\left(r_{t}\right)=0$ and $var(r_{t})=\sigma^{2}$.

$y_{t}=x_{t}+r_{t}$ (3)

In this general form the model can be used to describe many different
phenomena. Its main interest lies in the fact that it explicitly models the
pulse in terms of decay constant and fundamental frequency, as opposed to
traditional spectral analysis that focus on the frequency alone. This model,
being probabilistic in nature, naturally incorporates uncertainty into the
analysis, and also allows the inclusion of any prior information available
about the signal and the noise.

In this work we are particularly interested in the situation where the decay
constant $\lambda$ and the fundamental frequency $\omega$ are related by a
known, deterministic function. This is the case when the pulse is taken to
represent the acoustic emission of a gas bubble in water, in which case both
$\lambda$ and $\omega$ depend on the bubble’s mean (equilibrium) radius.
Directly estimating this radius can aid the design of signal detectors for
underwater gas leakages, and also serves as a tool for estimation of the
leaked gas flow.

Now taking the structural model of equation 3, and assuming Gaussian white
noise, the log-likelihood of the model can be written as

$\ell(y_{t}\mid
A,B,\omega,\sigma,\lambda,t_{0})=-\frac{N}{2}log\left(2\pi\sigma^{2}\right)-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left(y_{t}-x_{t}\right)^{2}$
(4)

Since our main interest lies in the estimation of $\lambda$, $\omega$ and
$t_{0}$, we proceed to marginalize all other factors from the likelihood. By
adopting uniform priors for the amplitudes $A$ and $B$, and a Jeffreys’ prior
for $\sigma^{2}$, this marginalization can be performed analitically,
resulting in the following expression for the marginalized likelihood:

$\displaystyle\begin{split}P(y_{t}\mid
t_{0},\omega,\lambda)\propto&\left(\sum_{i=1}^{N}d_{i}\right)^{-1/2}\left(\sum_{i=1}^{N}e_{i}\right)^{-1/2}\times\\\
&\left(\sum_{i=1}^{N}y_{i}^{2}-\frac{\left(\sum_{i=t_{0}}^{N}y_{i}d_{i}\right)^{2}}{\sum_{i=t_{0}}^{N}d_{i}^{2}}-\frac{\left(\sum_{i=t_{0}}^{N}y_{i}e_{i}\right)^{2}}{\sum_{i=t_{0}}^{N}e_{i}^{2}}\right)\end{split}$
(5)

with

$\displaystyle\begin{split}&d_{i}=e^{-\lambda(t-t_{0})}cos(\omega(t_{i}-t_{0}))\\\
&e_{i}=e^{-\lambda(t-t_{0})}sin(\omega(t_{i}-t_{0}))\end{split}$ (6)

This expression can further on be written in a concise manner as

$P(y_{t}\mid t_{0},\omega,\lambda)\propto K(\omega,\lambda,t_{0})\left[\lVert
y\rVert^{2}-\left(\frac{\langle y,d\rangle}{\lVert
d\rVert}\right)^{2}-\left(\frac{\langle y,e\rangle}{\lVert
e\rVert}\right)^{2}\right]$ (7)

Notice that this expression is independent of the form of both $d_{i}$ and
$e_{i}$, given that these functions do not depend on the marginalized
parameters.

This methodology is not new, as it has been originally proposed and studied by
(Bretthorst, 1988, 1990a, 1990b, 1990c) and applied to the analysis of Nuclear
Magnetic Ressonance data. Here we apply it to the analysis of the acoustic
emission of air bubbles in water.

##  III Acoustic emission of air bubbles

The study of the acoustic emission of air bubbles in liquid media dates back
to (Minnaert, 1933), is treated again by (Strasberg, 1956), and more recently
has been extensively studied by (Leighton _et al._ , 1996; Ainslie and
Leighton, 2011; Leighton _et al._ , 2012) and others.

In (Strasberg, 1956) a formula relating the bubble’s fundamental frequency of
oscilation to its radius is given as:

$\omega_{0}=\sqrt{\frac{3\gamma P_{0}}{\rho_{0}}}\left(\frac{1}{2\pi
R_{0}}\right)$ (8)

here $R_{0}$ is the bubble steady-state radius; $\rho_{0}$ is the liquid’s
density; $\gamma$ is the ratio of specific heats, and $P_{0}$ the static
pressure.

The model for the pressure variations on the surface of a single bubble over
time is of the same form as 2, given by:

$x_{t}=\left(\omega_{0}R_{0}\right)^{2}\frac{\rho_{0}}{r}\alpha_{0}e^{-\omega_{0}\lambda(t_{0}-t)}\mathbf{1}_{t>t_{0}}cos\left(\omega_{0}(t_{0}-t)\right)$
(9)

with $\alpha_{0}$ the initial wall amplitude of the bubble, and $r$ the
distance between the bubble and the sensor.

By assuming constant the liquid’s density, the distance between the bubble and
the sensor, and the initial amplitude, we rewrite equation 9 as

$x_{t}=C\cdot\left(\omega_{0}R_{0}\right)^{2}e^{-\omega_{0}\lambda(t-t_{0})}\mathbf{1}_{t>t_{0}}cos\left(\omega_{0}(t-t_{0})\right)$
(10)

Now including a phase parameter and reparameterizing, we see that model 10 can
be put in the form of 5 with

$\displaystyle\begin{split}&d_{i}=\left(w_{0}R_{0}\right)^{2}e^{-\omega_{0}\lambda(t-t_{0})}cos\left(\omega_{0}(t-t_{0})\right)\\\
&e_{i}=\left(w_{0}R_{0}\right)^{2}e^{-\omega_{0}\lambda(t-t_{0})}sin\left(\omega_{0}(t-t_{0})\right)\end{split}$
(11)

Using equation 8 that gives $\omega_{0}$ as a function of $R_{0}$,
$\omega_{0}=\Omega(R_{0})$, we can eliminate $\omega$ from the model:

$P(y_{t}\mid t_{0},R_{0},\lambda)\propto K(t_{0},R_{0},\lambda)\left[\lVert
y\rVert^{2}-\left(\frac{\langle y,d\rangle}{\lVert
d\rVert}\right)^{2}-\left(\frac{\langle y,e\rangle}{\lVert
e\rVert}\right)^{2}\right]$ (12)

with

$\displaystyle\begin{split}&d_{i}=\left(\Omega(R_{0})R_{0}\right)^{2}e^{-\Omega(R_{0})\lambda(t-t_{0})}cos\left(\Omega(R_{0})(t-t_{0})\right)\\\
&e_{i}=\left(\Omega(R_{0})R_{0}\right)^{2}e^{-\Omega(R_{0})\lambda(t-t_{0})}sin\left(\Omega(R_{0})(t-t_{0})\right)\end{split}$
(13)

Finally, (Strasberg, 1956) also provides an expression for the decay constant
as a function of the fundamental frequency. In our notation

$\lambda=0.014+1.1\cdot 10^{-5}\Omega(R_{0})$ (14)

By plugging this equation into our model we are able to eliminate one more
parameter, writing $\lambda=\Lambda(R_{0})$, and arriving at a final model
depending only on the bubble’s radius $R_{0}$ and its time of occurrence,
$t_{0}$.

This elimination depends crucially on the physical models relating the
frequency $\omega$ and decay factor $\lambda$ to the bubble’s radius. In this
formulation we accepted these models as groundtruth, and did not include any
uncertainty in these functions. It is possible, however, to augmentate the
model by including error terms in the equations for $\omega$ and $\lambda$,
and also to include prior distributions for all physical constants involved in
the model. For this work, however, we chose to adopt the simplest formulation
to allow a direct comparison of our model with the traditional tools of signal
analysis.

##  IV Estimating the bubble radius

An immediate use for this model is the estimation of $R_{0}$ and $t_{0}$ from
a given signal. This can be performed by first adopting prior distributions
for each parameter in the model, and obtaining a joint posterior by
application of Bayes’ theorem.

To illustrate this process we generate random data for a single pulse
representing a bubble starting at $t_{0}=0.5$ with radius $R_{0}=1mm$. We add
Gaussian white noise with $\sigma^{2}=0.1$.

By assuming flat, uninformative priors for all parameters, we obtain that the
(log) joint posterior is equal (up to a constant) to the (marginal) log-
likelihood from equation 12. The values of this log-posterior over a grid of
$500\times 500$ points in the $time$ $\times$ $radius$ plane appear as the
color scale in figure 1, where brighter colors are associated with high
posterior mass for a bubble with a given radius at a given time. The true
values of $t_{0}$ and $R_{0}$ are marked with a white cross.

In figure 1 we also included a spectrogram of the same data for comparison.

figure1a.png.4(a) Log-posterior figure1b.png.4(b) Spectrogram

Figure 1:  Log-posterior and spectrogram of simulated, single pulse data.

The first and most important difference between the visualization of the log-
posterior and the spectrogram is the noise. Since the log-posterior model
incorporates the noise term explicitly in the calculations (via the
marginalization of the noise power $\sigma^{2}$) it is able to better recover
the signal parameters. The spectrogram in itself is not a statistical
procedure, and thus does not account for the presence of noise. The
consequence is that signal and noise get blurred, as we see in the figure.

It will be also interesting to analyze a signal containing a sequence of
bubbles of various sizes. This will occurr in most practical situations with
higher gas flows. The original model we proposed in section II assumes a
single bubble; however, the presence of the parameter $t_{0}$ allows the
application of this same model for a longer signal. One can think of the usual
Short-Term Fourier Transform approach, where a local model is applied to small
sections of a long signal, possibly after applying a smoothing function in the
time domain. The smoothing, in our model, would be represented by a prior
distribution on $t_{0}$.

What we show in the next figures is the plot of the posterior log-density for
$R_{0}$, conditional on a fixed value for $t_{0}$. In this sense we are
looking for the radii with larger evidence of appearing in a pulse starting at
time $t$, for $t\in\left[t_{a},t_{b}\right]$.

When the gas flow is large enough and the bubbles start to cluster together,
the hydrodinamical models for the acoustic emission of a single bubble are no
longer exactly valid (Ainslie and Leighton, 2011). However, since the acoustic
emission of a gas bubble in water depends only on the forming process of the
bubble, spurious acoustic emissions will occur only in the events of
fragmentation of a bubble or the coalescence of two or more bubbles. For these
reasons we believe that our model can be applied succesfully even to the
analysis of higher flow bubble plumes.

To understand the behavior of the log-posterior in this case we simulate a
sequence of $10$ pulses occurring along a one second interval. The bubbles’
starting times and radii are simulated from uniform distributions between
$[0,1]$ and $[0.5,1.5]$, respectively. The evaluation of the posterior over a
grid of $500\times 500$ points and $R_{0}\in[0.2,2]$ and $t_{0}\in[0,1]$
appear in figure 2, where we included empty circles centered at the true
values of each bubble. Again the spectrogram estimated from the same signal is
included for comparison.

figure2a.png.4(a) Log-posterior figure2b.png.4(b) Spectrogram

Figure 2:  Log-posterior and spectrogram for simulated data.

Both figures 1 and 2 show that the model behaves as expected. In the multiple
bubbles case it appears that larger bubbles might be easier to detect (since
the log-posterior is higher around these bubbles), but local peaks in the log-
likelihood can be found near the radius and time of occurrence of all pulses.
By comparing the log-posterior plots with the corresponding spectrograms we
see that the log-posterior not only directly estimates the bubble’s radii but
also is more precise, since the model explicitly incorporates the background
noise as already observed.

These plots that evaluate the joint posterior of $R_{0}$ and $t_{0}$ over a
(not necessarily uniform) grid show some visual resemblance with the
traditional spectrogram plots usually applied to the analysis of acoustic
signals. Inspired by this similarity we call the plots of the log-posterior
_bubblegrams_ , since they are designed to respond to the acoustic emission of
air bubbles in water, representing the posterior probability associated with a
given radius at a given time. One advantage of the bubblegram over the
traditional spectrogram is that it does not suffer from the usual resolution
limitations of the spectrogram; to understand how that can be the case, one
can imagine that by the use of the probabilistic model one is effectively
introducing an _infinite zero-padding_ into the signal, thus allowing any
arbitrary level of resolution both in time and radius (frequency) domains.

##  V Experimental results

In this section we apply our model to the analysis of experimental data. The
experiments were conducted by the _Acoustics and Environment Laboratory_
(LACMAM in the portuguese acronym) at the Polytechnic School, University of
São Paulo.

We used sensors (hydrophones) developed by the laboratory itself; these
sensors have a frequency band of $5$ $Hz$ to $60$ $kHz$, with sensitivity of
$-157\pm 2$ $dB$ $rel$ $1\mu Pa$ (preamplified) or $-203\pm 2$ $dB$ $rel1\mu
Pa$ (not preamplified). A picture of one sensor appears in figure 3.

figure3.png.5

Figure 3:  Hidrophone developed by LACMAM.

Two experiments were conducted. The first was designed to measure the size of
single bubbles by using a high speed CCD camera with a $25$ $mm$ lens with
$F1.4$ opening. For this experiment a small aquarium with dimensions $45$ $cm$
(length) $\times$ $40$ $cm$ (width) $\times$ $25$ $cm$ (height) was used.
Three nozzles with diameters $2.5$, $4$ and $6$ $mm$ were positioned one at a
time under a $22$ $cm$ water column, and a flow of less than $1$ $l/min$ of
air was induced by the use of a flow controller connected to a compressed air
cylinder. A picture of the setup is seen on figure 4.

figure4.png.5

Figure 4:  Setup for measuring bubble sizes.

The second experiment was designed to obtain acoustic signals emitted by
bubble plumes in water. It was performed in a diving tank with a depth of $5$
$m$. The bubbles are created by the injection of compressed air into the tank,
with flow, pressure and exit diameter orifice controlled. Figures 5 and 6 show
the equipment used in this experiment. In figure 5 the black lines show the
position of the four sensors used, and the red circle indicates the
approximate position of the flow controller nozzle.

figure5.png.5

Figure 5:  Diving tank; see text for details.

figure6.jpeg.5

Figure 6:  Schematic drawing of the experimental setup.

###  V.1 Data acquisition and experimental results

####  V.1.1 Single bubble measurements

We recorded high speed ($1000$ frames per second) videos of the bubbles
generated in the aquarium, with a ruler attached to the aquarium wall. Figure
7 below shows snapshots of bubbles generated with the three different nozzle
diameters. The labels in each picture indicate the diameter of the nozzle used
in the air injection.

figure7.png.5

Figure 7: Images of single bubbles

The acoustic emissions were recorded with a sampling frequency of $48$ $kHz$.
To create the bubblegrams we first passed the signals through a Butterworth
bandpass filter (with pass range $200$ to $3000$ $Hz$). We then selected the
sections of the signal corresponding to the photographed bubbles. These
preprocessing steps were taken in order to accelerate further processing of
the signal.

The bubblegrams obtained from the preprocessed signals appear in figure 8.

figure8a.png.4(a)$2.5$ $mm$ nozzle figure8b.png.4(b) $4$ $mm$ nozzle
figure8c.png.4(c) $6$ $mm$ nozzle

Figure 8: Bubblegrams for a single bubble with multiple nozzle diameters

The Maximum Posterior (MAP) estimates for the bubbles radii were $0.6509$,
$0.4273$ and $0.5102$ $mm$ for the $2.5$, $4$ and $6$ $mm$ nozzles
respectively. By comparing this estimates to the pictures in figures 8 we find
that the estimates are accurate, considering the deformation of the bubbles in
each snapshot.

####  V.1.2 Bubble plumes

To obtain the acoustic emission of bubble plumes at higher flow rates the
diving tank was used. The nozzle of the pneumatic system was positioned at the
bottom of the tank, and several runs were conducted with different nozzle
diameters and gas flow rates. For each run, we open the flow controller and
let the gas escape for two minutes. Then we close the controller and wait for
one minute to open it again.

In this section we are interested in studying both the spectrogram and
bubblegrams for different values of the nozzle diameter and gas flow rates. We
analyze signals obtained with diameter $D\in\\{2.5,6,11\\}$ $mm$ and flows of
$F\in\\{2,5,10\\}$ $l\min$, ie, we analyze $9$ different sections of the
experimental signal.

To conduct our analysis we extract small sections ($2$ seconds) of the signals
and start by obtaining the spectrograms of each section. These spectrograms
appear in figure 9.

sgram_v2_d2_5.png4cm(a) $d=2.5$ $mm$; $f=2$ $l/min$ sgram_v2_d6.png4cm(b)
$d=6$ $mm$; $f=2$ $l/min$ sgram_v2_d11.png4cm(c) $d=11$ $mm$; $f=2$ $l/min$
sgram_v5_d2_5.png4cm$d=2.5$ $mm$; $f=5$ $l/min$ sgram_v5_d6.png4cm$d=6$ $mm$;
$f=5$ $l/min$ sgram_v5_d11.png4cm$d=11$ $mm$; $f=5$ $l/min$
sgram_v10_d2_5.png4cm$d=2.5$ $mm$; $f=10$ $l/min$ sgram_v10_d6.png4cm$d=6$
$mm$; $f=10$ $l/min$ sgram_v10_d11.png4cm$d=11$ $mm$; $f=10$ $l/min$

Figure 9: Full spectrograms for bubble plume data.

The spectrograms show that the signal’s energy is concentrated around smaller
frequencies (between $250$ and $1000$ $Hz$). However, for the reasons already
exposed, the graphs are noisy and do not allow easy or direct estimation of
the probable sizes of the bubbles.

The bubblegrams of figure 10, on the other hand, are much clearer and easier
to inspect. As we can see, higher flow rates are associated with bigger
bubbles.

bgram_v2_d2_5.png4cm(a) $d=2.5$ $mm$; $f=2$ $l/min$ bgram_v2_d6.png4cm(b)
$d=6$ $mm$; $f=2$ $l/min$ bgram_v2_d11.png4cm(c) $d=11$ $mm$; $f=2$ $l/min$
bgram_v5_d2_5.png4cm$d=2.5$ $mm$; $f=5$ $l/min$ bgram_v5_d6.png4cm$d=6$ $mm$;
$f=5$ $l/min$ bgram_v5_d11.png4cm$d=11$ $mm$; $f=5$ $l/min$
bgram_v10_d2_5.png4cm$d=2.5$ $mm$; $f=10$ $l/min$ bgram_v10_d6.png4cm$d=6$
$mm$; $f=10$ $l/min$ bgram_v10_d11.png4cm$d=11$ $mm$; $f=10$ $l/min$

Figure 10: Bubblegrams for bubble plume data.

The bubblegram is a useful representation of the signal, directly obtaining a
probability distribution over bubble sizes. From the figures we can observe
that higher gas flow rates create bubble plumes with more diverse sized
bubbles (we see that from the regions with relevant evidence, colored blue,
but far from the more concentrated regions that show the typical bubble size
in each case). Also the time between the formation of each bubble changes
according to the gas flow.

The gas flow appears to have a more evident effect on the bubble size
distribution than the nozzle diameter, as we can see by comparing figures in
the same row (i.e., signals obtained with the same flow) and then comparing
figures in the same column (signals obtained with the same nozzle diameter).
However both gas flow and nozzle diameter are positively correlated with
bubble sizes, even though the flow rate has a stronger effect.

By aggregating and normalizing the log-posterior for the bubble sizes we
arrive at the approximate density for the distribution of bubble sizes in the
entire signal. In figure 11 we plot these distributions for a fixed flow rate,
and varying nozzle diameter.

p_d2_5.png.4$d=2.5$ $mm$ p_d6.png.4$d=6$ $mm$ p_d11.png.4$d=11$ $mm$

Figure 11: Log-posterior density for bubble sizes, fixed flow

With these figures we are able to observe more directly the effect of both
nozzle diameter and flow rate in the estimated bubble size distribution.

Our results indicate that this relationship is complex. For instance, from
figure 11 we observe first of all that the smallest nozzle induces more
homogeneous bubbles populations, since the curves are close to being unimodal
(specially for the lowest flow, blue line). For this smallest flow, the mode
radius is also the smallest, around $0.25$ $mm$; the two bigger flows have
very similar curves.

Increasing the nozzle diameter we observe that the bubble populations start to
break down into subpopulations with different sizes (that is, the curves start
to show large deviations from unimodality). The curve for the smallest flow
(blue lines in all three plots) has always a prominent peak around $0.25$
$mm$, but as the nozzle size increases bigger bubbles can also be found
according to our analysis. The intermediate and higher flows that we tested
behave less monotonically; for a flow of $5$ $l/min$, when the diameter is set
at $6$ $mm$ (figure 11, red line) there appears a high peak around a radius of
$1.25$ $mm$, indicating that increasing the nozzle allows the formation of
bigger bubbles. However, when the diameter increases to $11$ $mm$, a
subpopulation of bubbles with radii around $0.75$ $mm$ appears, along with
another subpopulation of bubbles with raddi around $1.75$ $mm$. Therefore, for
the bubble plumes with intermediate flow, the bubbles increase in size with
nozzle diameter, but for the biggest diameter a population of small bubbles is
also formed.

The curves obtained from the analysis of the samples with biggest flow (green
line in all three subplots) all have an important peak around $1$ $mm$.
Increasing the nozzle diameter from $2.5$ to $6$ $mm$ causes the appearance of
subpopulations of smaller bubbles (peaks around $0.25$ and $0.7$ $mm$, green
line in figure 11). When the nozzle increases to $11$ $mm$ (figure 11) the
peak around smaller bubbles is much attenuated, and is replaced by peaks
around $1.25$ and $1.75$ $mm$.

##  VI Conclusion

We applied a known methodology of probabilistic (Bayesian) signal processing
to the analysis of the acoustic emission of air bubbles in water. Our model
directly expresses the acoustic signal in the time domain as a function of the
bubble radii and time of occurrence, thus providing a useful technique to
analyze acoustic signals containing the emission of bubble plumes.

Testing our model on laboratory generated data we found evidence that the
model is accurate, providing good estimates of the size of a single bubble
from its acoustic emission.

To analyze the acoustic emission of bubble plumes we propose to evaluate the
posterior of the single bubble model in small temporal segments of the signal
in the time domain. This is done in a similar way as the usual Short Time
Fourier Transform, and can be interpreted as imposing an informative prior
over the time of occurrence of a single bubble. This evaluation, performed
over a grid of values for both the time of occurrence and bubble radius,
generates a visual representation we called _bubblegram_ , that is also
similar to the usual spectrogram, but is built on a full probabilistic setup
that embodies a physical model relating the bubble size to the parameters of
an exponentially decaying sinusoidal. As such, the bubblegram provides a
representation of the signal that is more accurate than the spectrogram, and
can be helpful in the study of the bubble sizes in bubble plumes by passive
acoustic methods.

## References

  * Ainslie and Leighton (2011) Ainslie, M. A., and Leighton, T. G. (2011). “Review of scattering and extinction cross-sections, damping factors, and resonance frequencies of a spherical gas bubble.,” The Journal of the Acoustical Society of America 130, 3184–208.
  * Bretthorst (1988) Bretthorst, L. (1988). _Bayesian Spectrum Analysis and Parameter Estimation_ , 1a ed. (Springer-Verlag).
  * Bretthorst (1990a) Bretthorst, L. (1990a). “Bayesian analysis. i. parameter estimation using quadrature nmr models,” Journal of Magnetic Resonance 88, 533–551.
  * Bretthorst (1990b) Bretthorst, L. (1990b). “Bayesian analysis. ii. signal detection and model selection,” Journal of Magnetic Resonance 88, 552–570.
  * Bretthorst (1990c) Bretthorst, L. (1990c). “Bayesian analysis. iii. applications to nmr signal detection, model selection and parameter estimation,” Journal of Magnetic Resonance 88, 571–595.
  * Leighton _et al._ (1996) Leighton, T., Phelps, A., Ramble, D., and Sharpe, D. (1996). “Comparison of the abilities of eight acoustic techniques to detect and size a single bubble,” Ultrasonics 34(6), 661 – 667, http://www.sciencedirect.com/science/article/pii/0041624X96000534, doi: https://doi.org/10.1016/0041-624X(96)00053-4.
  * Leighton _et al._ (2012) Leighton, T. G., Baik, K., and Jiang, J. (2012). “The use of acoustic inversion to estimate the bubble size distribution in pipelines,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 468(2145), 2461–2484, https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2012.0053, doi: 10.1098/rspa.2012.0053.
  * Minnaert (1933) Minnaert, M. (1933). “On musical air-bubbles and the sounds of running water,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 16(104), 235–248, https://doi.org/10.1080/14786443309462277, doi: 10.1080/14786443309462277.
  * Strasberg (1956) Strasberg, M. (1956). “Gas bubbles as sources of sound in liquids,” The Journal of the Acoustical Society of America 28(1), 20–26, https://doi.org/10.1121/1.1908212, doi: 10.1121/1.1908212.


# Multi-focus averaging for multiple scattering suppression in optical
coherence tomography

Lida Zhu  1 Shuichi Makita  1 Junya Tamaoki  3 Antonia Lichtenegger  1,2
Yiheng Lim  1 Yiqiang Zhu  1 Makoto Kobayashi  3 and Yoshiaki Yasuno1
1Computational Optics Group, University of Tsukuba, Tsukuba, Ibaraki, Japan
2Center for Medical Physics and Biomedical Engineering, Medical University of
Vienna, Vienna, Austria
3Department of Molecular and Developmental Biology, Institute of Medicine,
University of Tsukuba, Japan
<EMAIL_ADDRESS>
###### Abstract

Multiple scattering is one of the main factors that limits the penetration
depth of optical coherence tomography (OCT) in scattering samples. We propose
a method termed multi-focus averaging (MFA) to suppress the multiple-
scattering signals and improve the image contrast of OCT in deep regions. The
MFA method captures multiple OCT volumes with various focal positions and
averages them in complex form after correcting the varying defocus through
computational refocusing. Because the multiple-scattering takes different
trajectories among the different focal position configurations, this averaging
suppresses the multiple-scattering signal. Meanwhile, the single-scattering
takes a consistent trajectory regardless of the focal position configuration
and is not suppressed. Hence, the MFA method improves the signal ratio between
the single- and multiple-scattering signals and improves the image contrast. A
scattering phantom and a postmortem zebrafish were measured for validation of
the proposed method. The results showed that the contrast of intensity images
of both the phantom and zebrafish were improved using the MFA method, such
that they were better than the contrast provided by the standard complex
averaging method. The MFA method provides a cost-effective solution for
contrast enhancement through multiple-scattering reduction in tissue imaging
using OCT systems.

††journal: opticajournal

https://optics.bk.tsukuba.ac.jp/COG/

##  1 Introduction

Deep tissue imaging has long been of strong interest in biomedical optics[1].
Optical coherence tomography (OCT) provides three-dimensional images of a
sample with high penetration and high resolution, and has become widely used
for the non-invasive imaging of biological samples. OCT has been successfully
adapted to investigations of various tissues, such as human skin[2], ex vivo
brain tissue[3, 4], and in vitro tumor spheroids[5, 6, 7].

Typically, an OCT image is mainly reconstructed from the single-back-
scattering (SS) signals of the sample. The SS signals are retrieved from the
photons that have been back-scattered once at the sample object and hence
directly carry the object information. However, when imaging in scattering
samples such as tissues, the multiple-scattering (MS) signals harm the image
contrast especially in deep regions[8]. The MS-originated signals at a
position in the image are not scattered at the corresponding position in the
sample and thus do not convey the correct object information. MS limits the
resolvable imaging depth of OCT in tissues and hampers the visualization of
deep microstructures[9, 10].

In standard OCT with a confocal configuration, most of the MS signals can be
rejected by the confocal gating. However, the residual MS signals still
disturb the imaging. Techniques for further reducing the MS effect have been
explored. One approach is modulating the probing light to decorrelate the MS
signals using wavefront manipulation devices, such as a spatial light
modulator (SLM) [11, 12] or a deformable mirror [13]. However, this approach
may result in high system complexity and cost. A solution with inexpensive
optics that can be easily implemented in the standard OCT scheme is more
preferable.

In this paper, we propose a method termed “multi-focus averaging” (MFA),to
reduce the MS signals and improve the image contrast in deep regions of
scattering samples. The MFA method uses a low-cost electrical tunable lens to
decorrelate the MS signals over multiple volumetric acquisitions. Here, the
depth position of the focus is modulated among the multiple volumes. A
computational refocusing is then adopted to cancel out the different defocus
in the focus-modulated OCT volumes. A scattering phantom and a postmortem
zebrafish were used for validation. The phantom results show that the MFA
method achieved a 4.5-dB enhancement of the signal-to-background ratio (SBR).
The postmortem zebrafish results show that the MFA method better visualized
anatomic structures in deep regions, suggesting that this method can be used
for biological imaging.

##  2 Principle and core methods

###  2.1 Principle of the MFA method

The MFA method is a combination of sequential volumetric OCT imaging with
different defocus, computational refocusing of each volume, and complex
averaging of the volumes. In an OCT volume of a scattering sample, the complex
en face OCT signal $S(x,y;z,z_{d})$ at an imaging depth of $z$ comprises two
components, namely $S_{SS}(x,y;z,z_{d})$ and $S_{MS}(x,y;z,z_{d})$ as

$S(x,y;z,z_{d})=S_{SS}(x,y;z,z_{d})+S_{MS}(x,y;z,z_{d}),$ (1)

where $x$ and $y$ are the lateral positions of the scanning locations, and
$z_{d}$ is the amount of defocus at the imaging depth $z$. $S_{SS}$ and
$S_{MS}$ are the signal components originating from SS and MS photons,
respectively. According to the formulation of Ralston et al.[14], the $S_{SS}$
component can be expressed as a convolution of a defocus-free OCT SS signal
and a depth- and defocus-dependent quadratic phase function $\phi$ as

$S(x,y;z,z_{d})=S_{SS}(x,y;z,0)*\exp\left[i\phi(x,y;z,z_{d})\right]+S_{MS}(x,y;z,z_{d}),$
(2)

where $S_{SS}(x,y;z,0)$ is the defocus-free SS OCT signal.

Figure 1:  Schematics of trajectories of SS and MS photons with different
focus depths; i.e., different defocus amounts. No matter the focal position,
the SS photon is scattered only once and hence has the same path length. This
results in the consistent phase of the SS signal after computational
refocusing. Meanwhile, the trajectories of MS photons are scrambled by changes
in the focal position. Hence, the phase of the MS signal is randomized even
after the computational refocusing.

The SS photons experience only a single scattering event, and hence, they take
consistent trajectories among the multiple measurements with different
defocus, as depicted by the SS trajectories in Fig. 1. Hence, once the defocus
is corrected by computational refocusing, the phase of SS components become
consistent among the volumes.

Meanwhile, the MS photons can be scattered by different combinations of
scatterers if the defocus changes the MS trajectories as shown in Fig. 1.
Hence, the trajectories of the MS photons are altered by the defocus. This
alteration of the trajectories results in an inconsistent phase of $S_{MS}$
components for the different defocus, even after computational refocusing.
These properties of the signal components are formulated as follows.

Computational refocusing can be described as the deconvolution of the
quadratic phase function $\phi$. After the computational refocusing, the OCT
signal becomes

$S^{\prime}(x,y;z,z_{d})=S_{SS}(x,y;z,0)+S_{MS}(x,y;z,z_{d})*\exp\left[-i\phi(x,y;z,z_{d})\right],$
(3)

where the deconvolution is represented as the convolution with the complex
conjugate of $\phi$ (see the Appendix for details of this representation). No
matter the original defocus amount, the computational refocusing provides the
consistent defocus-free SS component. In addition, the signal strength of the
SS component is enhanced by the refocusing. Meanwhile, the phase of the MS
component depends on the original defocus amount because the MS trajectories
depend on the defocus amount. Furthermore, because the scattering events of
the MS photons occur at shallower depths than the SS components, the MS signal
may have a defocus amount different from that of the SS signal. Hence, the
computational refocusing tailored for the SS component does not enhance the
signal strength of the MS component.

In our MFA method, an MS-signal-suppressed OCT volume is obtained by complex
averaging the multiple volumes with different defocus amounts after the
computational refocusing as

$\overline{S^{\prime}}(x,y;z)=S_{SS}(x,y;z,0)+\frac{1}{N}\sum_{j=0}^{N-1}S_{MS}(x,y;z,z_{d,j})*\exp\left[-i\phi(x,y;z,z_{d,j})\right],$
(4)

where $N$ is the number of volumes being averaged and $z_{d,j}$ is the defocus
amount of the $j$-th volume. As the phase of the MS term (the second term) is
unpredictable and practically random, the amplitude of this term is reduced by
the averaging. Meanwhile, the phase of SS term (the first term) is consistent
among the volumes, and hence, its amplitude is not reduced.

In practice, a larger difference in the defocus amounts among the volumes may
cause higher mutual randomness (i.e., higher decorrelation) of the phases of
the MS components, and a better suppression of MS components can thus be
achieved. This hypothesis is experimentally validated in Section 4.2.

###  2.2 Implementation of the MFA method

####  2.2.1 OCT setup

In this study, we used a polarization sensitive OCT (PS-OCT) system, which has
been described in Refs. [15, 16]. However, we only used a single-polarization-
channel OCT image (i.e., a conventional non-polarization sensitive OCT image)
to demonstrate the MFA method. The light source was a sweeping laser source
(AXP50124-8, Axsun Technologies, MA) with a central wavelength of
1.31-$\muup$m and a scanning rate of 50 kHz. The objective (LSM03, Thorlabs,
NJ) had an effective focal length of 36 mm and gave a depth of focus (DOF) of
0.36 mm in air. The lateral and axial resolutions were respectively 18
$\muup$m and 14 $\muup$m in tissue.

An electrical tunable lens (ETL, EL-10-30-CI-NIR-LD-MV, Optotune Switzerland
AG, Switzerland) was set on the sample arm to modulate the focal position. It
is noteworthy that the ETL can be easily aligned and applied to standard OCT
systems.

####  2.2.2 Core methods of the MFA

In our MFA method, multiple complex OCT volumes are acquired with different
focal positions. For each volume, the bulk phase error is estimated and
corrected using a smart-integration-path [17] method.

The defocus is corrected in each volume by applying computational refocusing,
where a series of phase-only deconvolution filters are applied at all depths
to correct the defocus. Details of the refocusing method can be found in
Section 2.2 of Ref. [18].

When modulating the focal position using an ETL, the deformation of the lens
affects the optical path length of the probe beam. This introduces an axial
image shift, which sometimes reaches a few tens of microns among the focus-
modulated volumes. In our MFA method, the axial shift is estimated and
corrected using a sub-pixel intensity cross-correlation method[19]. In this
process, B-scans from different volumes at the same location are extracted,
and four-fold up-sampling along the depth direction is conducted. Intensity-
based cross-correlations among the extracted B-scans are performed to estimate
shifts with resolution of 1/4 of a pixel along the axial direction. As this
image shift is consistent within each volume, the estimated shifts are used
for co-registration of all B-scans in each volume.

The inter-volume phase offsets are also corrected. The offsets are computed
for each A-line using an intensity-weighted phase difference averaging. In
this averaging, one volume is used as a reference, and the product of the
reference and the complex conjugate of another volume is computed. The product
is then averaged in complex form along the depth direction within a certain
depth range, which was 30 pixels with sufficient intensity in the present
implementation. The phase of the averaged signal is the phase offset. The
phase offsets are then computed for each A-line and corrected.

Finally, all volumes are averaged in complex form to obtain an “MFA volume.”
The OCT intensity volume is then generated from the MFA volume.

###  2.3 Validation study design

####  2.3.1 Samples

Figure 2:  Schematic (a) and photograph (b) of a scattering phantom. The
phantom comprises glass slips A and D, a scattering layer B, which is a
mixture of polystyrene micro-particles and ultrasound gel, and a glass plate C
embedded in the scattering layer. The glass plate C provides a scattering-free
area. A postmortem zebrafish at 30 days post fertilization (dpf) was used as a
biological sample. The sample is shown in a color photograph (c) and a wide-
field-of-view OCT intensity projection (d). The orange box around the belly
region denotes the measurement area in the validation study of the MFA method.

Two types of sample were measured to validate the performance of the MFA
method. One was the scattering phantom illustrated in Fig. 2(a), which
comprises glass slips (A, D), a scattering layer (B), and a glass plate buried
at the bottom of the scattering layer (C). The glass slip at the top (A) was
tilted to prevent a strong reflection from the surface. The scattering layer
(B) was a mixture of 0.025 mL of 10%-concentration polystyrene micro-particles
(diameter of 10 $\muup$m, 72968-10ML-F, Sigma-Aldrich) and 0.6 mL of
ultrasound gel (Pro Jelly, Jex, Japan). The glass plate (C) had a thickness of
approximately 0.12 mm and was buried at the bottom of the scattering layer to
create a space without scattering. The thickness of the scattering layer above
the glass plate was approximately 1.5 mm. A photograph of the sample from the
top is shown in Fig. 2(b).

A 30-dpf wild-type postmortem female zebrafish [Fig. 2(c)] was imaged to
demonstrate the measurement of a biological sample. The zebrafish was
anesthetized by tricaine and sacrificed by low-temperature treatment (placed
in ice for 2 minutes) and then placed in a petri dish and immersed in saline
solution for measurement. A piece of black tape was placed on the bottom of
the petri dish to prevent strong specular reflection. A wide-field-of-view en
face OCT-intensity projection of the zebrafish is shown in Fig. 2(d), where
the orange box indicates the measured area.

The postmortem zebrafish measurement was performed following the animal
experiment guidelines of the University of Tsukuba. The measurement protocol
was approved by the Institutional Animal Care and Use Committee (IACUC) of
University of Tsukuba.

####  2.3.2 Measurement protocol

In both the phantom and zebrafish measurements, the samples were placed on a
linear translation stage and measured. Seven volumes with different focal
positions were acquired for averaging. The interval between focal positions is
referred to as the focus shifting step ${\Delta z}$ and was set at 0.12 mm
(approximately $1/3$ of the DOF). The overall shifting distance of the focus
is denoted as $D=(N-1){\Delta z}$ and was 0.72 mm in the study. These
parameters were determined empirically, and details of the parameter selection
strategy are given in Section 4.2. For comparison, another set of raw volumes
were acquired and averaged without focus modulation. This averaged volume
without focus modulation is referred as the “complex averaging” volume .

In the measurements, the focal positions were set in the deep regions of the
samples to enhance the light collection efficiency. The lateral scanning range
was 1.5 mm $\times$ 1.5 mm sampled with 512 $\times$ 512 A-lines, which gave
an isotropic lateral pixel separation of 2.93 $\muup$m (around 1/6 of the
lateral optical resolution). The acquisition time for each volume was
approximately 6.5 s.

####  2.3.3 Quantitative image contrast analysis

We defined the SBR to quantify the image contrast in the scattering phantom
measurement. An en face image at a depth slightly beneath the top surface of
the buried glass plate was extracted. Several scatterers with similar
intensities were manually selected to compute the SBR. The signal intensity of
each scatterer was defined as the averaged image intensity within a 3-pixel
$\times$ 3-pixel window centered at the scatterer. The signal intensity used
to compute the SBR was then defined as the mean of the signal intensities of
the selected scatterers. The background intensity was defined as the mean
intensity of pixels in the buried glass plate region (i.e., the scattering-
free region).

##  3 Results

###  3.1 Scattering phantom

Figure 3:  (a)–(c) and (d)–(f) show the B-scans and en face images of the
phantom, respectively. In (a)–(c), the orange arrowheads denote the
approximate focal positions, and the yellow boxes denote a deep region of the
scattering layer. In (d)–(f), the yellow boxes and the magnification insets
are the en face images of a small region in the scattering layer containing
several scatterers. In (d), the red circles denote the scatterers selected to
compute the signal intensity. The blue line indicates a region with 2 $\times$
200 A-lines in the glass plate region used to compute the background intensity
for the SBR analysis. (g) shows the intensity depth profiles that are averaged
by the A-lines at the en face locations denoted by the blue line in (d). The
0-mm depth refers to the top surface of the cover glass on the top, and the
green arrowhead denotes the depth in air where the intensity is 0 dB. The blue
dashed lines in (a) and (g) denote the depth where the en face images are
taken.

Figure 3(a)–(c) shows the OCT intensity B-scans of the single acquisition,
complex averaging, and MFA methods, respectively. The approximate focal
positions are denoted by the orange arrowheads. Multiple orange arrowheads in
Fig. 3(c) denote different focal positions of each image used to generate the
MFA image. Among the B-scans, the MFA image provides a lower background noise
in the scattering layer than the other images, especially in the deep regions
(denoted by the yellow dashed boxes).

Figure 3(g) presents the intensity depth profiles obtained using the three
methods. The depth intensity profile of the MFA (orange) shows a lower
intensity than the other profiles as depth increasing (the region with yellow
background). Since the MS signal is expected to be more pronounced in deeper
regions, this profile may indicate the MS suppression by MFA method.
Meanwhile, the single acquisition (black) and complex averaging (blue)
profiles have intensities similar to each other, which indicates that the
complex averaging method does not appreciably reduce MS. At the depths of the
buried glass plate region (blue background region) in Fig. 3(g), we see that
the MFA and complex averaging methods show appreciably lower signal intensity
than that of the single acquisition method. It may be because of the
suppression of measurement noise by the averaging. At a depth close to the top
surface of the buried glass plate (indicated by a dashed blue line), the
complex averaging curve has a lifted “tail” with intensity higher than that in
the air. As the glass plate region is scattering-free, this tail might be due
not to scattering in the glass but to the MS in the superior scattering part
of the phantom. A similar result was reported by Yadlowsky et al.[8]. This
lifted tail is reduced when using the MFA method, which further supports the
MS-suppression ability of MFA. These findings suggest that the complex
averaging method reduces the system noise, and the MFA method reduces both
system noise and MS signals.

Figure 3(d)–(f) show the en face images at the same depth indicated by the
blue dashed lines in Figs. 3(a) and (g). The complex averaging image has
slightly better contrast of the scatterers than the single acquisition image,
but this contrast is not as good as that of MFA image (magnifications in the
yellow boxes). The MFA method provides better image contrast than the complex
averaging method because the former benefits from the MS reduction.

The SBR is computed from Fig. 3(d)–(f) for quantitative comparison. Ten
scatterers are selected to compute the signal intensity [red circles in Fig.
3(d)]. The background intensity is computed from the pixels in the glass plate
region [blue line in Fig. 3(d)]. The computed SBRs are 30.12, 32.13, and 34.64
dB for single acquisition, complex averaging, and MFA images, respectively.
Adopting the current protocol, the MFA method provides an SBR improvement of
4.5 dB relative to the single acquisition, and 2.5 dB relative to the complex
averaging method.

###  3.2 Postmortem zebrafish

Figure 4:  Cross sectional and en face images of the postmortem zebrafish. In
(a)–(c), the orange arrowheads denote the approximate focal positions, and the
orange boxes and magnification insets are supposed to show notochord
structure. The red arrowheads in (a)–(c) and (g)–(i) may indicate the
myosepta. In (d)–(f), yellow dashed areas denote the belly region, and blue
ellipses mark fine structures observed in the muscle region. The green
arrowheads in (d)–(f) and (j)–(l) may indicate the outer layer of the swim
bladder. In (g)–(i), yellow circles denote several hyper-scattering spots
observed in the belly region. In (j)–(l), green ellipses denote unidentified
structures that are better contrasted in the MFA image. At all depths, the MFA
method has the images with best contrast among the three methods.

Figure 4(a)–(c) shows the B-scans of the sample. The MFA images better
visualize structural boundaries [orange magnification inset in Fig. 4(c)],
which are smeared by strong noise in Figs. 4(a) and (b). This low-scattering-
intensity structure is expected to be the notochord according to its anatomic
features because a similar transparent structure of the notochord has been
observed at 5 dpf[20].

Figure 4(d)–(f) shows the en face images at depth Z1. The MFA image provides
better contrast for some fine structures in the muscle region (within blue
dashed ellipses) and some hyper-scattering spots in the belly region (yellow
dashed areas).

At depth Z2 [Fig. 4(g)–(i)], the contrast between the myosepta (red
arrowheads) and the surrounding muscle tissues is again highest for the MFA
method among the three methods. Some hyper-scattering spots in the belly
region (yellow dashed circles) are sharp and recognizable in the MFA image,
whereas they are blurred in the complex averaging image. These hyper-
scattering spots are almost unrecognizable in the single acquisition image.

At a deeper depth Z3, a thin hyper-scattering tissue is noted in the MFA image
[green arrowheads in Fig. 4(l)] but difficult to recognize in the single
acquisition and the complex averaging images [Figs. 4(j) and (k)]. This tissue
is also observed at the upper depth Z1 [green arrowheads in Fig. 4(d)–(f)]. A
similar structure has been visualized in adult zebrafish by OCT[21], and it is
expected to be the outer layer of the swim bladder. Unidentified structures
elongated along the dorsal side of the zebrafish [green ellipses in Fig.
4(j)–(l)] are also visualized with the highest contrast in the MFA image.

These findings suggest that the MFA method is applicable to biological
samples, reduces the MS signal, and enables better visualization of
structures.

##  4 Discussions

###  4.1 Validation of the induced defocus and its correction by MFA method

Figure 5:  (a) and (b) are B-scans at the same location of a phantom volume
without and with refocusing, respectively. (c) shows the FWHMs of the spatial
autocorrelation functions of the linear-scaled images, which are considered to
be proportional to the speckle size. The first and second rows present the
results without and with computational refocusing, respectively. Each column
presents the results for a volume measured with different focal positions. The
blue and red plots show the FWHMs for the fast-scanning (X) and slow-scanning
(Y) directions. The B-scans (a) and (b) were taken from the volumes
corresponding to the plots highlighted by green boxes. The used volumes are
identical to those used for the MFA images in Fig. 3.

To validate the defocus correction over multiple acquisitions, we compared the
lateral resolutions of a set of seven OCT volumes with different focal
positions. These volumes are the same phantom volumes used for the MFA method
demonstration in Section 3.1 (Fig. 3). Figures 5(a) and (b) show
representative B-scans of the phantom without and with refocusing,
respectively. Note that Fig. 5(b) is identical to Fig. 3(a).

The lateral resolution at each depth was evaluated by the speckle size of the
en face intensity image, which was estimated using a linear-intensity-based
auto-correlation. The speckle size was defined as the full width at half-
maximum (FWHM) of the auto-correlation function. This estimation was performed
along both fast-scanning (denoted as “X”) and slow-scanning (denoted as “Y”)
directions. A depth range of 1.1 mm (150 pixels) beneath the surface of the
scattering layer was used for this validation.

Without refocusing, the seven volumes with different defocus have speckle
sizes (X and Y FWHMs) varying along the depth [Fig. 5(c), first row]. After
applying refocusing to each volume, the FWHMs become almost constant at each
depth [Fig. 5(c), second row]. Note that the horizontal lines in the plot are
not regression or theoretical curves but plotted dots appearing as lines
because the FWHMs are highly constant. Most of the FWHMs converge on a value
of 14.6 $\muup$m, whereas the system optical lateral resolution (defined as
the focus radius where the intensity falls to $1/e^{2}$ of the maximum) is 18
$\muup$m. This suggests that computational refocusing corrects the defocus
over a depth range exceeding 1 mm, regardless of how much defocus is applied
in the measurement.

###  4.2 Focus shifting protocol optimization

####  4.2.1 Experimental optimization of MFA parameters

Table 1:  Summary of parameters used in the parameter-optimization experiment.
Parameters ${\Delta z}$, $N$, and $D$ are the focus shifting step, number of
averaged volumes, and overall focus shifting distance, respectively. The DOF
corresponds to 0.36 mm of the axial distance in air.

Measurement configuration |  ${\Delta z}$ [$\times$ DOF] | $N$ |  $D=(N-1){\Delta z}$ [$\times$ DOF]
---|---|---|---
#1 | 1 | 1, 2, 3, 4, 5 | 0, 1, 2, 3, 4
#2 | 1/2 | 1, 2, 3, …, 5, …, 9 | 0, 0.5, 1, …, 2, …, 4
#3 | 1/3 | 1, 2, 3,…, 13 | 0, 0.33, 0.36,…, 4
#4 | 1/4 | 1, 2, 3,…, 13 | 0, 0.25, 0.5,…, 3
#5 | 1/6 | 1, 2, 3,…, 13 | 0, 0.17, 0.33,…, 2

Figure 6:  SBR plot against the number of averaged volumes $N$ and the overall
focus shifting distance $D$. The type of plotted dot indicates the focus
shifting step ${\Delta z}$. The plot was produced to select optimal parameters
of $N$, $D$, and ${\Delta z}$. The measured sample is the phantom shown in
Fig. 3, and the SBR was computed by the process described in Section 3.1. The
black and green plots are projections showing the relationships between the
SBR and $D$ and between the SBR and $N$, respectively.

To find the optimal protocol, we experimentally explored the relationship
between the number of averaged volumes $N$, the overall focus shifting
distance $D$, the focus shifting step ${\Delta z}$, and the SBR. In this
experiment, we performed five sets of sequential volume measurements with
different ${\Delta z}$ as summarized in Table 1. For each set of measurements,
we extracted subsets to examine several configurations of $N$. $D$ is
determined from ${\Delta z}$ and $N$ as $D=(N-1){\Delta z}$. The sample used
for this measurement was the phantom shown in Fig. 3.

As there are two independent parameters, we plot the SBRs in a three-
dimensional space of the SBR, $D$, and $N$ as shown in Fig. 6. The color
gradient of the plotted points represents the SBR. The SBR takes a minimum
value when $N$ and $D$ are close to zero and increases with both $N$ and $D$.
The black plot is the projection along the $N$-axis that shows the
relationship between the SBR and $D$. We see that the SBR saturates at $D$
approximately twice the DOF and starts to drop at $D$ exceeding approximately
4 times the DOF. The green plot shows another projection along the $D$-axis
that gives the relationship between the SBR and $N$. It is seen that the
datasets with different ${\Delta z}$ saturate at different $N$. For example,
the datasets of configuration #1 (${\Delta z}=1\times\mathrm{DOF}$) and
configuration #4 (${\Delta z}=1/4\times\mathrm{DOF}$) saturate at
approximately $N$ = 4 and $N$ = 8, respectively. These results suggest that
both the parameters $N$ and $D$ contribute to the SBR improvement in the MFA
method.

Figure 7:  En face images of the zebrafish measured with different MFA
parameters at the same depth and the same lateral position. Among the four
sets of parameter configurations, configuration (c) gives best contrast. This
configuration was thus selected as the optimal parameter set in the present
study.

The optimal parameters were selected based on the above results and used in
the experiments described in the previous sections. We first selected $D$ as
twice the DOF, where the projection curve starts saturating (black plots in
Fig. 6). We then performed measurements of the zebrafish sample with ${\Delta
z}$ of 1, 1/2, 1/3, and 1/4. As $D$ was fixed to be twice the DOF, $N$ was 3,
5, 7, and 9, respectively. The imaging results (Fig. 7) show that the
parameters correpond to ${\Delta z}=1/3$, $D=2$, $N=7$ (third column) seems to
provide best contrast of the hyper-scattering spots in the belly region (see
the magnified insets within yellow boxes). We therefore selected this
parameter set as our optimal protocol and used it to perform the measurements
presented in the results section.

####  4.2.2 Interpretation of the SBR drop

In the protocol optimization, we note that the SBR drops when $D$ becomes too
large. A possible reason of this SBR drop is the uncorrected confocality;
i.e., the signal amplitude drops when the sample plane is far from the focal
plane. When $D$ exceeds a few DOFs, the focus shifts far from the depth of
interest. As a result, the amplitudes of both SS and MS signals decrease even
with the computational refocusing, whereas the amplitude of the measurement
noise is unchanged regardless of the focal position. The solution to overcome
this limit remains an open issue and requires further investigation.

###  4.3 Advantages of electrical tunable lens (ETL)

Table 2:  Comparison of examples of wavefront manipulation devices. Although
the ETL controls only defocus, it is the cheapest among the devices.

Device | ETL | Deformable mirror | SLM | Deformable membrane
---|---|---|---|---
Type | EL-10-30-CI-NIR-LD-MV | DM97-15 | PLUTO NIR2 | uDM2
Manufacturer | Optotune Switzerland | Alpao France | HoloEye Photonics Germany | Dyoptyka Ireland
Manipulable order of wavefront | Defocus only | Till 4-th order | Very high order | Not applicable
Response time | 15 ms | 1 ms (open-loop) |  $\sim$ 30 ms | <1 ms
Approximate price | $ 600 | $ 30,000 | $ 10,000 | $ 3,500

Figure 8:  Optical design schematics of the sample arm without (a) and with
(b) the ETL. (c) Color photograph of the sample arm with the ETL. GS:
galvanometric scanner, OBJ: objective, ETL: electrical tunable lens, and CM:
collimator. Since the ETL is a refractive optical element and easy to align,
switching between the two configurations without and with the ETL takes only a
few minutes.

Other wavefront manipulation devices have recently been used to reduce MS for
improving the OCT image contrast. One such device is a deformable mirror,
which was used by Liu et al.in aberration-diverse OCT. Here, the astigmatism
was introduced to the illumination beam, and its astigmatic angle was
modulated to decorrelate the MS among the measurements [13]. In addition, a
deformable membrane mirror has been combined with a full-field SS-OCT, and
fast three-dimensional volumetric cross-talk-free imaging has been
demonstrated [22, 23]. Borycki et al.used a spatial phase modulator for
digital aberration correction and MS reduction [24].

An advantage of the ETL over these wavefront manipulation devices is the low
cost of the ETL. A simple comparison of the aforementioned devices is shown in
Table 2. Another advantage of the ETL is that the ETL is easily integrated
with standard OCT systems because it is a refractive (not reflective) device.
The optical design schematics of the sample arm without and with the mounting
of the ETL in our system are shown in Fig. 8(a) and (b), respectively. Figure
8(c) is a photograph of our sample arm, where the ETL (shown by red dashed
box) was inserted between the collimator and the mount of the galvanometric
scanner. In our implementation, the replacement of the ETL and the realignment
of the system usually takes only a few minutes. This simple implementation of
the ETL makes our probe arm switchable between the standard OCT mode and MFA-
compatible mode.

###  4.4 Implementation in PS-OCT

PS-OCT is a functional extension of OCT that provides additional polarization
contrast for the quantitative evaluation of the optical properties of tissues
[25, 26, 27]. It is also known that MS affects the polarization measurement of
PS-OCT [26]. Adie et al.reported that the presence of MS randomized the
polarization states[28], leading to so-called depolarization. Several
researchers suspected that MS generates artifacts in accumulative or local
phase retardation measurements in the imaging of biological samples [29, 30].

However, the MS issue in PS-OCT has rarely been addressed. One example was
demonstrated by Gao et al.in FF-OCT, who used a Muller-matrix-based MS
subtraction process to remove the MS-induced local phase retardation[31].

The MFA method is applicable also to PS-OCT. Our preliminary study showed that
in a scattering phantom, the measured degree-of-polarization uniformity (DOPU)
decreases in deep regions, which can be mitigated by adopting the MFA
method[32]. We are currently working on adapting the MFA method to PS-OCT
imaging. This adaption may improve the accuracy of quantitative polarization
measurements in biological tissues.

###  4.5 In vivo imaging

The ETL enables fast and accurate focus modulation with a mechanical response
time of 15 ms [33]. However, the current MFA method captures multiple volumes
for averaging, which is time consuming and a challenging task for in vivo
imaging.

Nevertheless, the long acquisition time can be overcome by several means. One
possible solution is to use ultra-high-speed acquisition OCT[34, 35, 36] and a
motion tracking module to reduce the effect of sample motion. Another possible
solution is to use B-scan-based focus modulation instead of modulating the
focus for each volume. In other words, we can capture multiple B-scans at the
same location with different focal positions, applying a one-dimensional
version of computational refocusing [37], and average the focus-corrected
B-scans to generate an MS-reduced B-scan. Finally, the three-dimensional
volume can be obtained by repeating this process for a volume. In this case,
the acquisition time for generating an MS-reduced B-scan can be reduced if
compared with the current MFA method. A proof-of-concept experiment showed
that this method achieves a level of MS suppression similar to that achieved
by the current MFA method [38].

##  5 Conclusion

We have developed a new method termed “MFA”, to suppress the MS signals and
improve the image contrast in OCT. A scattering phantom was measured to
validate the MS reduction of the MFA method. The contrast improvement of the
MFA image was compared with single acquisition and complex averaging images
quantitatively. A postmortem zebrafish was measured to demonstrate the ability
of the MFA method in biological imaging. We expect this proposed method will
help better visualize the anatomic features. MFA method has great potential of
being adopted in in vivo imaging and other imaging modalities to reduce the
multiple-scattering effect.

## Appendix

##  Appendix A Deconvolution of the quadratic phase function

In Eq. (4), the deconvolution is represented as a convolution with the complex
conjugate of the quadratic phase function of the defocus amount, $\phi$. This
representation is adopted because the Fourier transform of the quadratic phase
function can be approximated by a phase-only function [37] and $\phi$ is an
even function. In this case, the deconvolution is expressed by the convolution
of the complex conjugate of $\phi$ as follows.

In general, the convolution of two functions $f(x)$ and $g(x)$ is expressed as
a multiplication in Fourier space as

$\mathcal{F}\left[{f(x)*g(x)}\right]=F(\nu)G(\nu),$ (5)

where $F(\nu)$ and $G(\nu)$ are the Fourier transforms of $f(x)$ and $g(x)$,
respectively. $\nu$ is the Fourier pair of $x$, and $x$ is not a specific
variable of space but a general variable.

Meanwhile, the deconvolution of $g(x)$ from $f(x)$ is expressed as the
multiplication $F(\nu)G^{-1}(\nu)$ in the Fourier space. When $G(\nu)$ is a
phase-only function, as in the case of the Fourier transform of the quadratic
phase function $\phi$, $G^{-1}(\nu)=G^{*}(\nu)$, where the superscript $*$
denotes the complex conjugate. Hence, the deconvolution of $g(x)$ from $f(x)$
can be written as the convolution $f(x)*g^{*}(-x)$, where we use a general
property of the Fourier transform, namely the inverse Fourier transform of
$G^{*}(\nu)$ is $g^{*}(-x)$.

When $g(x)$ is an even function, as in the case of the quadratic phase
function $\phi$, the deconvolution by $g(x)$ can be expressed as the
convolution with $g^{*}(x)$ as $f(x)*g^{*}(x)$. Hence, the refocusing can be
expressed as in Eq. (3).

## Acknowledgments

Lida Zhu is supported by the China Scholarship Council through the Chinese
Government Graduate Student Overseas Study Program. The authors also thank
Pradipta Mukherjee (University of Tsukuba) for providng careful review of the
manuscript.

## Funding

Core Research for Evolutional Science and Technology (JPMJCR2105); Japan
Science and Technology Agency (JPMJMI18G8); Japan Society for the Promotion of
Science (18H01893, 21H01836, 22K04962); China Scholarship Council
(201908130130); Austrian Science Fund (J4460, Schrödinger grand).

## Disclosures

L. Zhu, Makita, Lichtenegger, Lim, Y. Zhu, Yasuno: Yokogawa Electric Corp.
(F), Sky technology (F), Nikon (F), Kao Corp. (F), Topcon (F). Tamaoki,
Kobayashi: None.

## Data Availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * [1] S. Gigan, “Optical microscopy aims deep,” Nat. Photonics 11, 14–16 (2017).
  * [2] J. Olsen, J. Holmes, and G. B. E. Jemec, “Advances in optical coherence tomography in dermatology—a review,” Journal of Biomedical Optics 23, 040901–040901 (2018).
  * [3] C. Kut, K. L. Chaichana, J. Xi, S. M. Raza, X. Ye, E. R. McVeigh, F. J. Rodriguez, A. Quiñones-Hinojosa, and X. Li, “Detection of human brain cancer infiltration ex vivo and in vivo using quantitative optical coherence tomography,” Science Translational Medicine 7, 292ra100–292ra100 (2015).
  * [4] A. Lichtenegger, M. Muck, P. Eugui, D. J. Harper, M. Augustin, K. Leskovar, C. K. Hitzenberger, A. Woehrer, and B. Baumann, “Assessment of pathological features in alzheimer’s disease brain tissue with a large field-of-view visible-light optical coherence microscope,” Neurophotonics 5, 035002–035002 (2018).
  * [5] Y. Huang, S. Wang, Q. Guo, S. Kessel, I. Rubinoff, L. L.-Y. Chan, P. Li, Y. Liu, J. Qiu, and C. Zhou, “Optical coherence tomography detects necrotic regions and volumetrically quantifies multicellular tumor spheroids,” Cancer Res. 77, 6011–6020 (2017).
  * [6] I. A. El-Sadek, A. Miyazawa, L. T.-W. Shen, S. Makita, S. Fukuda, T. Yamashita, Y. Oka, P. Mukherjee, S. Matsusaka, T. Oshika, H. Kano, and Y. Yasuno, “Optical coherence tomography-based tissue dynamics imaging for longitudinal and drug response evaluation of tumor spheroids,” Biomed. Opt. Express 11, 6231–6248 (2020).
  * [7] I. A. El-Sadek, A. Miyazawa, L. T.-W. Shen, S. Makita, P. Mukherjee, A. Lichtenegger, S. Matsusaka, and Y. Yasuno, “Three-dimensional dynamics optical coherence tomography for tumor spheroid evaluation,” Biomed. Opt. Express 12, 6844–6863 (2021).
  * [8] M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple scattering in optical coherence microscopy,” Appl. Opt. 34, 5699–5707 (1995).
  * [9] C. Dunsby and P. M. W. French, “Techniques for depth-resolved imaging through turbid media including coherence-gated imaging,” J. Phys. D: Appl. Phys. 36, R207–R227 (2003).
  * [10] D. A. Boas, C. Pitris, and N. Ramanujam, _Handbook of Biomedical Optics_ (CRC Press, 2016).
  * [11] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Smart optical coherence tomography for ultra-deep imaging through highly scattering media,” Sci. Adv. 2, e1600370 (2016).
  * [12] M. Wojtkowski, P. Stremplewski, E. Auksorius, and D. Borycki, “Spatio-temporal optical coherence imaging — a new tool for in vivo microscopy,” Photonics Lett. Pol. 11, 44–49 (2019).
  * [13] S. Liu, M. R. E. Lamont, J. A. Mulligan, and S. G. Adie, “Aberration-diverse optical coherence tomography for suppression of multiple scattering and speckle,” Biomed. Opt. Express 9, 4919–4935 (2018).
  * [14] T. Ralston, D. Marks, F. Kamalabadi, and S. Boppart, “Deconvolution methods for mitigation of transverse blurring in optical coherence tomography,” IEEE Trans. Image Process. 14, 1254–1264 (2005).
  * [15] E. Li, S. Makita, Y.-J. Hong, D. Kasaragod, and Y. Yasuno, “Three-dimensional multi-contrast imaging of in vivo human skin by Jones matrix optical coherence tomography,” Biomed. Opt. Express 8, 1290–1305 (2017).
  * [16] A. Miyazawa, S. Makita, E. Li, K. Yamazaki, M. Kobayashi, S. Sakai, and Y. Yasuno, “Polarization-sensitive optical coherence elastography,” Biomed. Opt. Express 10, 5162–5181 (2019).
  * [17] K. Oikawa, D. Oida, S. Makita, and Y. Yasuno, “Bulk-phase-error correction for phase-sensitive signal processing of optical coherence tomography,” Biomed. Opt. Express 11, 5886–5902 (2020).
  * [18] L. Zhu, S. Makita, D. Oida, A. Miyazawa, K. Oikawa, P. Mukherjee, A. Lichtenegger, M. Distel, and Y. Yasuno, “Computational refocusing of Jones matrix polarization-sensitive optical coherence tomography and investigation of defocus-induced polarization artifacts,” Biomed. Opt. Express 13, 2975–2994 (2022).
  * [19] D. Hillmann, T. Bonin, C. Lührs, G. Franke, M. Hagen-Eggert, P. Koch, and G. Hüttmann, “Common approach for compensation of axial motion artifacts in swept-source OCT and dispersion in Fourier-domain OCT,” Opt. Express 20, 6761–6776 (2012).
  * [20] R. Haindl, A. J. Deloria, C. Sturtzel, H. Sattmann, W. Rohringer, B. Fischer, M. Andreana, A. Unterhuber, T. Schwerte, M. Distel, W. Drexler, R. Leitgeb, and M. Liu, “Functional optical coherence tomography and photoacoustic microscopy imaging for zebrafish larvae,” Biomed. Opt. Express 11, 2137–2151 (2020).
  * [21] D. Yang, M. Hu, M. Zhang, and Y. Liang, “High-resolution polarization-sensitive optical coherence tomography for zebrafish muscle imaging,” Biomed. Opt. Express 11, 5618 (2020).
  * [22] E. Auksorius, D. Borycki, and M. Wojtkowski, “Crosstalk-free volumetric in vivo imaging of a human retina with Fourier-domain full-field optical coherence tomography,” Biomed. Opt. Express 10, 6390–6407 (2019).
  * [23] P. Stremplewski, E. Auksorius, P. Wnuk, Ł. Kozoń, P. Garstecki, and M. Wojtkowski, “In vivo volumetric imaging by crosstalk-free full-field OCT,” Optica 6, 608–617 (2019).
  * [24] D. Borycki, E. Auksorius, P. Węgrzyn, P. Węgrzyn, and M. Wojtkowski, “Computational aberration correction in spatiotemporal optical coherence (STOC) imaging,” Opt. Lett. 45, 1293–1296 (2020).
  * [25] J. F. de Boer, C. K. Hitzenberger, and Y. Yasuno, “Polarization sensitive optical coherence tomography — a review [Invited],” Biomed. Opt. Express 8, 1838–1873 (2017).
  * [26] B. Baumann, “Polarization sensitive optical coherence tomography: A review of technology and applications,” Appl. Sci. 7, 474 (2017).
  * [27] Y. Yasuno, “Multi-contrast Jones-matrix optical coherence tomography — the concept, principle, implementation, and applications,” IEEE J. Sel. Top. Quantum Electron. 29, 1–18 (2023).
  * [28] S. G. Adie, T. R. Hillman, and D. D. Sampson, “Detection of multiple scattering in optical coherence tomography using the spatial distribution of Stokes vectors,” Opt. Express 15, 18033–18049 (2007).
  * [29] S. J. Matcher, C. P. Winlove, and S. V. Gangnus, “The collagen structure of bovine intervertebral disc studied using polarization-sensitive optical coherence tomography,” Phys. Med. Biol. 49, 1295–1306 (2004).
  * [30] A. Lichtenegger, P. Mukherjee, L. Zhu, R. Morishita, K. Tomita, D. Oida, K. Leskovar, I. A. El-Sadek, S. Makita, S. Kirchberger, M. Distel, B. Baumann, and Y. Yasuno, “Non-destructive characterization of adult zebrafish models using jones matrix optical coherence tomography,” Biomedical Optics Express 13, 2202–2223 (2022).
  * [31] W. Gao and S. Liu, “A method of removing retardance induced by scattering of collagenous fiber bundles,” J. Innov. Opt. Health Sci. 14, 2140003 (2020).
  * [32] L. Zhu, S. Makita, Y. Lim, J. Tamaoki, Y. Zhu, M. Kobayashi, and Y. Yasuno, “Multiple scattering suppression by multi-focus averaging in Jones-matrix optical coherence tomography,” in _Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XXVII,_ vol. 12367 J. A. Izatt and J. G. Fujimoto, eds. (SPIE, 2023), pp. 58–60.
  * [33] Optotune, “El-10-30-c lens,” https://www.optotune.com/el-10-30-c-lens. Accessed September 15, 2022.
  * [34] T. Klein, W. Wieser, L. Reznicek, A. Neubauer, A. Kampik, and R. Huber, “Multi-MHz retinal OCT,” Biomed. Opt. Express 4, 1890–1908 (2013).
  * [35] A. Kumar, L. M. Wurster, M. Salas, L. Ginner, W. Drexler, and R. A. Leitgeb, “In-vivo digital wavefront sensing using swept source OCT,” Biomed. Opt. Express 8, 3369–3382 (2017).
  * [36] J. P. Kolb, W. Draxinger, J. Klee, T. Pfeiffer, M. Eibl, T. Klein, W. Wieser, and R. Huber, “Live video rate volumetric OCT imaging of the retina with multi-MHz A-scan rates,” PLoS One 14, e0213144 (2019).
  * [37] Y. Yasuno, J.-i. Sugisaka, Y. Sando, Y. Nakamura, S. Makita, M. Itoh, and T. Yatagai, “Non-iterative numerical method for laterally superresolving Fourier domain optical coherence tomography,” Opt. Express 14, 1006–1020 (2006).
  * [38] Y. Zhu, L. Zhu, S. Makita, Y. Lim, and Y. Yasuno, “B-scan-wise multi-focus averaging method to suppress multiple scattering signal,” in _Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XXVII,_ vol. 12367 J. A. Izatt and J. G. Fujimoto, eds. (SPIE, 2023), pp. 175–178.