# Parallel implementations of random time algorithm for chemical network
stochastic simulations

Chuanbo Liu  Jin Wang  State Key Laboratory of Electroanalytical Chemistry,
Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Jilin,
People’s Republic of China  Department of Chemistry, Physics and Applied
Mathematics, State University of New York at Stony Brook, Stony Brook, USA

###### Abstract

In this study, we have developed a parallel version of the random time
simulation algorithm. Firstly, we gave a rigorous basis of the random time
description of the stochastic process of chemical reaction network time
evolution. And then we reviewed the random time simulation algorithm and gave
the implementations for the parallel version of next reaction random time
algorithm. The discussion of computational complexity suggested a factor of
$M$ (which is the connection number of the network) folds time consuming
reduction for random time simulation algorithm as compared to other exact
stochastic simulation algorithms, such as the Gillespie algorithm. For large-
scale system, such like the protein-protein interaction network, $M$ is on
order of $10^{8}$. We further demonstrate the power of random time simulation
with a GPGPU parallel implementation which achieved roughly 100 folds
acceleration as compared with CPU implementations. Therefore the stochastic
simulation method we developed here can be of great application value for
simulating time evolution process of large-scale network.

###### keywords:

Random time , Stochastic simulations , Parallel algorithm , GPGPU

††journal: arxiv

##  1 Introduction

Chemical reaction network time evolution is intrinsically stochastic. Starting
by Gillespie [1], a lot of algorithms had been developed for stochastic
simulation of the time evolution of chemical reaction network, such like the
First Reaction Method (FRM), Direct Method (DM). Some improved methods had
also been developed to accelerate the simulation processes. Such like the
Optimized Direct Method (ODM), computes $\sum a_{j}$ incrementally. The Next
Reaction Method (NRM), use one single random number for each simulation step.
Others like Sorting Direct Method (SDM), Logarithmic Direct Method (LDM),
Partial-propensity Stochastic Simulation Algorithm (PSSA), PSSA-composition
rejection (PSSA-CR), and Sorting Partial-propensity Direct Method (SPDM).
Focused on different aspects of the simulation steps.

Development of hardware, especially for the General Purpose Graphics
Processing Units (GPGPU), give the opportunity for developing parallel
algorithm that can take advantage of the large number of thread blocks (TBs).
Various data structure and implementations were also been developed for
acceleration of the stochastic simulation on fine-grid or coarse-grid level.
These implementations can achieve about 10 folds of simulation acceleration as
compared to mature CPU based simulation toolkit, such like StochKit. We here
demonstrated a new implementation of random time algorithm with GPGPU for
stochastic simulations. By carefully arranging the data distribution on global
and local memories we can accelerate the stochastic simulation with a data
parallel manner for roughly 100 folds. In this paper, we first briefly
reviewed the methods of random time simulation, then discussed the relative
computational complexity. At last, we demonstrated the methods with an
oscillation predator-prey model.

##  2 A brief introduction to the random time approach

In this part, we follow [2]. But we are trying to give a rigorous generalized
approach.

Consider $N\geqslant 1$ chemical species in a well-mixed system. Also consider
$M\geqslant 1$ reactions, labeled with $k$, the stoichiometric number of
$j$-th chemical species in $k$-th reaction is $\sigma_{jk}$. The simplest
stochastic models for describing the system is a continuous time Markov chain.
The state is represented by the molecule number of chemical species
$X(t)=\\{X_{1},X_{2},\ldots,X_{N}\\}$ and reactions modeled as possible
transitions of the chain. It was shown by Gillespie [6, 1] that in a well
mixed system, the probability for a specific reaction takes place is governed
by the propensity function $a_{k}(X(t),t)$ when consider time-inhomogeneous
chemical reaction networks. The propensity function is

$a_{k}(X(t),t)=c_{k}(t)h_{k}(X(t))$ (1)

where

$h_{k}(X(t))=\prod_{j}\binom{X_{k}(t)}{\sigma_{jk}}=\prod_{j}\frac{X_{k}(t)!}{\sigma_{jk}!(X_{k}(t)-\sigma_{jk})!}$
(2)

when considering the time evolution , and the state vector for $k$-th jump is
$v_{k}$. With Markov property, the conditional probability can be expressed
as:

$\mathbf{P}\\{\text{$k$-th reaction fired once in }(t,t+\Delta
t)|\mathcal{F}_{t}\\}=a_{k}(X(t),t)\Delta t$ (3)

where $\mathcal{F}_{t}$ is the $\sigma$-algebra representing the information
about the system that is available at time t [7]. The consequence of Eq. 3 is
the possibility of two reactions occurred at the same is on the order of
$(\Delta t)^{2}$, which means it is very unlikely to happen. This can be shown
by simply calculate the joint probability as

$\displaystyle\mathbf{P}\\{R_{k}(t+\Delta t)-R_{k}(t)\geqslant 2\\}$ (4)
$\displaystyle=\mathbf{P}\\{R_{k}(t+\Delta t)-R_{k}(t)\\}^{2}$
$\displaystyle\sim\mathbf{P}\\{R_{k}(t+\Delta t)-R_{k}(t)\geqslant
1,R_{j}(t+\Delta t)-R_{j}(t)\geqslant 1\\}$
$\displaystyle=\mathbf{P}\\{R_{k}(t+\Delta
t)-R_{k}(t)\\}\cdot\mathbf{P}\\{R_{j}(t+\Delta t)-R_{j}(t)\\}$
$\displaystyle\sim O((\Delta t)^{2})$

The system evolution can be described equivalently as a counting process by
replacing the random variables from chemical species molecule number to number
of times that a specific reaction fired. If $R_{k}(t)$ is the number of times
that the $k$-th reaction has fired up to time $t$, then the state at $t$ that
originated from $X(0)$ is given by

$X(t)=X(0)+\sum_{k=1}^{M}R_{k}(t)v_{k}$ (5)

This counting process is a generalized Poisson process $G(\lambda(t),t)$ in
the sense that the arrival rate $\lambda$ is a function of time $t$. The total
arrival times for generalized Poisson process $G$ can be calculated as

$R(t)=\int_{0}^{t}dG(\lambda(s),s)$ (6)

The Poisson distribution with parameter $\lambda$ in the time interval $(s,t]$
describe the arrival times of events of a Poisson process $N$ in this time
interval.

$\mathbf{P}\\{N(\lambda,s)-N(\lambda,t)=k\\}=\frac{(\lambda(s-t))^{k}}{k!}e^{-\lambda(s-t)}$
(7)

Therefore the probability of $k$ times events happened in the time interval
$(t,t+\Delta t]$ is

$\displaystyle\mathbf{P}\\{R(t+\Delta t)-R(t)=k|\mathcal{F}_{t}\\}$ (8)
$\displaystyle=\mathbf{P}\\{\int_{t}^{t+\Delta
t}dG(\lambda(s),s)=k|\mathcal{F}_{t}\\}$
$\displaystyle=\mathbf{P}\\{G(\lambda(t),t+\Delta
t)-G(\lambda(t),t)=k|\mathcal{F}_{t}\\}$
$\displaystyle=\frac{(\lambda(t)\Delta t)^{k}}{k!}e^{-\lambda(t)\Delta t}$
$\displaystyle=\frac{\left(\int_{t}^{t+\Delta
t}\lambda(s)ds\right)^{k}}{k!}e^{-\int_{t}^{t+\Delta t}\lambda(s)ds}$
$\displaystyle=\mathbf{P}\left\\{Y\left(\int_{t}^{t+\Delta
t}\lambda(s)ds\right)-Y\left(\int_{0}^{t}\lambda(s)ds\right)=k|\mathcal{F}_{t}\right\\}$

where $Y$ is unit rate Poisson process. The trick here is $\lambda$ is not
changed in the infinitesimal time interval $(t,t+\Delta t]$. Also it can be
noticed

$\lim_{\Delta t\rightarrow 0}\int_{t}^{t+\Delta t}\lambda(s)ds=\lim_{\Delta
t\rightarrow 0}\lambda(t)\Delta t$ (9)

With Eq. 7, the result is obvious. Eq. 8 suggests

$R(t)=Y\left(\int_{0}^{t}\lambda(s)ds\right)$ (10)

Next we can establish the relationship between $\lambda$ and the propensity
function $a_{k}$. We can now rewrite the left side of Eq. 3 as

$\displaystyle\mathbf{P}\left\\{R_{k}(t+\Delta
t)-R_{k}(t)=1|\mathcal{F}_{t}\right\\}$ (11)
$\displaystyle=\left(\int_{t}^{t+\Delta
t}\lambda_{k}(s)ds\right)e^{-\left(\lambda_{k}(s)ds\right)}$
$\displaystyle=\left(\int_{t}^{t+\Delta t}\lambda_{k}(s)ds\right)+O((\Delta
t)^{2})$

Compared to the right side of Eq. 3, we can have

$\int_{t}^{t+\Delta t}\lambda_{k}(s)ds=a_{k}(X(t),t)\Delta t$ (12)

which means,

$\lambda_{k}(s)=a_{k}(X(t),t)$ (13)

Therefore, from Eq. 5, Eq. 10 and Eq. 13, the system evolution equation is,

$X(t)=X(0)+\sum_{k=1}^{M}Y_{k}\left(\int_{0}^{t}a_{k}(X(s),s)ds\right)v_{k}$
(14)

So we have represented the chemical reaction network evolution process as an
increment counting process. Every infinitesimal time frame of this stochastic
process can be further decomposed to $M$ independence unit rate Poisson
processes. Eq. 14 can be rewritten to have the same formula with Eq. 5 by
introducing the “internal time” for each chemical reaction channel. The
internal time $T_{k}(t)$ is defined as

$T_{k}(t)=\int_{0}^{t}a_{k}(X(s),s)ds$ (15)

And Eq. 5 is becoming

$X(t)=X(0)+\sum_{k=1}^{M}Y_{k}\left(T_{k}(t)\right)v_{k}$ (16)

This is where the random time notion comes from.

##  3 Random time simulation algorithm

In order to perform a complete stochastic simulation, two things must be
cleared firstly,

  1. 1.

how much time passed before one of the stochastic processes, $Y_{k}$, fires;

  2. 2.

which $Y_{k}$ fires at that later time;

If we view the state of chemical reaction network as points in multi-
dimensional phase space, the time evolution of this model can be viewed as
random walked in spatial correlated non-homogeneous space. Then the two
questions above is the fundamental questions about space and time.

We will first determine the firing time problem. If we denote $Q(t,s)$ as the
probability of no reaction occurred in the time interval $(s,t]$, then by the
Markov property we have,

$Q(t,s+\Delta s)=Q(t,s)Q(s,s+\Delta s)$ (17)

According to Eq. 3, the probability of no reaction occurred in the time
interval $(s,s+\Delta s)$ is the prod of all independent reaction channels.

$Q(s,s+\Delta s)=\prod_{k=1}^{M}\left(1-a_{k}(X(s),s)\Delta s\right)$ (18)

Therefore,

$\displaystyle Q(t,s+\Delta s)$
$\displaystyle=Q(t,s)\prod_{k=1}^{M}\left(1-a_{k}(X(s),s)\Delta s\right)$ (19)
$\displaystyle=Q(t,s)\left(1-\sum_{k=1}^{M}a_{k}(X(s),s)\Delta s+O((\Delta
s)^{2})\right)$

By taking the limit of $\Delta s\rightarrow 0$, we arrive at the ordinary
differential equation,

$\frac{dQ(t,s)}{ds}=-Q(t,s)\sum_{k=1}^{M}a_{k}(X(s),s)$ (20)

Integral in the time interval $(t,s]$, and notice $a_{k}$ is not changed in
the time interval $(t,s]$, we have

$Q(t,s)=\exp\left(-\sum_{k=1}^{M}a_{k}(X(t),t)(s-t)\right)$ (21)

Eq. 21 tells us if a random number $r$ is uniformly distributed on $[0,1]$,
then the time interval is

$\Delta t=\ln(1/r)$ (22)

In the random time representation of the stochastic model of chemical reaction
network, for every infinitesimal time frame the stochastic process can be
described by a counting process and further be decomposed into unit rate
Poisson processes. All these unit rate Poisson processes are independent and
remain stationary until some chemical reaction channel is fired. For $k$-th
channel, follow the same argument from Eq. 17 to Eq. 20, the internal time has
the same distribution formula as Eq. 22.

$\Delta T_{k}(t)=\ln\left(1/r_{k}\right)$ (23)

where $r_{k}$ is a random number uniformly distributed on $[0,1]$. So if a set
of random number is given $\left\\{r_{1},r_{2},\ldots,r_{M}\right\\}$, the
corresponding real time for $k$-th chemical channel can be calculated from Eq.
15.

$\Delta T_{k}(t)=\int_{t}^{t+\Delta t_{k}}a_{k}(X(s),s)ds$ (24)

Since all the decomposed processes are unit rate Poisson processes, event is
expected to arrive in $\Delta T_{k}(t)$ for $k$-th chemical reaction channel.
The real time passed when $k$-th chemical reaction channel arrived is $\Delta
t_{k}$. Accordingly the first fired chemical channel is the one that has the
smallest internal time. First arriving wins all. And the real time interval
between two consequence reactions is

$\Delta t=min_{k}\\{\Delta t_{k}\\}$ (25)

Of course the firing chemical reaction channel is $k$-th channel. Therefore we
have answered the two main questions for stochastic simulations of the
chemical reaction network using random time representation.

##  4 Next reaction random time algorithm

According to the discussions in the former section, the simulation procedure
by applying the random time algorithm can be listed as Table 1.

Table 1: Simulation procedure of random time algorithm initialize time $t=0$
with state $X(0)=\\{X_{1}(0),X_{2}(0),\ldots,X_{N}(0)\\}$
---
while $(t<t_{max})$ do:
generate $M$ random numbers $r_{k}\sim U(0,1)$, where $k=1,\ldots,M$
set internal time interval $\Delta T_{k}=\ln(1/r_{k})$ for all $k=1,\ldots,M$
compute $\Delta t_{k}$ for all $k=1,\ldots,M$ by solving $\int_{t}^{t+\Delta
t_{k}}a_{k}(X(s),s)ds=\Delta T_{k}$
select reaction having $\Delta t=min_{k}\\{\Delta t_{k}\\}$, assume $j$-th
channel is triggered
update time $t=t+\Delta t$
update state $X(t+\Delta t)=X(t)+v_{j}$
end while

One thing about this simulation algorithm is it demands too many random
numbers. For every time step, it needs $M$ random numbers to determine the
time interval and firing reaction channel. Usually a stochastic simulation
must be running for a long time to collect sufficient data points for
statistical analysis. And the random number generators currently developed can
only give a pseudorandom number sequence which will fail more likely in long
sequences. The stochastic simulation algorithm that developed by Gillespie [6]
only demands two random numbers in each time step. If the next reaction
algorithm is used, only one random number is demanded in one simulation loop.

However, it is possible to develop an algorithm that only need one random
number for each simulation step. It can be noticed the only random number that
is consumed in one simulation step is the one that trigger the firing of a
specific reaction channel. Other random numbers are untouched. This means the
random numbers that do not trigger reaction channel firing can be reused. As a
matter of fact, this discussion suggests that the internal time of unfired
reaction channel is not changed after the end of this time frame. That is to
say, if the triggered reaction channel is the $j$-th channel, we can continue
the simulation by only refresh one internal time.

$\displaystyle\Delta T_{j}=\ln\left(1/r\right)$ (26) $\displaystyle\Delta
T_{k\neq j}=\Delta T_{k\neq j}$

where $r$ is a random number uniformly distributed in $[0,1]$.

This random time reuse algorithm is much like the Next Reaction Algorithm of
stochastic simulation, so it is called the Next Reaction Random Time Algorithm
(NRRTA). The simulation procedure of NRRTA is shown as Table 2.

Table 2: Simulation procedure of Next Reaction Random Time Algorithm

initialize time $t=0$ with state
$X(0)=\\{X_{1}(0),X_{2}(0),\ldots,X_{N}(0)\\}$
---
generate $M$ random numbers $r_{k}\sim U(0,1)$
set internal time intervals $\Delta T_{k}=\ln(1/r_{k})$ for all $k=1,\ldots,M$
while $(t<t_{max})$ do:
compute $\Delta t_{k}$ for all $k=1,\ldots,M$ by solving $\int_{t}^{t+\Delta
t_{k}}a_{k}(X(s),s)ds=\Delta T_{k}$
select reaction having $\Delta t=min_{k}\\{\Delta t_{k}\\}$, assume $j$-th
channel is triggered
update time $t=t+\Delta t$
update state $X(t+\Delta t)=X(t)+v_{j}$
generate a random number $r\sim U(0,1)$
refresh internal time interval: $\Delta T_{j}=\ln\left(1/r\right)$, $\Delta
T_{k\neq j}=\Delta T_{k\neq j}$
end while

Notes: the boundary conditions can be achieved by restrict the firing
channels. Those channels that are forbidden for firing is assigned a frozen
internal time.

##  5 Implementations

Table 3: Parallel next reaction random time algorithm. Initialize. Compute
$a_{k}$, set $\Delta t_{k}=ln(1/r_{k})/a_{k}$, set $p_{k}=a_{k}$
---
while $(t<t_{max})$ do:
Scan for $\Delta t_{\mu}=min_{k}\\{\Delta t_{k}\\}$
Update species for reaction $\mu$, set $t=t+\Delta t_{\mu}$
Recalculate $a_{k}$ for every reaction, set $\Delta t_{\mu}=ln(1/r)/a_{\mu}$,
set $\Delta t_{k}=p_{k}/a_{k}(\Delta t_{k}-\Delta t_{\mu})$
end while
Table 4: Model parser. // simple model of Lotka predator oscillation
---
# init
A = 1000
B = 1000
# reactions
1 A - -$>$ 2 A : 10
1 A 1 B - -$>$ 2 B : 0.01
1 B - -$>$ NULL : 10
# end
time = 2
steps = 10000

The most crucial and time consuming step with next reaction random time
algorithm is the selection of the minimum inner time $\Delta t_{\mu}$. For
memory efficiency and algorithm simplicity, the algorithm was modified as
outlined in Table 3.

The most accepted model file formate for current system biology is the SBML
(System Biology Markup Language). However, there is a lack of collective
graphic view of the model file, making the edition of these files not
convinient for ordinary use. We used a different model file formate which can
give both insight of the while network and also provide all the informations
for GSSA simulations. A typical predator oscillation system can be expressed
as in Table 4.

Since biology network follow the scale-free (or power-law) topography, the
reactant matrix, stochiometric matrix and species update matrix are
intrinsically sparse. Memory consuming is a potential issue for handling large
network which could consist of over $10^{2}$ nodes, in other words, over
$10^{4}$ reactions. Sparse matrix is the natural choice for biological network
simulation. For achieving $O(1)$ accessibiligy, we chose the COO (coordinate
formated matrix) formate for sparse matrix storage and computation.

In a large network, sequentially update of all the reactions after a firing
process is not a wise idea. Therefore the denpendency graph was proposed for
the efficient update of the network state. The idea was first published by
(Michael A. Gibson and Jehoshua Bruck, 2000, JPC). By following the network
connection through the affected species, a dependency graph can be constructed
as illustrated in the following graph. However, for massive parallel
computation, the network state is updated concurrently, so no dependency graph
is used in the cuda algorithm.

The key property for cuda programming is the memory arrangement, since memory
access latency is about 20-30 times longer than arithmetic operations.
Although memory access can be hidden partially by parallel execution of warps
by SMs, cooperate with the GPU memory caches is still vital. The utilizing of
registers and shared memory contribute to more efficient latency hidding for
parallel algorithm. For detail of this part, refer to (Professional CUDA C
Programming, 2014). In the original cuda programming set, computing hierarchy
is arranged as the grid-block 2 layer structure. While for natural setup for
parallel simulation of biology network, an intermediate layer called the
”bundle” is used for holding an independent trajectory. The reason for the
need of this intermediate layer is the limiting of the threads that could
reside in a particular block. Also for best performance, the threads per block
should better be $2^{n}$ for efficient reduction. Given this considerations,
the simulation architecture is arranged as the following.

  * 1.

1 trajectory corresponds to 1 bundle

  * 2.

bundle consists of a batch of blocks

  * 3.

1 thread in a block corresponds to 1 reaction related computation

The best pesudo-random number generator we discover is the Mersenne Twister
(MT) algorithm. With proper parameter values, MT can generate sequence with a
period as long as $2^{19,937}$ and extremely good statistical properties.
However, when implement MT in GPU architecture, the algorithm operation speed
was limited by the bandwidth of memory. For solving this issue, we had
implemented an modified version of the GPU MT algorithm. The generation speed
was increased but due to the limit of the size of shared memory, maximum 256
threads were allowed in a single block, and totally 256 blocks were allowed in
a single run. There are possibilities that the implementation can be improved,
but that involved with complete rewrite the algorithm in a self-made version
which demand much effort and time. So we changed the idea of a pure GPU
implementation to a heterogeneous version. The CPU were responsible for the
generation of random numbers and then passed to the global memory of GPU. This
implementation works fine for small trajectory number and short simulation
time. But when simulation trajectory number is large or simulation time is
very long, the global memory would be exhuasted. The situation can be solved
with different streams of GPU computation. Depending on the computation
latency of simulation, several streams can be performed to simulate different
trajectories while at the same time performing memory copy from CPU memory
heap to GPU global memory.

We find shared memory are vital for GPU implementation in various ways. One
problem of using shared memory is the memory size limitation. By using sparse
representation and limiting the biology network to maximum 2 reactants,
matrices were able to be fitted into the shared memory of each block. In most
biology network, the reactants is less than 2. Even for reactions involved
with 3 reactants, the reaction can also be transformed into an equivalence 2
reactants formate. With this setup, the loading from global memory can be
performed coasleced, also bank conflicing is avoided to some extent.

The most time consuming step in the GPU GSSA algorithm is the reduction step
for finding the next reaction id and fire time. For classic reduction
algorithm, the algorithm were optimized for $2^{n}$ elements reduction. For
biology network, the reaction number can be arbitrarily valued, therefore the
classic reduction algorithm should be modified to efficiently handle arbitrary
number of element. Another complexity is the reduction across different
blocks. Blocks within a bundle should be reduced together, while blocks belong
to different bundles should not be reduced. Therefore a two level reduction is
needed. With these considerations in mind, we developed an optimized version
for block reduction and this kernel can also be applied to bundle reduction as
well. This part forms the loop unroll template function which is the core part
of the reduction computation.

Figure 1:  Parallel next reaction random time stochastic simulation speedup of
predator-prey model. Speed up of parallel version of next reaction random time
stochastic simulation algorithm with predator-prey model. The computation was
done using a graphic card of Nvidia RTX 1060 with capability 6.1, global
memory size 6GB. 10000 trajectories were run simultaneously, with 100000 steps
for each trajectory.

The performance of the GPU version of GSSA depends on the capability of the
graphic card, also on the hierarchy computation units of cuda kernel. In our
current setup, when using Nvidia RTX 1060 with capability 6.1, global memory
size 6GB, about 400-500 folds speed-up can be achieved with single-precision
floating point number as shown in Figure 1. The speed up can be expected to
further increase when using a more powerful GPU.

##  6 Discussions

The major time consuming step for Gillespie’s stochastic simulation algorithm
is the calculation of the time interval between two consecutive reactions as
well as the searching for firing reaction channel. The searching step in
random time simulation algorithm is solved by finding the reaction channel
that has the minimal expected real time. However the time consumed by
calculating the time interval is not much accelerated by the random time
simulation algorithm as compared with Gillespie’s simulation algorithm. For a
time-inhomogeneous chemical reaction system the time interval for Gillespie’s
stochastic algorithm can be calculated from Eq. 1 and Eq. 22.

$\ln\left(1/r\right)=\sum_{k=1}^{M}\int_{t}^{t+\Delta
t}c_{k}(s)\prod_{j}\frac{X_{k}(t)!}{\sigma_{jk}!(X_{k}(t)-\sigma_{jk})!}ds$
(27)

While for random time stochastic simulation algorithm, the time interval can
be obtained by solving the following integral equation group.

$\displaystyle\ln\left(1/r_{k}\right)=\prod_{j}\frac{X_{k}(t)!}{\sigma_{jk}!(X_{k}(t)-\sigma_{jk})!}\int_{t}^{t+\Delta
t_{k}}c_{k}(s)ds,$ (28) $\displaystyle k=1,\ldots,M$ $\displaystyle\Delta
t=min_{k}\left\\{\Delta t_{k}\right\\}$

For large-scale chemical reaction network, the only hope for solving Eq. 27 is
through numerical integration. While for the random time simulation algorithm,
the time evolution is separated from the summation of propensity functions, so
integral can be calculated analytically for well behaved time-dependent
reaction rates.

If we assume the number of molecules in the system is on order of $N$, the
number of reaction channel is on order of $M$, the stoichiometric number is on
order of $H$, the number of reaction species for one single reaction is on
order of $Q$ and the integral steps is on order of $I$. In order to find
$\Delta t$, we apply binary tree searching to find the next reaction, the
element operation number is on order of $R$. So the number of element
operations for calculate the propensity function is on order of $2NQ+3Q$. And
the number of element operations for solving Eq. 27 is on order of
$\left[MI(2NQ+3Q)+M-1\right]R$. For Eq. 28, with well behaved time-dependent
perturbubation, the integral can be solved analytically, so the calculation is
reduced to a single element operation. The number of element operations for
solving Eq. 28 is on order of $[2NQ+3Q+2]M$. For large-scale chemical reaction
network, like $10^{4}$ chemical species, $M$ is on the order of $10^{8}$.
While if we consider a real biology chemical reaction network, $N$ is
typically on the order of $10^{3}$. The stoichiometric number for a biology
system is typically less than $10$, and reaction species involved in one
biochemical reaction is also less than $10$. The integral steps is related to
the calculation accuracy of the results since $I=\Delta t/\delta s$, here we
assume it to be on order of $10^{4}$. So the number of element operations for
solving Eq. 27 is on order of $10^{16}$, while for Eq. 28 the number is
reduced to an order of $10^{10}$.

But the accelerations can not be the same of computed above. The main reason
for not achieving the high speed include the memory limitation of the GPU
hardware, the limit number of the computation core and also the limitation of
total thread numbers. On the other side, for Eq. 27, the integration of
propensity function can be made parallel. In this circumstances, the summation
of contributions of propensity from all the reaction channels demands a
reduction algorithm, which followed by a searching process that locate the
next reaction. However, the situation is very different in solving Eq. 28.
Since no summation is demanded, the summation and searching are done within
one single reduction process. This means at the most time consuming step,
parallel random time algorithm can be approximate 2 times faster than parallel
Gillespie’s stochastic simulation algorithm. Hence random time based
algorithms are more appreciate for parallel computations.

##  7 ACKNOWLEDGEMENTS

Chuanbo Liu thanks supports by Natural Science Foundation of China,
No.32000888, Jin Wang thanks the supports from grant no. NSF-PHY 76066 and
NSF-CHE-1808474, Ministry of Science and Technology of China,
No.2016YFA0203200, Projects of Science and Technology Development, Jilin,
China, No.20180414005GH, Projects of Instrument and Equipment Development,
Chinese Academy of Sciences, No.Y928041001.

## References

  * [1] Gillespie,D.T. (1977) Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem., 81(25), 2340-2361.
  * [2] Anderson,D.F. and Kurtz,T.G. (2011) Continuous Time Markov Chain Models for Chemical Reaction Networks. In Koeppl,H., Setti,G., di Bernardo,M., Densmore,D. (eds), Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology. Springer New York, New York, NY, pp. 3-42.
  * [3] Dittamo, C., and Cangelosi, D. (2009). Optimized Parallel Implementation of Gillespie’s First Reaction Method on Graphics Processing Units. In 2009 International Conference on Computer Modeling and Simulation (pp. 156-161). IEEE. https://doi.org/10.1109/ICCMS.2009.42
  * [4] Komarov, I., and D’Souza, R. M. (2012). Accelerating the Gillespie Exact Stochastic Simulation Algorithm Using Hybrid Parallel Execution on Graphics Processing Units. PLoS ONE, 7(11), e46693. https://doi.org/10.1371/journal.pone.0046693
  * [5] Tian, H., and Burrage, K. (2005). Parallel implementation of stochastic simulation for large-scale cellular processes. In Proceedings - Eighth International Conference on High-Performance Computing in Asia-Pacific Region, HPC Asia 2005, pp. 621-626.
  * [6] Gillespie, Daniel T. (1992). A rigorous derivation of the chemical master equation. Physica A, 188:404-425.
  * [7] Kolmogorov, A. N. (1956). Foundations of the theory of probability. Chelsea Publishing Co., New York. Translation edited by Nathan Morrison, with an added bibliography by A. T. Bharucha-Reid.


# Adaptive Background Matting Using Background Matching

Jinlin Liu
<EMAIL_ADDRESS>
###### Abstract

Due to the difficulty of solving the matting problem, lots of methods use some
kinds of assistance to acquire high quality alpha matte. Green screen matting
methods rely on physical equipment. Trimap-based methods take manual
interactions as external input. Background-based methods require a pre-
captured, static background. The methods are not flexible and convenient
enough to use widely. Trimap-free methods are flexible but not stable in
complicated video applications. To be stable and flexible in real
applications, we propose an adaptive background matting method. The user first
captures their videos freely, moving the cameras. Then the user captures the
background video afterwards, roughly covering the previous captured regions.
We use dynamic background video instead of static background for accurate
matting. The proposed method is convenient to use in any scenes as the static
camera and background is no more the limitation. To achieve this goal, we use
background matching network to find the best-matched background frame by frame
from dynamic backgrounds. Then, robust semantic estimation network is used to
estimate the coarse alpha matte. Finally, we crop and zoom the target region
according to the coarse alpha matte, and estimate the final accurate alpha
matte. In experiments, the proposed method is able to perform comparably
against the state-of-the-art matting methods.

###### Abstract

Due to the difficulty of solving the matting problem, lots of methods use some
kinds of assistance to acquire high quality alpha matte. Green screen matting
methods rely on physical equipment. Trimap-based methods take manual
interactions as external input. Background-based methods require a pre-
captured, static background. The methods are not flexible and convenient
enough to use widely. Trimap-free methods are flexible but not stable in
complicated video applications. To be stable and flexible in real
applications, we propose an adaptive background matting method. The user first
captures their videos freely, moving the cameras. Then the user captures the
background video afterwards, roughly covering the previous captured regions.
We use dynamic background video instead of static background for accurate
matting. The proposed method is convenient to use in any scenes as the static
camera and background is no more the limitation. To achieve this goal, we use
background matching network to find the best-matched background frame by frame
from dynamic backgrounds. Then, robust semantic estimation network is used to
estimate the coarse alpha matte. Finally, we crop and zoom the target region
according to the coarse alpha matte, and estimate the final accurate alpha
matte. In experiments, the proposed method is able to perform comparably
against the state-of-the-art matting methods.

##  1 Introduction

The matting problem aims at estimating the precise alpha matte from the input
images, which is of great value in lots of image and video editing
applications. Considering the diversity of the real world scenes and the
complex details of the target objects, it is challenging to estimate accurate
alpha matte automatically. Mathematically, The matting formula is the
combination of the foreground $F$, background $B$ and alpha matte $\alpha$.
The given image $I$ can be calculated by the following method [32]:

$I_{z}=\alpha_{z}F_{z}+(1-\alpha_{z})B_{z}\,,\quad\alpha_{z}\in[0,1]\,.$ (1)

where $I$ is the input image, $F$ is the foreground, $B$ is the background,
$\alpha$ is the alpha matte, $z$ is the pixel of the input image $I$. In this
formula, only the input image $I$ is observed, all other variables are
unknown. As a result, solving such a problem is difficult and requires extra
information. Lots of methods have been proposed to achieve this goal. We
classify them to three categories: trimap-based, trimap-free and background-
based methods.

Trimap-based methods [33, 9, 23, 13, 31, 4, 30] refer to methods using trimap
as external input. Trimap contains background, foreground and unknown regions,
among which, only the unknown regions are meant to be estimated. In recent
years, trimap-based methods are capable of producing very accurate alpha matte
estimation results using deep networks. Whereas, the quality of trimap
generated by manual interaction in real applications may vary significantly,
which will lead to some unstable results. Also generating trimap is not easy
and can be time-consuming especially for video kind applications. It is nearly
not possible to generate trimap of consistent quality frame by frame. In
recent years, lots of trimap-free methods [22, 24, 34, 5, 28, 35] are proposed
to get rid of the limitations of trimap in real applications. The progress are
very promising. Trimap-free methods require no extra input and estimate
accurate alpha-matte by training on paired image to alpha matte datasets. Both
accurate semantic information and the detailed edge should be estimated
precisely by trimap-free methods. As a result, lots of training datasets are
required by these methods so as to perform well in complicated real world
applications [27]. In addition, for video kind applications, trimap-free
methods usually suffer from inconsistent estimation results among frames.
Background-based methods [27, 19] use extra background as input, and estimate
alpha matte from the input image and the given background. These kinds of
methods are very useful when applying to videos. The user only needs to
capture the background first. Background-based methods are able to produce
consist and stable results. In this scene, background-based methods are much
easier to use then trimap-based methods and are more stable then trimap-free
methods. Whereas, background-based methods require that the background should
not change during applications, which means that these methods will not work
when the background are dynamic or change frame by frame.

In summary, trimap-based methods require manual interaction and are nearly not
operable in video applications. Though trimap-free methods are convenient,
they require large amount of datasets and are difficult to generate stable
results in complicated scenes. Background-based method are stable and easy to
use in real applications, but they are limited to static backgrounds. Thus we
propose a new adaptive background matting method, which is stable, easy to use
regarding to all kinds of backgrounds.

In the proposed pipeline, the user first captures normal videos freely, moving
the camera as they want. Then, the foreground person steps out the scene, the
user captures the background according to the content of the captured video,
roughly covering the background regions. Differently from previously preposed
background-based methods, the proposed method does not require the background
to be static (demonstrated in Figure 1). Then we learn from the dynamic
backgrounds for accurate matting. The proposed method contains three parts:
background matching network (BMN), robust semantic estimation network (REN)
and accurate alpha estimation network (AEN). In the first step, for each frame
of the video, we identify the corresponding background frame in the background
frames using the background matching network. After that, the robust semantic
estimation network estimates the coarse alpha matte. With the coarse alpha
matte, we crop the input to remove blank regions, zoom it to a large size, and
use the alpha matte estimation network to estimate the final alpha matte.

The proposed method breaks the limitation of background-based method while
keeping the advantages, and is much more easily to use in real world
applications. In experiments, the proposed method is able to perform
comparably against the state-of-the-art matting methods. The main
contributions of this work are:

  * •

We propose a novel adaptive background matting method. To our best knowledge,
this is the first method that uses dynamic background video for accurate
matting.

  * •

We propose an effective pipeline to estimate alpha matte from the input video
and dynamic background video. BMN is designed to find the best matched
background from dynamic background frames. We estimate accurate alpha matte
progressively from coarse to fine using REN and AEN, making the proposed
method robust against possible background misregistration in reality.

Figure 1: The proposed adaptive background matting method uses dynamic
background video rather than static image for accurate matting. The camera is
not limited to be static, allowing the users to capture videos much more
freely. The user first captures their videos freely, moving the cameras as
they want. Then, the foreground person steps out the scene, the user captures
the background according to the content of the captured video, roughly
covering the background regions. The proposed method finds the best matched
background for each frame automatically.

##  2 Related work

#### Trimap-based methods.

Lots of traditional methods [7, 8, 10, 12, 15, 16, 26, 1, 2, 3, 6, 11, 17, 18,
29] take trimap as external input. Most of these methods solve 1 using some
prior assumptions. More recently, deep learning networks are widely used [33,
9, 23, 13, 31, 4, 30]. [33] first proposed deep image matting using
convolutional networks. [14] proposed to estimate the foreground and alpha
matte at the same time. Further, [30] incorporated semantic classification in
their matting networks. Lots of excellent works have been proposed recently
and the results are promissing. Whereas, the requirement of trimap limits the
applications of these methods. For video kind applications, it is nearly
impossible to generate trimap of consistent quality frame by frame.

#### Trimap-free methods.

Considering the limitation of trimap-based methods, lots of trimap-free
methods have been developed. Trimap-free methods require large training
dataset. Considering the difficulty of collecting matting datasets, most works
focus on human-related matting [22, 5, 28, 35]. [22] proposed to use coarse
dataset to boost the performance of human matting. [5] collected lots of human
dataset and proposed an end-to-end human matting method. Trimap-free methods
are convenient to use in real world applications, but they require large
amount of datasets and are difficult to generate stable results in complicated
scenes [27].

#### Background-based methods.

The aid of pre-captured background helps improve the quality of matting
significantly [27]. [27] proposed to use context switching block to combine
different inputs. In addition, they used unlabelled data to improve the
performance. Recently, [19] introduced a more efficient method to process high
resolution input in real time. As pointed out by [27], background-based method
requires the background to be static with tiny camera motion. This is not in
accord with the reality, as we usually move the camera to record more scenes.

Figure 2: The pipeline of the proposed method. The proposed method contains
three parts: background matching network (BMN), robust semantic estimation
network (REN) and accurate alpha estimation network (AEN). BMN aims at finding
the most matched background of the current frame from all background frames by
estimating the similarity values. In reality, the matched background from BMN
may exist misregistration, color and position difference with the current
frame. we estimate the alpha matte progressively. REN estimates the coarse
alpha matte. With the coarse alpha matte, we remove blank region first and
zoom the target area, then estimate the final alpha matte using AEN.

##  3 Proposed Approach

The proposed method aims at estimating alpha matte using dynamic backgrounds
and input videos. There are three parts in the proposed method: background
matching network (BMN), robust semantic estimation network (REN) and accurate
alpha estimation network (AEN). BMN identifies the corresponding background in
the background frames, trying to find the most matched background of the
current frame. In reality, the matched background from BMN may exist
misregistration, color and position difference with the current frame.
Therefore, we estimate the alpha matte progressively, estimating coarse alpha
matte first and then refining the coarse alpha matte. REN estimates the coarse
alpha matte from the input frame and the matched background from BMN. AEN
finally refines the coarse alpha matte using the coarse alpha matte, the input
frame and the matched background. The whole flowchart is demonstrated in
Figure 2.

###  3.1 Background Matching Network

In the proposed method, the background is dynamic and we have to identify the
corresponding background with the current input from all background frames.
Given the current input image or frame, the target is to find the best-matched
background. In real applications, the captured video and dynamic background
contain consecutive frames, and the difference among frames is small. Instead
of simply classifying the candidate background frame to matched or not matched
categories, we try to estimate the similarity value. The background matching
network takes two inputs: the input frame, and the background frame. The
network outputs the estimated similarity value of the two input images. The
structure of BMN is consisted of several convolutional layers, down-sampling
layers, and full-connected layers.

During the training process, we have the current frame, the corresponding
background, i.e. the ground truth background, and the ground truth alpha
matte. We use lots of miss-matched background and simple transformation of the
ground truth background to train BMN. We define the similarity value $S$ as
one minus the mean $L_{1}$ difference between the non-foreground region, i.e.
the background region,

$\displaystyle
S(x,y,\alpha)=1-\frac{\sum_{i}(|x_{i}-y_{i}|_{1}*(1-\alpha_{i}))}{\sum_{i}(1-\alpha_{i})},$
(2)

where $i$ represents the $i$th pixel, $x$ is the current frame, $y$ is the
input background frame, $\alpha$ is the corresponding alpha matte, and
therefore the non zero elements of $1-\alpha$ indicates the background region.

With the definition of the similarity value, we use the following loss
function to train BMN,

$\displaystyle\mathcal{L}_{BMN}=|B(x,y)-S(x,y_{gt},\alpha_{gt})|_{1},$ (3)

where $B$ represents the output of BMN, $y_{gt}$ is the ground truth
background, $\alpha_{gt}$ is the ground truth alpha matte.

During testing, for all background frames, we use BMN to calculate the
similarity value with the current frame. The background frame with the highest
similarity value is the best-matched background.

###  3.2 Robust Semantic Estimation Network

Though BEN finds the best-matched background with the current frame from all
background frames, it is possible that the found background is not perfectly
matched with the current frame. There might exist misregistration, color or
position difference in reality. Instead of estimating the alpha matte end-to-
end in one step, we choose a progressive way to estimate the alpha matte from
coarse to fine. In addition, the foreground region is usually much smaller
compared with the large background region in most videos. It is reasonable to
estimate a coarse alpha matte first, and then crop and zoom the target
foreground region accordingly so as to refine the details. Therefore, in the
first step, we use robust semantic estimation network (REN) to estimate the
coarse alpha matte.

The structure of REN is well displayed in Figure 2, which is basically a
feature pyramid network [21]. Low level features are integrated with the high
level features through skip connection. Multi-size features are employed to
predict the coarse alpha matte with different resolutions during training
process. As the target of REN is to estimate the coarse semantic information,
we use binary cross entropy loss to train REN.

$\displaystyle C(x,y)=-x*log(y)-(1-x)*log(1-y),$ (4)

where $x$ is the output and $y$ is the ground truth. With multi-scale alpha
matte output, we correspondingly calculate multi-scale loss,

$\displaystyle\mathcal{L}_{REN}=C(R(x),\alpha_{gt})+C(R^{\frac{1}{2}}(x),\alpha^{\frac{1}{2}}_{gt})+C(R^{\frac{1}{4}}(x),\alpha^{\frac{1}{4}}_{gt}),$
(5)

where $R(x)$ is the full size output of REN, $\alpha_{gt}$ is ground truth
alpha matte, $\frac{1}{2}$ and $\frac{1}{4}$ superscript represents the ratio
compared with the input size.

###  3.3 Accurate Alpha Estimation Network

We observe that the foreground region is usually a small part of the whole
image, and therefore we crop the target foreground region with slight board
and zoom it to a fixed large size ($640\times 640$ in experiments) to better
refine the details, including the input frame, the matched background, and the
coarse alpha matte. The cropping operation is based on the coarse alpha matte
estimated by REN. Accurate alpha estimation network (AEN) contains encoders
and decoders, low level features are concatenated with the high level features
so as to predict fine detailed alpha matte. Other than the alpha matte, we
simultaneously estimate the foreground residual similarly with [27, 19]. We
use $L_{1}$ loss to train AEN. The alpha matte loss is,

$\displaystyle\mathcal{L}_{\alpha}=|A(x)^{\alpha}-\alpha_{gt}|_{1},$ (6)

The foreground loss is,

$\displaystyle\mathcal{L}_{fgr}=|(A(x)^{fgr}+x)*\alpha_{gt}-{fgr}_{gt}*\alpha_{gt}|_{1},$
(7)

where $A(x)^{\alpha}$ is the output alpha matte, $A(x)^{fgr}$ is the output
foreground residual, $\alpha_{gt}$ is the ground truth alpha matte,$fgr_{gt}$
is the ground truth foreground. The final loss function is the sum of alpha
loss and foreground loss.

$\displaystyle\mathcal{L}_{AEN}=\mathcal{L}_{\alpha}+\mathcal{L}_{fgr}.$ (8)

###  3.4 Implementation details

#### Training details.

For speed consideration, both BMN and AEN are trained at a low resolution. The
input and the background are downscaled to $192\times 320$. When training AEN,
all images are cropped according to the coarse alpha matte and then zoomed to
$640\times 640$. BMN is trained independently with BMN and AEN. REN and AEN
are trained respectively for some epochs and then we co-train REN and AEN.

#### Background matching acceleration.

The proposed method manages to find the best-matched background from all
background frames by estimating the similarity of the current frame and all
candidate background frames. It can be very time-consuming when the background
video is long. Normally, the background videos contain consecutive frames, and
the difference among neighboring frames is small. Instead of finding the best
matched background from all background frames, we reduce the candidate
background frames first using fixed interval sampling. For example, if we set
the sampling interval to 8, only 12.5% of the background frames will be used,
accelerating the background matching process significantly.

Figure 3: The dynamic background matting dataset. Both the foreground and
background are dynamic. The background frame miss-matches the composed frame
frame-wisely. (a) and (b) well simulate the user capturing the video freely
and then capturing the background regions randomly in reality.

##  4 Dynamic Background Matting Dataset

Different from previous works, we aims to solve the dynamic background
problem, requiring the background to change frame by frame. To be more
realistic, the foreground object should move frame by frame at the same time,
to simulate the scene captured by a moving camera. VideoMatte240K [19] dataset
contains 482 high quality videos, with groundtruth alpha matte and foreground.
We collect 102 dynamic videos from the internet and manually compose the
videos from VideoMatte240K with the 102 collected dynamic videos to create our
dynamic background dataset. In our dataset, both the foreground and the
background change frame by frame. After composing the videos, we reverse the
background videos to serve as the captured backgrounds videos. In reality, the
captured background is random, and not matched with the freely captured video
frame by frame. This reversion operation ensures that the captured background
video miss-matches the composed video frame-wisely.

478 foreground videos from VideoMatte240K dataset, coupled with 82 collected
dynamic videos, compose our training dataset. Each video is combined with 20
background videos, making up 9560 dynamic videos. The left 4 foreground videos
and 12 dynamic background videos are used to create the testing dataset,
making up 48 testing videos, 22933 testing images. There is no overlap between
the training and testing dataset regarding to both the foreground videos and
the background videos. Some samples of our dataset is displayed in Figure 3.

Figure 4: Qualitative evaluation on real images from [19]. RVM [20] suffers
from inaccurate semantic estimation due to the lack of large training dataset.
FBA [9], RTBGM [19] and the proposed perform stably and output excellent hair
details. The proposed method performs better when the foreground is close to
the background.

z

##  5 Experiments

###  5.1 Evaluation results

#### Evaluation metrics.

The most widely used metrics to measure the matting results are mean square
error (MSE), sum of absolution difference (SAD), gradient and connectivity
error [33, 5]. The gradient and connectivity is able to measure the quality of
alpha matte similarly to human perception [25]. We calculate all four four
metrics for quantitative comparison purpose. As non-trimap-based methods need
to stimate per-pixel alpha matte, we calculate all metrics using the whole
image.

#### Baselines.

We categorize matting methods to three categories: trimap-based methods,
trimap-free methods, and background-based methods. Therefore, we compare the
proposed method with all three kinds methods, FBA [9], RVM [20] and RTBGM
[19]. FBA is the representative of trimap-based methods. RVM is one of the
latest trimap-free video matting methods. RTBGM represents the latest
background-based methods. As RTBGM requires static background, we test RTBGM
using ground truth backgrounds.

#### Performance comparison.

We compare all fours methods quantitatively on synthetic testing dataset and
qualitatively using real images from [19]. Using the synthetic testing
dataset, we calculate the evaluation metrics corresponding to all four methods
(Table 1). FBA performs the best with the given GT trimap. Though GT
background is used for RTBGM and the proposed method needs to identify proper
background from video frames, the metrics of RTBGM and the proposed method is
close, all of which are significantly better than trimap-free method RVM. In
Figure 4, we display the qualitative results on real images. The trimap for
FBA is generated manually using Photoshop. Generally, FBA, RTBGM and the
proposed method perform stably on all images and output excellent hair
details. Trimap-free method RVM suffers from inaccurate semantic estimation.
Trimap-free methods has no extra assistance and require large amount of
training dataset so as to estimate accurate semantic information. With respect
to the scene when the foreground is close to the background, FBA and RTBGM
wrongly segment part of foreground pixels to background. In comparison, the
proposed method output accurate alpha matte.

#### Running time.

All methods are tested using 720P images, Nvidia P100 graphic card. The
running time is listed in Table 1. FBA [9] takes 377$ms$, RVM [20] takes about
12.1$ms$, RTBGM [19] takes 24.1$ms$. For our method, the matting time
(REN+AEN) is 34.8ms, which is also very fast. Note that here we only compare
the matting time of the four methods. In reality, the background matching
network takes 2.83ms per-image to estimate the similarly value. So it will
take another $2.83n\ ms$ to find the best matched ground from $n$ background
frames.

Table 1: The quantitative results on testing dataset.

Method |  FBA [9] |  RVM [20] |  RTBGM [19] | ours
---|---|---|---|---
SAD($1e^{-3}$) | 4.04 | 5.84 | 4.46 | 4.45
MSE($1e^{-3}$) | 0.285 | 1.25 | 0.467 | 0.362
Gradient($1e^{-3}$) | 0.304 | 0.629 | 0.320 | 0.395
Connectivity($1e^{-3}$) | 0.979 | 2.623 | 1.19 | 1.23
Running time(ms) | 377 | 12.1 | 24.1 | 34.8 (REN+AEN)

Figure 5: Comparison of RTBGM [19] and the proposed method taking in the same
background with respect to different sampling intervals (background accuracy).
The difference between the best-matched background and the groundtruth
backgrounds increases with the growth of sampling intervals (from up to
bottom, the sampling interval is 1, 2, 4, 8, 16, 32, 64 respectively). Using
the same matched background, the proposed method performs more stably than
RTBGM. The performance of RTBGM becomes worse significantly when the sampling
interval is greater than 8 (the fifth row). Whereas, the proposed method
performs stably even when the sampling interval is 32 (the second row from the
bottom). Table 2: The quantitative results when using different background
sampling intervals.

Background sampling interval | 1 | 2 | 4 | 8 | 16 | 32 | 64
---|---|---|---|---|---|---|---
Difference between the matched and GT background | 2.13 | 4.40 | 6.87 | 10.04 | 13.71 | 18.01 | 23.88
Matched bk + RTBGM [19] |  SAD($1e^{-3}$) | 4.52 | 4.92 | 5.54 | 7.27 | 11.08 | 24.09 | 66.65
MSE($1e^{-3}$) | 0.499 | 0.849 | 1.41 | 3.03 | 6.63 | 19.2 | 61.1
Gradient($1e^{-3}$) | 0.336 | 0.376 | 0.420 | 0.496 | 0.609 | 0.807 | 1.15
Connectivity($1e^{-3}$) | 1.23 | 1.98 | 2.12 | 3.93 | 7.69 | 21.1 | 63.8
Matched bk + Ours |  SAD($1e^{-3}$) | 4.45 | 4.48 | 4.51 | 4.63 | 4.98 | 6.12 | 12.3
MSE($1e^{-3}$) | 0.362 | 0.382 | 0.411 | 0.520 | 0.833 | 1.93 | 7.92
Gradient($1e^{-3}$) | 0.395 | 0.401 | 0.407 | 0.419 | 0.439 | 0.476 | 0.580
Connectivity($1e^{-3}$) | 1.23 | 1.27 | 1.29 | 1.42 | 1.77 | 2.95 | 8.97

#### Ablation study.

The proposed method is composed of background matching and background matting.
To measure the quality of the proposed background matting method, we compare
RTBGM and the proposed using the same input background, i.e. the matched
background. In the proposed method, we propose to reduce the number of
candidate background frames by fixed interval sampling to accelerate the
background matching process. Considering lots of frames being removed by the
sampling strategy, it is possible that the found best-matched background is
not perfectly matched, existing some difference with the input frame. This
misregistration may impact the final performance. Therefore, we also test the
performance of the proposed method when the matched background exists
misregistration with the input frame.

In Figure 5, we display the results of RTBGM and the proposed method taking in
the same background with respect to different sampling intervals. In addition,
we draw the difference map between the best-matched background and the
groundtruth background. From the difference map, we can observe that when the
sampling interval is large, the found best matched background shows larger
difference with the groundtruth background. Even so, the proposed method
performs consistently well regarding to different sampling intervals,
benefiting from our robust semantic estimation network. In comparison, RTBGM
fails to perform well when the difference is large. In Table 2, we calculate
all four metrics to compare the results quantitatively. We also calculate mean
$L_{1}$ difference between the matched background and the groundtruth
background. The difference increases with the growth of sampling intervals.
The performance of RTBGM becomes worse significantly when the sampling
interval is greater than 8. In contrast, the proposed method performs stably
when the sampling interval is less than 32.

###  5.2 Applying to real world video

The greatest advantage of the proposed method is that the user now is not
limited to capture static background, which makes it much easier to use in
real world applications than previous background-based works. The user first
captures their videos freely, moving the cameras as they want. Then the user
captures the background, roughly covering the background regions. In Figure 6,
we demonstrate the results on real world videos. The background video lasts
for 12s, containing 344 frames. The sampling interval is set to 8. For
background matting process, we search from $344/8=43$ frames to find the best
matched background. The full inference time is $2.83\times 43+34.8=156ms$ per-
frame. The proposed method is able to find roughly matched background from
background frames. Though the best-matched background exhibits some
misregistration with the input frame, the results are very promising.

Also, we conduct user study using 35 real videos, comparing with RTBGM. The
proposed method is proved to be robust with background variations. Whereas
RTBGM is designed for static background and is sensitive to background
variations (Figure 5). Our method predominates the user study in experiments.
81% participants consider that the results of the proposed method are better
than RTBGM.

Figure 6: Applying to real world video. For fast background matting, the
sampling interval is set to 8. The proposed method is able to find roughly
matched background from background frames. Though the best-matched background
exhibits some misregistration with the input frame, the proposed method
outputs promising results.

##  6 Limitations and Discussions

The first step of the proposed method is to find the matched background frame
from the captured background video. Using Nvidia P100 graphic card, the
inference time of background matching network is about 2.83ms. Even so, it can
be time-consuming when the captured background video is long. If we set the
background sampling interval to 8, it takes seconds if we have to search for
the best matched background from more than 2826 frames (about 112s long
video). For applications where the running time is important, we recommend to
take background videos no more than 100 seconds. Whereas, if the running time
is not the limit, long videos are also applicable.

In the proposed pipeline, the user first captures normal videos, then the
background. There might be a consideration whether the users have to take the
background very carefully, trying to duplicate the video capturing process. As
we have shown in Figure 5, the proposed method is robust against background
misregistration. Therefore, the user only needs to cover the background region
roughly.

##  7 Conclusion

The limitations of background-based methods motivate us to propose a dynamic
background matting method. To achieve this goal, we use background matching
network to find the best-matched background frame by frame from dynamic
backgrounds. For acceleration purpose, we use fixed interval sampling to
reduce the background matching time. Then, robust semantic estimation network
is used to estimate the coarse alpha matte. Finally, we crop and zoom the
target region according to the coarse alpha matte, and estimate the final
accurate alpha matte. The proposed method performs stably and comparably
against the state-of-the-art methods. The proposed method allows the user to
move their camera when capturing videos in real applications. Considering the
limitations of static background matting methods (not flexible), green screen
matting (not convenient), trimap-based (not convenient) and trimap-free (not
stable) matting methods, the proposed method (stable, flexible and convenient)
probably reaches the best balance between convenience and matting quality.

## References

  * [1] Yağiz Aksoy, Tae-Hyun Oh, Sylvain Paris, Marc Pollefeys, and Wojciech Matusik.  Semantic soft segmentation.  ACM Transactions on Graphics (TOG), 37(4):72, 2018.
  * [2] Yagiz Aksoy, Tunc Ozan Aydin, and Marc Pollefeys.  Designing effective inter-pixel information flow for natural image matting.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 29–37. IEEE, 2017.
  * [3] Xue Bai and Guillermo Sapiro.  A geodesic framework for fast interactive image and video segmentation and matting.  In The IEEE International Conference on Computer Vision (ICCV), pages 1–8. IEEE, 2007.
  * [4] Shaofan Cai, Xiaoshuai Zhang, Haoqiang Fan, Haibin Huang, Jiangyu Liu, Jiaming Liu, Jiaying Liu, Jue Wang, and Jian Sun.  Disentangled image matting.  In The IEEE International Conference on Computer Vision (ICCV). IEEE, 2019.
  * [5] Quan Chen, Tiezheng Ge, Yanyu Xu, Zhiqiang Zhang, Xinxin Yang, and Kun Gai.  Semantic human matting.  In Proceedings of the 26th ACM international conference on Multimedia, pages 618–626. ACM, 2018.
  * [6] Qifeng Chen, Dingzeyu Li, and Chi-Keung Tang.  Knn matting.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 35(9):2175–2188, 2013.
  * [7] Yung-Yu Chuang, Brian Curless, David H Salesin, and Richard Szeliski.  A bayesian approach to digital matting.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 264–271. IEEE, 2001.
  * [8] Xiaoxue Feng, Xiaohui Liang, and Zili Zhang.  A cluster sampling method for image matting via sparse coding.  In The European Conference on Computer Vision (ECCV), pages 204–219. Springer, 2016.
  * [9] Marco Forte and François Pitié.  $f$, $b$, alpha matting.  arXiv preprint arXiv:2003.07711, 2020.
  * [10] Eduardo SL Gastal and Manuel M Oliveira.  Shared sampling for real-time alpha matting.  In Computer Graphics Forum, volume 29, pages 575–584. Wiley Online Library, 2010.
  * [11] Leo Grady, Thomas Schiwietz, Shmuel Aharon, and Rüdiger Westermann.  Random walks for interactive alpha-matting.  In Proceedings of VIIP, volume 2005, pages 423–429, 2005.
  * [12] Kaiming He, Christoph Rhemann, Carsten Rother, Xiaoou Tang, and Jian Sun.  A global sampling method for alpha matting.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2049–2056. IEEE, 2011.
  * [13] Qiqi Hou and Feng Liu.  Context-aware image matting for simultaneous foreground and alpha estimation.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4130–4139, 2019.
  * [14] Qiqi Hou and Feng Liu.  Context-aware image matting for simultaneous foreground and alpha estimation.  In The IEEE International Conference on Computer Vision (ICCV). IEEE, 2019.
  * [15] Jubin Johnson, Ehsan Shahrian Varnousfaderani, Hisham Cholakkal, and Deepu Rajan.  Sparse coding for alpha matting.  IEEE Transactions on Image Processing (TIP), 25(7):3032–3043, 2016\.
  * [16] Levent Karacan, Aykut Erdem, and Erkut Erdem.  Image matting with kl-divergence based sparse sampling.  In The IEEE International Conference on Computer Vision (ICCV), pages 424–432. IEEE, 2015.
  * [17] Anat Levin, Dani Lischinski, and Yair Weiss.  A closed-form solution to natural image matting.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 30(2):228–242, 2007.
  * [18] Anat Levin, Alex Rav-Acha, and Dani Lischinski.  Spectral matting.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 30(10):1699–1712, 2008.
  * [19] Shanchuan Lin, Andrey Ryabtsev, Soumyadip Sengupta, Brian L Curless, Steven M Seitz, and Ira Kemelmacher-Shlizerman.  Real-time high-resolution background matting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8762–8771, 2021.
  * [20] Shanchuan Lin, Linjie Yang, Imran Saleemi, and Soumyadip Sengupta.  Robust high-resolution video matting with temporal guidance.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 238–247, 2022.
  * [21] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2117–2125, 2017.
  * [22] Jinlin Liu, Yuan Yao, Wendi Hou, Miaomiao Cui, Xuansong Xie, Changshui Zhang, and Xian-sheng Hua.  Boosting semantic human matting with coarse annotations.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8563–8572, 2020.
  * [23] Hao Lu, Yutong Dai, Chunhua Shen, and Songcen Xu.  Indices matter: Learning to index for deep image matting.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 3266–3275, 2019.
  * [24] Yu Qiao, Yuhao Liu, Xin Yang, Dongsheng Zhou, Mingliang Xu, Qiang Zhang, and Xiaopeng Wei.  Attention-guided hierarchical structure aggregation for image matting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13676–13685, 2020.
  * [25] Christoph Rhemann, Carsten Rother, Jue Wang, Margrit Gelautz, Pushmeet Kohli, and Pamela Rott.  A perceptually motivated online benchmark for image matting.  In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 1826–1833. IEEE, 2009.
  * [26] Mark A Ruzon and Carlo Tomasi.  Alpha estimation in natural images.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), volume 1, pages 18–25. IEEE, 2000.
  * [27] Soumyadip Sengupta, Vivek Jayaram, Brian Curless, Steven M Seitz, and Ira Kemelmacher-Shlizerman.  Background matting: The world is your green screen.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2291–2300, 2020.
  * [28] Xiaoyong Shen, Xin Tao, Hongyun Gao, Chao Zhou, and Jiaya Jia.  Deep automatic portrait matting.  In The European Conference on Computer Vision (ECCV), pages 92–107. Springer, 2016.
  * [29] Jian Sun, Jiaya Jia, Chi-Keung Tang, and Heung-Yeung Shum.  Poisson matting.  In ACM Transactions on Graphics (ToG), volume 23, pages 315–321. ACM, 2004.
  * [30] Yanan Sun, Chi-Keung Tang, and Yu-Wing Tai.  Semantic image matting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11120–11129, 2021.
  * [31] Jingwei Tang, Yagiz Aksoy, Cengiz Oztireli, Markus Gross, and Tunc Ozan Aydin.  Learning-based sampling for natural image matting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 3055–3063, 2019.
  * [32] Jue Wang, Michael F Cohen, et al.  Image and video matting: a survey.  Foundations and Trends® in Computer Graphics and Vision, 3(2):97–175, 2008.
  * [33] Ning Xu, Brian Price, Scott Cohen, and Thomas Huang.  Deep image matting.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2970–2979, 2017.
  * [34] Yunke Zhang, Lixue Gong, Lubin Fan, Peiran Ren, Qixing Huang, Hujun Bao, and Weiwei Xu.  A late fusion cnn for digital matting.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR). IEEE, 2019.
  * [35] Bingke Zhu, Yingying Chen, Jinqiao Wang, Si Liu, Bo Zhang, and Ming Tang.  Fast deep matting for portrait animation on mobile phone.  In Proceedings of the 25th ACM international conference on Multimedia, pages 297–305. ACM, 2017.