key: cord-0112137-6bxr1wmv
authors: Gao, Yuan; Jin, Guangzhen; Liu, Jian-Guo
title: Inbetweening auto-animation via Fokker-Planck dynamics and thresholding
date: 2020-05-18
journal: nan
DOI: nan
sha: e5a6774846de80474db96ace30d593ad1b9dc83d
doc_id: 112137
cord_uid: 6bxr1wmv

We propose an equilibrium-driven deformation algorithm (EDDA) to simulate the inbetweening transformations starting from an initial image to an equilibrium image, which covers images varying from a greyscale type to a colorful type on plane or manifold. The algorithm is based on Fokker-Planck dynamics on manifold, which automatically cooperates positivity, unconditional stability, mass conservation law, exponentially convergence and also the manifold structure suggested by dataset. The thresholding scheme is adapted for the sharp interface dynamics and is used to achieve the finite time convergence. Using EDDA, three challenging examples, (I) facial aging process, (II) coronavirus disease 2019 (COVID-19) invading/treatment process, and (III) continental evolution process are conducted efficiently.

Inbetweening auto-animation is to automatically generate animations (motions) given a starting and end images. The classical method for auto-animation use detailed kinematic equations for each object in the starting images, which is precise but time consuming due to case by case c.f. [CBE + 15, ZTG + 18].

Instead of analyzing the detailed kinematic equation for each object, we aim to propose an efficient and universal algorithm for inbetweening auto-animation based on Fokker-Planck dynamics on manifold and thresholding. We call this algorithm equilibrium-driven deformation algorithm (EDDA).

EDDA regards the end image as an equilibrium state of a Fokker-Planck equation and the inbetweening motion is driven by an underlying potential force determined by the equilibrium. This viewpoint is especially useful when the detailed physical process is not clear or hard to describe. For instance, the inbetweening motion of aging process, tumor growth, pneumonia invading for coronavirus disease 2019 or the formation of current continents/oceans starting from pangaea.

We first consider a Fokker-Planck equation in a flat domain Ω ⊂ R with a unique equilibrium π and no-flux boundary condition in Section 2 and then we propose an efficient solver for this Fokker-Planck equation in Section 3. The numerical solver for this part is based on structured grids and finite volume method [EGH00] . An unconditionally stable explicit time discretization is introduced, which automatically enjoys positivity, mass conservation law, exponentially convergence and also efficiency. For a Fokker-Planck equation on a closed manifold, we propose a similar efficient solver based on point clouds and the associated Voronoi tessellation in Section 4. The Voronoi tessellation automatically gives the manifold information and can be used to approximate surface gradient/divergence in the Fokker-Planck equation. Based on this, an analogue unconditionally stable explicit time discretization is introduced.

To realize the end image (the equilibrium) at finite time and the sharp dynamics of the inbetweening motion, we combine the explicit-time-discretization of the Fokker-Planck equation with the thresholding dynamics. When the equilibrium image has a sharp interface, the scheme adapting thresholding step converges faster than the purely Fokker-Planck iteration and the relative error reaches machine accuracy at finite time.

In Section 5, we apply EDDA proposed for either structured grids on Ω ⊂ R , or for pointclouds which suggests an underlining manifold to conduct three challenging and important examples: (I) facial aging process, (II) COVID-19 invading/treatment process, and (III) continental evolution process. In Example (I), inbetweening facial aging process at each time is simulated and potentially reveals the detailed changes of different part of human face over time. In Example (II), the inbetweening evolution of COVID-19 pneumonia invading before treatment and the fading away after treatment are simulated, which shows a good agreement with computerized tomography (CT) scans and also reveals promising application in the studying of pathology for COVID-19. In Example (III), the Fokker-Planck dynamics and thresholding are combined together to simulate the continental drifting process, which may suggest a new explanation for the formation of the current five continents of the world. From those examples from quite different research fields, EDDA are shown to be a very efficient and universal method with enormous potential applications in other fields of science and industry.

Suppose Ω ⊂ R is a closed subset in R . Assume the end image on Ω is described by a equilibrium density π(x) : Ω → R. The value of ρ indicates the gray level of the image for a grayscale image. In the case of Red-Green-Blue (RGB) image, we use three separate densities to indicate the RGB levels of the image separately. Then with π ∝ e −φ , the Fokker-Planck equation is given by

We consider the following natural no-flux boundary condition

(2.1) can be recast as the relative entropy formulation

See Section 4 for a Fokker-Planck equation on a d dimensional smooth closed submanifold of R . Now we state the ergodicity of Fokker-Planck equation (2.1). Assume

(2.5) π > 0, π ∈ C 1 (Ω).

Let L 2 (Ω; 1 π dx) be the weighted Sobolev space. Define the Fokker-Planck operator for (2.1) as L * : D(L * ) ⊂ L 2 (Ω; 1 π dx) → R with D(L * ) := {u ∈ H 2 (Ω; 1 π dx); ∂ n u π = 0 on ∂Ω} (2.6)

This can be regarded as the adjoint operator of the generator L = − 1 π ∇ · (π∇) of Fokker-Planck equation (2.1). One can see L * is self-adjoint operator in L 2 (Ω; 1 π dx) with compact resolvent (λI + L * ) −1 for λ large enough. Thus L * has only discrete spectrum without finite accumulation points. Furthermore, since π > 0, for ρ ∈ D(L),

Therefore, we conclude 0 is the simple principal eigenvalue of L * with the corresponding eigenfunction π, which leads to the spectral gap of L * in L 2 (Ω; 1 π dx), i.e. (2.8)

Thus due to (ρ − π) dx = 0, we have the following Poincare's inequality

Therefore, multiplying (2.1) by ρ π − 1, by (2.3) we have

which gives the ergodicity that

We present the numerical method based on structured grids for a Fokker-Planck equation on C ij = ((i − 1)∆x, i∆x) × ((j − 1)∆y, j∆y), i = 1, · · · , N, j = 1, · · · , M.

Then the cell centers (x i , y i ) are

)∆y, i = 1, · · · , N, j = 1, · · · , M.

We use ρ i,j to approximate the value of ρ(x i , y j ) and take π i,j = π(x i , y j ). Then the continuoustime finite volume scheme is

for i = 1, · · · , N, j = 1, · · · , M with the no-flux boundary condition (2.3). We assume the equilibrium π satisfies (3.4) π 0,j = π 1,j , π N +1,j = π N,j j = 1, · · · , M, π i,0 = π i,1 , π i,M +1 = π i,M i = 1, · · · , N, then the no-flux boundary condition (2.3) is reduced to (3.5) ρ 0,j = ρ 1,j , ρ N +1,j = ρ N,j j = 1, · · · , M,

Denote ρ k i,j as the value of ρ at t k = k∆t with time step size ∆t. Now we introduce an unconditionally stable explicit time discretization for (3.3)

for i = 1, · · · , N, j = 1, · · · , M with the no-flux boundary condition (3.5).

We now further simplify (3.6) as (3.7)

Then (3.7) can be rewritten as (3.9)

Denote h := max{∆x, ∆y}. From (3.9), we recast the scheme using a rescaled generator operator

Now we state the the positivity, maximal principle, mass conservation law and ergodicity of the scheme (3.7) as follows.

Proposition 3.1. Let π i,j = π(x i , y j ) > 0. Let ∆t be the time step and consider the explicit scheme (3.6) for the numerical solution ρ k i,j with (3.5). Assume the initial data ρ 0 > 0 satisfies

Then we have (i) positivity preserving property

(iii) the unconditional maximal principle for

(v) the exponential convergence

where µ 2 is the second eigenvalue of A defined in (3.26).

Proof. For (i), from (3.9), since π > 0, we know ρ k i,j > 0 implies ρ k+1 i,j > 0. To prove (ii), taking summation in (3.9), we have the (3.17)

The second term in the RHS of (3.17) is

where we used the no-flux boundary condition (3.4). Similarly, the third term in the RHS of (3.17) is

One can shift index for the last two terms in the RHS of (3.17) similarly. Therefore, using the no-flux boundary condition (3.5), we have the mass balance

To prove (iii), directly taking maximum in the RHS of (3.7) implies

which leads to (4.17).

To prove (iv), subtract (1+∆tλ i,j ) from both sides of (3.7) and then multiply by sgn ρ k+1

i,j π i,j − 1 . Thus using same argument with (iii), we have

which implies (4.18). Now we prove (v). Recall (3.10), i.e.

(3.23)

By shifting index and no-flux boundary condition (3.5) we have (3.24)

From this, one know

By the Perron-Frobenius theorem, µ 1 = 1 is the simple, principal eigenvalue of A with the ground state u * i,j ≡ 1 and other eigenvalues µ i of A satisfy |µ i | < µ 1 . Notice also the mass conservation for initial data u 0 = ρ 0 π satisfying (4.16), i.e., (3.27) u 0 − u * , u * (1+∆tλ)π = 0. Since also A is self-adjoint operator in the weighted space l 2 ((1 + ∆tλ)π), we can express u 0 using

c u , u is the eigenfunction corresponding to µ .

Therefore, we have

3.1. Thresholding for sharp dynamics. In this section, we combine the thresholding scheme with the Fokker-Planck dynamics to generate the inbetween motions with sharp interface, i.e., the density is described by linear combinations of two characteristic functions. In the computations later, one will see that the thresholding scheme also helps to achieve the finite time convergence to the sharp equilibrium density. Notice the dynamics of the Fokker-Planck equation is invariant when replacing ρ by cρ. Therefore, the initial density shall be adjust based on the mass conservation law (3.11). After this initial adjustment, assume initial data ρ 0 i,j ∈ {ρ 0 s , ρ 0 b }, which takes alternatively the value ρ 0 s , ρ 0 b . Assume the equilibrium is π i,j ∈ {π s , π b } which takes alternatively the value π s , π b . To combine the thresholding scheme with the Fokker-Planck dynamics, we need to choose the threshold ξ k at each step to conserve (3.13) as follows:

Step 1. Given ρ k i,j ∈ {π s , π b }, compute the explicit Fokker-Planck scheme (3.7) to updatẽ ρ k+1 i,j ∈ [π s , π b ] for any i = 1, 2, · · · , N and j = 1, 2, · · · , M .

Step 2. Choose threshold ξ k+1 and define (3.31) ρ k+1 i,j := π s χ {i,j;ρ k+1 i,j ≤ξ k+1 } + π b χ {i,j;ρ k+1 i,j >ξ k+1 } such that ρ k+1 satisfies (3.13).

In

Step 2, ξ k+1 can be found using bisection such that

4. EDDA based on point-clouds: Fokker-Plank equation on N Suppose (N , d N ) is a d dimensional smooth closed submanifold of R 3 . Assume the end image on N is described by a equilibrium density ρ ∞ (x) : N → R. Then the Fokker-Planck equation is given by

This can be recast as the relative entropy formulation

4.1. Construction of Voronoi tessellation and the upwind scheme on manifold N .

In this section, we construct an upwind scheme based on Voronoi tessellation for manifold N , which automatically gives a positive-preserving upwind scheme for the Fokker-Planck (4.1). Suppose (N , d N ) is a d dimensional smooth closed submanifold of R 3 and d N is induced by the Euclidean metric in R 3 . Q := {y i } n i=1 are sampled from the equilibrium density π = ρ N ∞ . Define the Voronoi cell as for any j = 1, · · · , n. If Γ ij = ∅ or i = j then we set |Γ ij | = 0. Let χ C i be the characteristic function such that χ C i = 1 for y ∈ C i and 0 otherwise. For i = 1, · · · , n,

is the piecewise constant probability distribution on N provided n i=1 ρ i |C i | = 1 and ρ i ≥ 0. Let π i be the approximated equilibrium density at y i satisfying n i=1 π i |C i | = 1. If ρ approx (y) = i ρ i χ C i (y) is an approximation of density ρ N (y), then ρ i is an approximation of the density ρ N (y i ).

Define the associated adjacent grids as Then using the finite volume method and the divergence theorem on manifold, we have

where n is the unit outward normal vector field on ∂C i . Based on this, we introduce the following upwind scheme. For i = 1, · · · , n,

We now interpret the upwind scheme as the forward equation for a Markov process with transition probability P ij (from j to i) and jump rate λ j

where (4.9)

|Γ ij |, i = 1, 2, · · · , n;

One can see it satisfies i P ij = 1 and the detailed balance property (4.10) P ij λ j π j |C j | = P ji λ i π i |C i |.

We refer to [GLW20] for the ergodicity of this Markov process.

In practice, instead of the |C i |, Γ ij in (4.8), one shall use the approximated coefficientsC i and Γ ij because we do not know the exact metric information of the manifold based only on point clouds. We omit the algorithm of finding the approximatedC i andΓ ij and refer to [GLW20] . 4.2. Unconditional stable explicit time stepping and exponential convergence. Now we propose an unconditionally stable explicit time discretization for the upwind scheme (4.7), which enjoys several good properties as the scheme (3.6), such as maximal principle, mass conservation law and exponential convergence.

Let ρ k i be the discrete density at discrete time k∆t. To achieve both stability and efficiency, we introduce the following unconditional stable explicit scheme We give the following proposition for several properties of scheme (4.11). The proof of this proposition is similar to Proposition 3.1 so we omit it.

Proposition 4.1. Let ∆t be the time step and consider the explicit scheme (4.11). Assume the initial data satisfies

Then we have (i) the conversational law for g k+1

(4.16)

(iv) the exponential convergence

where µ 2 is the second eigenvalue of I +∆tB (in terms of magnitude), i.e. µ 2 = 1−gapB∆t and gapB is the spectral gap ofB.

Thresholding for sharp dynamics. Assume the initial density is adjusted based on the mass conservation law (4.15). We now give the sharp dynamics by combining the Fokker-Planck equation on manifold with the thresholding scheme. Assume initial data ρ 0 i ∈ {ρ 0 s , ρ 0 b }, which takes alternatively the value ρ 0 s , ρ 0 b . Assume the equilibrium is π i ∈ {π s , π b } which takes alternatively the value π s , π b . Similar to Section 3.1, we choose the threshold ξ k at each step to conserve (4.15) as follows.

Step 1. Given ρ k i ∈ {π s , π b }, compute the explicit scheme (4.11) to updateρ k+1 i ∈ [π s , π b ] for any i = 1, 2, · · · , n.

Step 2. Choose threshold ξ k+1 and define such that ρ k+1 satisfies (4.15). Here ξ k+1 can be found using bisection such that 

In this section, three numerical examples are carried out to examine the capability and efficiency of the equilibrium-driven deformation algorithm (EDDA), which are the RGB colored facial aging transformation, the pneumonia of COVID-19 invading and fading away on CT scan images and the continental evolution process. 5.1. Example I: RGB colored facial aging transformation. In this example, we have two images with the same size in the RGB color model showing a ladys face at two different age, and employ the model to simulate the transformation from one image (initial) to another image (equilibrium), which will illustrate the facial aging process with time. The strategy is to define each image as three matrices, each matrix containing the value of a color mode (R or G or B). Then the transformation between the two images is completed by applying the inbetweening auto-animation three times based on Fokker-Planck dynamics (2.1).

The two images are extracted from [Posml] and are both 355 pixels in width and 575 pixels in height, which means a total of 204125 pixel points in each image. The initial image data is first adjusted to meet the mass conservation law (3.11). Time step ∆t is set to 0.01 and the total number of iterations is set to 10000 thus the final iteration time T = 100. The horizontal resolution ∆x and ∆y are both 10 −4 . We use the unconditional stable explicit time stepping scheme (3.6) and the no-flux boundary condition (3.5) to the Fokker-Planck equation (2.1) in domain Ω.

The relative root mean square errors (3.16) for the three color-modes are illustrated in Fig.  1 separately in semiology plot. Except for the different descend rates for the three colors, all simulated errors have the exponential convergence rates, which is consistent with the analysis in Proposition 3.1.

In order to see the transformation process between the two images, the images after iteration step 40, 100, 200, 400, 1000, 2000, 4000 and 10000 are shown and compared with the initial and the equilibrium images in Fig. 2 . The transformation process between two images are fast in the beginning (e.g. before step 200) and relatively slow after then. The transformation process in Fig. 2 clearly reveals the potential changes in different parts of the ladys face and hair with time. After 10000 steps of iterations, the updated image is nearly the same with the equilibrium except for the hair color.

Example II: COVID-19 pneumonia invading and fading away on CT scan images. In this section we will focus on an example based on the COVID-19 pneumonia invading and fading away process in a patients lung reflected on CT scan images and try to show the possible COVID-19 pneumonia growth dynamics with time before and after the treatment. In order to fulfill the task, two parts of simulations are presented. In the first part, two CT scan images taken on a patients lungs in the beginning (January 23th) and severe state (February 2nd) of the disease [ZL20] are selected to be the initial and equilibrium state, respectively; see Fig. 3 (left) . In the second part, two scan images at the severe state (February 2nd) and after a few-days treatment (February 9th) are selected to be the initial and equilibrium state, respectively; see Fig. 3 (right) . Each CT scan image can be represented by a gray scale image matrix thus the same method in Example I can be applied. The CT scan images are all cropped to 461 pixels in width and 370 pixels in height, which means a total of 170570 pixel points in each image. The time step ∆t is 0.01 and the total number of iterations is 6000 thus the final iteration time T = 60. The resolutions are ∆x = ∆y = 10 −4 .

After 6000 iterations, the relative root mean square errors from two parts of simulations both decrease with an exponential rate, as is shown in Fig. 4 . Moreover, the image evolution after step 20, 50, 100, 200, 500, 1000, 5000, 10000 (see Fig. 5 ) clearly demonstrate the pneumonia invading process into the patients lungs caused by COVID-19 in a few days (upper group of figures in Fig. 5 ) and the pneumonia fading away from the lungs after a stem cell treatment is applied to the patient (lower group of figures in Fig. 5 ), indicating a potential success of this treatment [ZL20] . We can further compare the evolution process with the real CT scan images taken on January 30th (see Fig. 3 ) and find satisfactory agreements, which indicates promising applications in this field. 5.3. Example III: Continental evolution process with thresholding for sharp dynamics. In this section, we try to reveal the evolution process of continentals in the world from Pangaea supercontinent (250 million years ago) to the current globe. In order to clearly distinguish the sharp dynamics evolving the continentals and the oceans, the thresholding scheme (4.20) and the Fokker-Planck dynamics (4.11) are combined to generate the inbetweening motions with the above sharp interfaces. The numerical experiment is carried out as follows.

Step (I). A group of points is selected on a unit sphere to be the dataset points. With the Centroidal Voronoi Tessellation (CVT) method on the unit sphere [Ren97, DGJ03] , the Voronoi cells on the unit sphere are generated and the locations of dataset points are adjusted accordingly to ensure the uniformity of these cells. Thus, the distributions of continentals and oceans derived from Pangaea period and current globes topography are described by two values (i.e. π s and π b ) on the Voronoi cells and are set to be the initial and equilibrium states, respectively. The Voronoi cell area C i , i = 1, · · · , n, with the total number of dataset points n, is computed and the Voronoi face Γ ij is determined by the geodesic length of the neighboring arc between cell i and j.

Step (II). Update the density at each point using the explicit scheme (4.11) linear Fokker-Planck equation.

Step (III). After several linear iteration steps, the threshold is selected following the steps in Section 4.3 and the thresholding scheme is applied to update the data.

The computations for Step (II) and

Step (III) will be looped until reaches the total iteration steps. Besides, a simulation case which only evolves the linear Fokker-Planck equation is carried out as the comparison.

For example III, we select a total of 3000 dataset points and generate the Voronoi cells on the unit sphere via the CVT approximation algorithm; see Fig. 6 The standard deviation for all the cell areas is 3.2 × 10 −4 , which means the nearly uniform distributions of data points on the sphere. The values at continental cells and the ocean cells are set to 0.9 and 0.1, respectively. The time step ∆t is set to 0.05 and the total number of linear iterations before the (k + 1)th thresholding adjustment is set to 2k, k = 1, · · · , N t , where N t = 50 is the times of the thresholding adjustments. The threshold ξ k is determined by bisection method such that (4.21) is satisfied. Here, the bisection domain limitation criterion is set to 10 −6 . The total number of iterations for the comparison simulation is set to be 10000. Table 1 . Comparison of the needed steps to reach the criterions between linear method and threshold method. Fig. 7 shows the temporal variations of the relative root mean square errors for the numerical example in the first 1200 time steps. The error from the thresholding method generally have a descend trend although with some abrupt increase due to the thresholding adjustments. The error decreases to nearly zero (less than the machine accuracy) after the 30th thresholding adjustment (a total of 960 time steps), which indicates the data is updated to the equilibrium. As a comparison, the error of the simulation via the linear method (red line in Fig. 7) , which leads an exponential convergence rate, is smaller than that from thresholding method before the 960th time step (black circle in Fig. 7) and is larger after then. In order to further compare the efficiency of the two methods, we calculated the total time steps needed for the error to reach the criterions and listed them in Table 1 . The comparisons clearly reveal the efficiency of the thresholding adjustment in the application of the sharp dynamics, especially when the criterion is small.

The continental evolution on the sphere after the 2th, 5th, 10th, 30th thresholding adjustment is illustrated and compared with the initial and equilibrium states in Fig. 8 . After several steps of thresholding adjustment, the sharp shapes of continentals quickly move from the initial Pangaea supercontinent towards the equilibrium state of current continentals. The distributions of continentals and oceans reaches the equilibrium state after the 28th thresholding adjustment, Figure 8 . The evolutions of continental movements on the unit sphere with the parameter ∆t = 0.05 and the total number of the thresholding adjustments is N t = 50. The continental evolution on the sphere after the 2th, 5th, 10th, 30th thresholding adjustment are illustrated and compared with the initial and equilibrium states. The black dots and polygons in each subplot illustrate the point clouds and the Voronoi cells, respectively. The orange and blue patch indicate the land and ocean, respectively. TH is short for thresholding step. The formation of the Antarctic is revealed at the bottom (southern part) of the globe (black arrow in TH 5). Note that the globes are shown in the same view angle so the Antarctic continental is out of view in the last two subplots. exactly the same as the current distributions. Although the evolution of the continental movements is simulated with the data-driven model, some potential dynamics of continental drifting such as the Antarctic formation can be noticed in the evolutions, which may contribute to the detailed explanation of the continental drifting theory.

We propose an efficient and universal equilibrium-driven deformation algorithm (EDDA) to simulate the inbetweening transformations given an initial and equilibrium. The algorithm automatically cooperates positivity, unconditional stability, mass conservation law, exponentially convergence and also the manifold structure suggested by dataset. Using EDDA, three challenging examples, (I) facial aging process, (II) COVID-19 invading/treatment process, and (III) continental evolution process are conducted efficiently. EDDA is shown to be a very efficient and universal method with enormous potential applications in other fields of science and industry.

Fully automatic generation of anatomical face simulation models

Voronoi-based finite volume methods, optimal voronoi meshes, and pdes on the sphere

Thierry Gallouët, and Raphaèle Herbin. Finite volume methods. Handbook of numerical analysis

Data-driven efficient solvers and predictions of conformational transitions for langevin dynamics on manifold in high dimensions

Can your face reveal how long you'll live?

Algorithm 772: Stripack: Delaunay triangulation and voronoi diagram on the surface of a sphere

Transplantation of ace mesenchymal stem cells improves the outcome of patients with covid-19 pneumonia