Optimization algorithm of the velocity field determining in image processing

The paper proposes a new algorithm for constructing the velocity field, which is based on the study of the integral functional on the ensemble of trajectories. The resulting analytical representation of the variation of the integral functional gives us the gradient of the investigated functional. It allows to find the desired parameters using gradient methods, which determine the velocity field. This approach allows both optical flow and non-optical flow construction. The proposed algorithm can be used in the analysis of various images, in particular in radionuclide image processing.


Introduction
There is a large number of image processing algorithms corresponding to different purposes: image restoring, improving image quality, object recognition and motion analysis, contour detection etc. All these types of processing are important for radionuclide images, in particular image analysis based on the velocity field construction. Radionuclide image processing based on the velocity field determining was considered in papers [Kotina and Ovsyannikov, 2018;Kotina and Pasechnaya, 2014;Kotina and Ploskikh, 2012;Kotina, Ploskikh, and Babin, 2013].
The problem of the velocity field determining is considered in the literature by many authors in various formulations [Anandan, 1989;Barron, Fleet, and Beauchemin, 1994;Black and Anandan, 1996;Bruhn, Weickert, and Schnorr, 2005;Fleet and Weiss, 2005;Horn and Schunck, 1981;Sun D., Roth S., and Black, 2010]. The most famous formulation of the problem using the concept of optical flow assumes the constancy of brightness along the trajectories of the system under consideration [Horn and Schunck, 1981]. As an image property, not only the brightness but its gradient, Laplacian or Hessian are also considered [Papenberg et al, 2006]. In this formulation, functionals of quality (energy function) are constructed, which also include in addition spatial smoothness assumptions for the required velocity field [Horn and Schunck, 1981;Black and Anandan, 1996]. The task of minimization of the constructed functionals is reduced to solving the corresponding Euler-Lagrange equations. These equations can be reduced to sparse linear systems of large order, which are solved by block iterative methods. In works [Kotina and Ovsyannikov, 2018;Ovsyannikov and Kotina, 2012] the approaches for obtaining the velocity field in case of non-optical flow are suggested and corresponding Euler-Lagrange equations are obtained. In the work [Kotina, 2010], an optimization algorithm for constructing the velocity field based on discrete systems is proposed. This paper proposes a new algorithm for the velocity field construction, which is based on the variation of the integral functional in the problems of trajectory ensembles control presented in works [Ovsyannikov, 1990;Ovsyannikov, 2012]. It is assumed that the brightness (density) along the trajectories can change. The velocity field is considered as a function depending on the vector of unknown parameters. In this paper, illustrations of the work of the proposed algorithm are presented.

Problem Statement
Let us consider a system of differential equationṡ (1) Here t -time, x -spatial coordinate vector, x ∈ R n , u -parameter vector, u ∈ R r , fsufficiently smooth vector-function. We assume that the movement is described by the system (1). Let us introduce the distribution density ρ = ρ(t, x), which plays the role of the mass or charge density in various problems of mechanics and electrodynamics. In our case, it is quantitative characteristics of the image (brightness) depending on the spatial coordinates and time, or the intensity of radiopharmaceutical distribution in the processing of radionuclide studies [Kotina and Ovsyannikov, 2018;Kotina and Pasechnaya, 2014]. The equation which for a given vector function f (t, x, u) determines the change in the density function in space over time has the form [Ovsyannikov, 1990]: Initial condition for the equation (2) is where ρ 0 (x) -given function. We assume that function ρ = ρ(t, x) is given and we have to reconstruct function f (t, x, u), which determines the required velocity field. Here parameter vector u = (u 1 , u 2 , . . . , u r ) is considered as unknown.
We consider the problem of reconstruction of parameter vector u as an optimization problem. Therefore we use the methods of optimization of charged particle beam dynamics presented in works [Ovsyannikov, 1980;Ovsyannikov, 1990;Ovsyannikov, 2012]. Let M 0 ∈ R n be the set of initial values for the system (1). We assume that set M 0 is a closed set and has nonzero Lebesgue measure. The solutions of the system (1) we denote by The set of these solutions is called trajectory beam (or beam), coming from the set M 0 for a given parameter vector u. We denote by M t,u cross-section of the trajectory beam at time t for the fixed vector u: Let ρ 0 (x) -a known density (brightness) at the moment t = 0, which determines an image. Let us further assume that we know the density (brightness)ρ(x), characterizing changed image in time ∆t, here ∆t is small enough. Let us fix time T = ∆t. The problem is to find the vector of parameters u so that at the time T the density distribution calculated by the equation (2) with condition (3) would coincide with the density distributionρ(x), that is, Let us formulate the optimization problem. For that introduce a functional here M T,u -cross-section of trajectory beam for a moment T , g(x, ρ) -non-negative continuously differentiable on x and ρ function. For instance, as function g(x, ρ) the following function could be taken: whereρ(x) -given density in R n . Let us note, that the moment T is fixed here, but it also could be varied. We will solve the following problem: to find a minimum of the functional (7) when u ∈ U ⊂ R r , U is a compact convex subset of R r . Solving the problem of minimization of the functional (7) and determining the vector u, we solve the problem of function f (t, x, u) reconstruction, i.e. we determine the velocity field given by the formula (1).

Variation and gradient of functional
Let us follow the work [Ovsyannikov, 1990] and write down the variation of the functional (7) as here vector-function ψ(t, x) and function λ(t, x) are auxiliary functions, which satisfy along the trajectories of system (1) the following equations: with final conditions The determination of variation (9) is based on the use of the transformation of trajectory beams along crosssections [Ovsyannikov, 1980;Ovsyannikov, 1990;Ovsyannikov, 2012] using the variational equations for (1), (2). Let us assume that function f is differentiable with respect to u. Then, taking into account convexity of set U and using (9), we obtain an expression for the gradient of the functional (7) The complete output is presented in works [Ovsyannikov, 1980;Ovsyannikov, 1990]. The resulting expression (14) could be used for implementation of the optimization algorithm for determining velocity field. By using the functional gradient expression (14) it is possible to construct various methods of directed search of vector u. Let us note, that methods of parametric optimization and the optimal parameters obtaining on the basis of gradient techniques are used for many different problems [Andrianov and Edamenko, 2014;Aminov and Ovsyannikov, 2015;Shalymov, Fradkov, Liubchich, and Sokolov, 2017;Zavadskiy and Kiktenko, 2014].

The Velocity Field Determining Algorithm
Let us note that it follows from equation (2) that the derivative along the trajectories of system (1) satisfies equation Solution of the equation (15) with initial condition can be represented as Under images processing the form of the function f (t, x, u) is unknown. Therefore, the function f (t, x, u) can be considered as a function represented by a segment of some series, for example, of the Taylor series. The coefficients of the series are the required parameter vector u. In particular, at the first stage of constructing the velocity field we can consider function f (t, x, u) as a linear vector-function, i.e.
is a vector. Parameter vector u consists of the components of matrix A and vector C, i. e. u = (a 11 , a 12 , . . . , a nn , c 1 , . . . , c n ) * . The superscript * denotes the transposition of the vector. As before, we assume that u ∈ U , where U is a compact set. By obtaining u we determine the system (18) and the velocity field. Let us consider a particular case of n = 2 -the planar image case. In this case, the system (18) is a second order differential equation system. Let us write out the gradient of the functional (14) for this case: The equations (10), (11) will be the following: λ(T, x(T )) * = −(ρ(T, x(T )) −ρ(x(T ))) 2 + (ρ(T, x(T )) −ρ(x(T )))ρ(T, x(T )).
The equation (17) has a form Let us consider the iterative algorithm for determining the velocity field. We assume that M 0 , ρ 0 (x),ρ(x) are given and T (small enough) is fixed. The algorithm consists of the following main steps: 1. Define the vector of initial values for the unknown parameters u 0 = (u 0 1 , u 0 2 , . . . , u 0 n ). Let k = 0; 2. Start of k-th iteration of algorithm; 3. Calculate the value of integral functional J(u k ) using expression (7); 4. Obtain the values of the auxiliary functions ψ(t, x(t)) and λ(t, x, (t)) integrating equations (20) in reverse order from T to 0 under conditions (21), (22); 5. Calculate the gradient of the integral functional ∂J ∂u u=u k using the expression (14); 6. Obtain u k+1 with the gradient descent algorithm: , whereμ -parameter of gradient descent algorithm; 7. Check the stop conditions: achieving a given accuracy or a given number of iterations; 8. Let k = k + 1. Go to step 2.
The scheme of the algorithm for determining the velocity field is presented in Figure 1. The calculation of integrals in the described algorithm can be done as follows. Let us consider step 3. We take some finite partition where σ is a small positive value. For example, in a two-dimensional case, it would be a partition of the image into pixels. Let us represent the integral (7) in the form where g -integrand vector function calculated at the moment T at point x(T ) of the trajectory of the system (1) with a fixed parameter vector u k coming from the point x 0i ∈ e i , mes(x(T, e i , u k )) -measure of the set image e i according to the system (1) for a given parameter vector u k , at time T . For calculation of g ( it is necessary to integrate the system (1) from 0 to T for a given parameter vector u k with initial conditions x(0) = x 0i , i = 1, 2, . . . , N and also to find ρ(T, x(T, x 0i , u k )) according to the formula (17). The measure of the set x(T, e i , u k ) can be calculated taking into account the expression for the Jacobian of the transformation x(T ) = x(T, x 0 , u k ). Then we have Let us consider step 5. The integrals over the set M t,u in the formula (19) are calculated similarly to the calculation of integral in step 3 of the algorithm. Integral from 0 to T , due to the smallness of T , can be calculated by the method of middle rectangles. For the described algorithm any gradient descend method could be used. In this work the BFGS method (Broyden-Fletcher-Goldfarb-Shanno) [Nocedal and Wright, 2006] was used, which implementation is included in various mathematical packages. In this work the BFGS implementation from Octave package was used.
It should be noted that the described algorithm could be used in general case with different representations of function f (t, x, u). Linear representation was used as an example.

Implementation
In this section we consider test images and obtained in dynamic data acquisition mode radionuclide images. The radionuclide dynamic study is conducted as follows: the radiopharmaceutical is injected to the patient's body and then data collection begins. This mode Figure 2. The sequence of radionuclide images of hepatobiliary system study allow to observe the distribution of indicator (radiopharmaceutical) in studying organism system depending on time and spatial coordinates (2), i. e. radiopharmaceutical density distribution ρ = ρ(t, x), t ∈ [0, T ], x ∈ M , or taking into account the discrete nature of the obtained data we have the sequence of matrices ρ 1 (i, j), . . . , ρ N (i, j), i, j = 0, 1, . . . , n.
We will consider pairs of neighboring images of the sequence (24) and will define them as ρ 0 (i, j) andρ(i, j).
Note that for the analysis of the image data both the case of optical flow and the case of non-optical flow are interesting. Let us note that we have the optical flow, for example, if some motion of organ occurred and we have non-optical flow if there is a significant redistribution of the radiopharmaceutical in the considered system. By the algorithm proposed above we will determine velocity field at the nodes of a square grid with a step equal to the one pixel change in the distance along any axis.

Optical Flow
In case of div x f = 0 in the equation (15) we have optical flow, when it is assumed that the brightness remains constant along the trajectories.
Thus, if we put in the corresponding formulas above a 11 = −a 22 , then proposed algorithm could be considered as new algorithm for determining optical flow.  Figure 3 presents test images: initial image, initial image which is rotated relative to its center and contains some noise, obtained velocity field, which is shown in the initial image. Let us consider the work of the described algorithm on radionuclide images. There are images of the gallbladder in Figure 4, and we can see the shift of gallbladder from one frame to another. The first image shows the constructed velocity field, which determines the observed displacement.
In Figure 5 we also see images of the gallbladder, and we can see that some motion occurred. In the first image the obtained velocity field is drawn, which determines the observed motion. The figures in this subsection illustrate the results of the proposed algorithm work for the optical flow case.

Non-optical Flow
The case of a non-optical flow, when the brightness along the trajectories of the system changes, can be illustrated by examples of radionuclide images. We consider the image acquisition in the dynamic mode. This mode allows to observe the density distribution radiopharmaceutical depending on time and spatial coordinates in the considered organ or organism system [Kotina and Pasechnaya, 2014]. In Figure  6 we see a redistribution of the indicator (brightness) between two frames of the hepatobiliary system. There are pictures of a liver image with a visible redistribution of the radiopharmaceutical and the resulting velocity field drawn on the first image. In Figures 7, 8 we see examples of images of the gallbladder. In these frames, there are both the movement of the organ under investigation and the redistribution of the indicator (brightness) in it, the constructed velocity field reflects this.

Conclusion
The algorithm proposed in the article gives new possibilities for determining the velocity field, both for the case of an optical flow and for the case of a non-optical flow. The application of this algorithm can be useful in various areas of image processing, motion detection and correction, contouring and image analysis. This approach could have particular interest and development for the processing of radionuclide images, for obtaining both visual and quantitative information for the analysis of diagnostic images. It should also be noted that in this work we considered example of the linear a) b) Figure 5. a) Image of galbladder with velocity field. b) Image of galbladder after motion occured. approximation of the system (1). However, this approach allows us to consider any other approximations of the system (1). In this case, the gradient over all parameters specifying the approximation can be written out on the basis of formula (14). a) b) c) Figure 6. a) Initial image b) Image after redistribution of radiopharmaceutical c) Velocity field