An Optical Flow-Based Approach for Minimally Divergent Velocimetry Data Interpolation.

Three-dimensional (3D) biomedical image sets are often acquired with in-plane pixel spacings that are far less than the out-of-plane spacings between images. The resultant anisotropy, which can be detrimental in many applications, can be decreased using image interpolation. Optical flow and/or other registration-based interpolators have proven useful in such interpolation roles in the past. When acquired images are comprised of signals that describe the flow velocity of fluids, additional information is available to guide the interpolation process. In this paper, we present an optical-flow based framework for image interpolation that also minimizes resultant divergence in the interpolated data.

1. Introduction

Image interpolation is a fundamental problem encountered in many fields [1-9]. There are countless scenarios wherein images are acquired at resolutions that are suboptimal for the needs of specific applications. For example, biomedical images spanning a three-dimensional (3D) volume are often acquired with in-plane pixel spacings far less than the out-of-plane spacings between images. This can be the case with clinical images (e.g., from computed tomography (CT) and/or magnetic resonance (MR) imaging) as well as in vitro images acquired with modalities such as particle image velocimetry (PIV) [10-17]. In cases where motion estimation and registration are parts of an interpolation framework, hardware based approaches can offer solutions as well [18-24].

However, when acquired images are comprised of signals that describe the flow velocity of fluids, additional information is available to guide the interpolation process. Specifically, the flows of an incompressible fluid into and out of an interrogation volume must be equal according to conservation of mass [25]. Quantifying the deviation from zero net flow that is entering (or alternatively leaving) an interrogation volume (i.e., divergence) thus provides a means to direct interpolation in such a way as to reconstruct more physically accurate data.

Optical flow and/or other registration-based interpolators have proven useful in interpolating velocimetry data in the past [26-37]. Particle Image Velocimetry (PIV) is a technique that measures a velocity field in a fluid volume with the help of tracer particles in the fluid and specialized cameras [38, 39]. The default technique to determine the velocity field from the raw PIV data is a correlation analysis between two frames that were acquired by the cameras [40]. This technique can be extended to 3D as well. Optical flow-based approaches have been widely used in computer vision [41-44], and they have been appealing to researchers because of the flexibility of variational approaches. Regularizers can be used for different constraints in the energy functional to be minimized. In the conventional optical flow method there are two constraints, brightness and smoothness [45]. Optical flow-based methods have been promising in the area of fluid flow estimation in PIV [46-51]. For example, in [47], incompressibility of the flow is added as a constraint in the optical flow minimization problem. In [48], the vorticity transport equation, which describes the evolution of the fluid's vorticity over time, is used in physically consistent spatio-temporal regularization to estimate fluid motion.

Divergence and curl (vorticity) have been used in estimating optical flow previously [52-55]. In [52], the smoothness constraint is decomposed into two parts, divergence and vorticity, in this way, the smoothness properties of the optical flow can be tuned. In [56], both incompressibility and divergence-free constraints are used in the ill-posed minimization problem to calculate a 3D velocity field from 3D Cine CT images. In [54], a second-order div-curl spline smoothness condition is employed in order to compute a 3D motion field. In [55], a data term based on the continuity equation of fluid mechanics [25] and a second-order div-curl regularizer are employed to calculate fluid flow.

Here we present an optical-flow based framework for image interpolation that also minimizes resultant divergence in the interpolated data. That is, the divergence constraint attempts to minimize divergence in interpolated velocimetry data, not the divergence of the optical flow field. To our knowledge, using divergence in this way as a constraint in an optical-flow framework for image interpolation has not been investigated prior to the preliminary work presented in [57]. The method is applied to PIV, computational fluid dynamics (CFD), and analytical data and results indicate that the trade-off between minimizing errors in velocity magnitude values and errors in divergence can be managed such that both are decreased below levels observed for standard truncated sinc function-based interpolators, as well as pure optical flow-based interpolators. The proposed method thus has potential to provide an improved basis for interpolating velocimetry data in applications where isotropic flow velocity volumes are desirable, but out-of-plane data (i.e., data in different images spanning a 3D volume) cannot be resolved as highly as in-plane data.

The remainder of this paper is structured as follows. In Section 2, a definition of the term optical flow will be given and a canonical optical flow method will be briefly described. This will provide a basis for the following sections as most of the work described in this paper has been built on the described method. In Section 2, an optical flow-based framework for interpolating minimally divergent velocimetry data is described. The new method uses flow velocity data to guide the interpolation toward lesser divergence in the interpolated data. In Section 3, performance of the proposed technique is presented with experiments and simulations on real and analytical data. The results and performance of the proposed method are discussed and concluded in Section 4.

2. Methods

2.1. Optical Flow

Optical flow is the apparent motion of objects in image sequences that results from relative motion between the objects and the imaging perspective. In one canonical optical flow paper [45], two kinds of constraints are introduced in order to estimate the optical flow: the smoothness constraint and the brightness constancy constraint. In this section, we give a brief overview of the original optical flow algorithm (Horn-Schunck method) and the modified algorithm that was used in this project.

Optical flow methods estimate the motion between two consecutive image frames that were acquired at times t and t + δt. A flow vector for every pixel is calculated. The vectors represent approximations of image motion that are based in large part on local spatial derivatives. Since the flow velocity has two components, two constraints are needed to solve for it.

2.1.1. The Brightness Constancy Constraint

The brightness constancy constraint assumes that the brightness of a small area in the image remains constant as the area moves from image to image. Image brightness at the point (x, y) in the image at time t is denoted here as I(x, y, t). If the point moves by δx and δy in time δt, then according to the brightness constancy constraintThis can also be stated asIf we expand the left side of (2) with a Taylor series expansion, thenwhere the ellipsis (…) denotes higher order terms in the expansion. After canceling I(x, y, t) from both sides of the equationWe can divide this equation by δt, which leads toSubstitutingthe brightness constraint can be written in a more compact form:where I = ∂I/∂x, I = ∂I/∂y, and I = ∂I/∂t. In this form α and β represent the image velocity components and (I, I) represents the brightness gradients.

2.1.2. The Smoothness Constraint

Fortunately, points from an object that is imaged in temporally adjacent frames usually have similar velocities, which results in a smooth velocity field. Leveraging this property, we can express a reasonable smoothness constraint by minimizing the sums of squares of the Laplacians of the velocity components α and β. The Laplacians are

2.1.3. Minimization

Optical flow assumes constant brightness and smooth velocity over the whole image. The two constraints described above are used to formulate an energy functional to be minimized:Using variational calculus, the Euler-Lagrange equations can be determined for this problem. Those equations need to be solved for each pixel in the image. Iterative methods are suitable to solve the equations since it can be very costly to solve them simultaneously. The iterative equations that minimize (9) arewhere n denotes the iteration number and and denote neighborhood averages of α and β. More detailed information on the method can be found in [45].

2.2. Optical Flow with Divergence Constraint

2.2.1. Continuity Equation

According to the continuity equation in fluid dynamics, the rate of mass entering a system is equal to the rate of the mass leaving the system [25]. The differential form of the equation iswhere ρ is the fluid density, t is time, and is the velocity vector field. In the case of incompressible flow, ρ becomes constant and the continuity equation takes the form:This means that the divergence of the velocity field is zero in the case of incompressible flow. Figure 1 shows the change in flow velocity of a voxel.

Figure 1

Change in flow velocity of a sample voxel.

2.2.2. Symmetric Setup

For the new method, a symmetric interpolation setup is proposed as shown in Figure 2. In the figure, upper and lower slices are from the dataset and the interpolated slice is in the middle.In this section, I(x, y, t) denotes the velocity magnitude image and denotes the velocity vector components (i.e., V,V,V). If one approximates the expressions with Taylor expansion around the points (x, y), we getAfter substituting (14a) and (14b) into (13), terms can be arranged to obtain the new brightness constraint:

Figure 2

Illustration of the symmetric interpolation setup.

In the next step, we aim to minimize the divergence of the interpolated slice. Ideally, the divergence equation of the interpolated slice should be used:Since this information is unavailable, to generate the middle slice with as little divergence as possible, we can use the fact thatwhich leads to the following constraint by using the divergence expressions of the two outer slices, I(z − Δ) and I(z + Δ):Using Taylor expansion on (18) yieldsIn (19), we need the derivatives of V(z + Δ) and V(z − Δ) in the z-direction. Calculating these derivatives in the z-direction would require additional outer slices. To simplify this requirement, we can expand V(x + α, y + β, z + Δ) and V(x − α, y − β, z − Δ) around the points (x, y, z) and obtain the following:

Using (20a) and (20b) in (18), we obtain the new divergence constraint that does not require additional slices for the z-direction derivative, Combining (15), (21) and the optical flow smoothness constraint, we obtain the new energy functional that needs to be minimized,where Using variational calculus, the Euler-Lagrange equations can be determined for this problem. They need to be solved for each pixel in the image. The iterative equations that minimize the solutions are given bywhere n denotes the iteration number and and denote neighborhood averages of α and β. The coefficient expressions in (24a) and (24b) are given asThe numerical scheme to solve the Euler-Lagrange equations is based on the solution laid out in [45]. More detailed information on the steps of the derivation can be found in the appendix.

There have been several studies that attempt to improve the performance of optical flow techniques and computation schemes [41, 44, 58–64]. For example, in [59] non-linear convex penalty functions are used for the constraints in the optical flow energy functional. The approach uses numerical approximations to obtain a sparse linear system of equations from the highly nonlinear Euler-Lagrange equations. The resulting linear system of equations can be solved with numerical methods like Gauss-Seidel, which is similar to Jacobi method, or successive over-relaxation (SOR), which is a Gauss-Seidel variant. Another improvement to variational optical flow computation is presented in [60]. The approach uses a multigrid numerical optimization method and because of its speedup gains, it can be used in real-time. After all these advances, in [58], it was argued that the typical formulation of optical flow has changed little and most of the advances have been mainly numerical optimization and implementation techniques and robustness functions. This is also true for the proposed method as well. The optical flow portion of this interpolation framework is closely related to the Horn-Schunck method. The derived numerical scheme to solve the equations enhances this notion while its implementation is straightforward and simple. For example, setting the divergence coefficient γ to 0 in (24a) and (24b) reduces the solutions to Horn-Schunck solutions. The numerical scheme is also sufficient for velocimetry data because unlike in many other types of images, stark discontinuities are unexpected in velocimetry images at Reynolds numbers on the order of biomedical flows.

2.3. PIV Setup

The testing datasets were acquired using particle image velocimetry, an optical experimental flow measurement technique. PIV data acquisition and processing generally consists of the following steps: (1) computational modeling, (2) physical model construction, (3) particle image acquisition, (4) PIV processing, and (5) data analysis. The testing datasets were acquired for an in-vitro model of a cerebral aneurysm. Patient-specific computed tomography (CT) images were first segmented and reconstructed to obtain the computational cerebral aneurysm model as shown in Figure 3. The computational model was then translated into an optically clear, rigid urethane model using a lost-core manufacturing methodology. The physical model was connected to a flow loop consisting of a blood analog solution seeded with 8 μm fluorescent microspheres. Fluid flow through the physical model was controlled at specific flow rates (3, 4, and 5 mL/s). PIV was performed using a FlowMaster 3D Stereo PIV system (LaVision, Ypsilanti, MI), where the fluorescent particles were illuminated with a 532 nm dual-pulsed Nd:YAG laser at a controlled rate, while two CCD cameras captured the images across seven parallel planes (or slices) within the aneurysmal volume. A distance of 1 mm separated the planes. Two hundred image pairs, at each flow rate and slice, were acquired at 5 Hz. The image pairs were processed using a recursive cross-correlation algorithm using Davis software (LaVision, Ypsilanti, MI) to calculate the velocity vectors within region of interest (i.e., the aneurysm). Initial and final interrogation window sizes of 32 by 32 pixels and 16 by 16 pixels, respectively, were used. Detailed explanation of the experimental process can be found in [65]. A sample experimental model is shown in Figure 4.

Figure 3

Dimensions of the aneurysm.

Figure 4

Example flow slice from the PIV experiments.

The proposed algorithm was developed in MATLAB (Mathworks, Inc). Since the proposed algorithm has two separate terms for divergence and smoothness, different combinations of coefficients can be used for the terms. However, in order to get a clear idea about the performance of the method only one set of parameters were used in the simulations. The divergence term's coefficient γ was set to 150. From previous tests, it was seen that the proposed method performed better when a relatively large γ was used while keeping the smoothness coefficient λ small. The smoothness coefficient λ was set to 1. The same smoothness coefficient was also used for the Horn-Schunck based method. The iterations for both methods were set to 2000. Each PIV dataset used in testing had 7 slices. The slices were originally 154x121. They were cropped and zero-padded to reach 128x128. The size of the region where MSE and divergence were calculated is 110x110. Even though there are 7 slices in each dataset, only 3 slices were reconstructed from the datasets. These are slices 3, 4 and 5. Two different spacing steps were used between the slices. The first one is Δz=2 where the neighboring slices z-1 and z+1 were used to reconstruct the middle slice. The second one is Δz=4 where slices z-2 and z+2 were used for the interpolation; e.g., slices 1 and 5 were used to reconstruct slice 3. The method was tested against linear interpolation and an implementation of Horn-Schunck optical flow based interpolation.

2.4. Analytical Datasets

The method was tested with a 3D divergence-free analytical dataset and a CFD data set with turbulent flow. The analytical dataset is given below:Out-of-plane distance was kept much higher than the in-plane resolution. In order to assess the robustness of the proposed method, each velocity field was perturbed by Gaussian noise. The noise had zero mean and standard deviation of 10% of the maximum velocity in each velocity field.

2.5. Computational Fluid Dynamics (CFD) Simulations

The original computational aneurysm model was imported into ANSYS ICEM (ANSYS, Canonsburg, PA), where the inlet and outlets of the aneurysm model were extruded. After meshing was performed to discretize blood volumes into tetrahedrons, the final mesh was imported into ANSYS Fluent where the blood volume was modeled as an incompressible fluid with the same density and viscosity as the blood analog solution used in experiments. The vessel wall was assumed to be rigid, and a no-slip boundary condition was applied at the walls. A steady flat 4 ml/s flow profile was applied at the inlet of the model, and zero pressure boundary conditions were imposed at the outlets. The overall CFD approach has been described previously in [65, 66].

3. Results

Figure 5 shows divergence and MSE comparison graphs when Δz=2. The proposed method consistently achieves lower divergence values than the Horn-Schunck-based interpolation whereas the MSE values vary between better and worse values. On average, divergence values were 11% lower than the Horn-Schunck-based interpolation. In some cases, the proposed method achieves up to 20% lower divergence values.

Figure 5

Divergence and MSE comparisons when slice distance is 2 mm.

Figure 6 shows divergence and MSE comparison graphs when Δz=4. In this case, the proposed method consistently achieves lower divergence and MSE values than the other tested methods.

Figure 6

Divergence and MSE comparisons when slice distance is 4 mm.

Figure 7 shows original, noisy, and interpolated slices from the analytical dataset for comparison. In the figure, only V and V components were plotted to show the effect of the divergence term. In Figure 8, it can be seen that the proposed algorithm reduces divergence while the MSE is increased in the CFD dataset.

Figure 7

Plotted V and V components of the 3D analytical divergence-free vector field. (a) Original, (b) Gaussian noise added, (c) linear interpolation, (d) Horn-Schunck based interpolation, and (e) proposed method. Note that the proposed method is able to achieve a smoother velocity field in the corners of the interpolated data.

Figure 8

Divergence and MSE comparisons for the CFD dataset.

The graphs in Figure 9 show the behavior of the proposed method as the divergence coefficient γ increases linearly. In this simulation, the smoothness coefficient λ was kept constant (λ=1). The graphs are taken from the PIV dataset. The divergence graph profiles were consistent across different images and three datasets. The MSE graph profiles may differ slightly from the divergence graph profiles across different datasets, but MSE always increased with increasing γ. The coefficient values tested were from 0 to 2000. The profiles shown in the figure show that there needs to be a balance between the divergence and the smoothing terms. The graphs in the figure are consistent with profiles of other published ℓ2-based regularization methods [67, 68]. Figure 10 shows the behavior of the proposed method as γ and λ increase linearly.

Figure 9

Divergence and MSE profiles of the proposed method as γ is increased linearly while λ = 1.

Figure 10

Divergence and MSE profiles of the proposed method as γ and λ are increased linearly from 0.1 to 2500.

The computational cost of obtaining flow vectors with the proposed method is similar to that of the Horn-Schunck approach. Even though the iterative solutions of the proposed method employ several terms, these need to be computed only once and can be reduced to a simpler form that is similar to the Horn-Schunck solutions. Both approaches took around 0.1 seconds to obtain an optical flow field on a single core of an Intel dual core CPU (i7-6500U @ 2.50GHz).

Another parameter that could affect the divergence, MSE, and computational cost was the iteration number. We ran simulations with different iteration numbers and noted that the divergence and MSE results seem to stabilize after 200 iterations. Higher iteration number mostly had an effect on the computation time.

4. Discussion and Conclusions

A new optical flow-based framework for image interpolation that also reduces divergence is proposed. The new method uses flow velocity data to guide the interpolation toward lesser divergence in the interpolated data. In addition to the symmetric interpolation setup, the method introduces a new divergence term into the canonical optical flow method. The method is applied to PIV, analytical, and CFD data. The method was tested against linear interpolation and the Horn-Schunck optical flow method since it uses a similar formulation as the Horn-Schunck method. The proposed method applies a symmetric interpolation setup and considers a new divergence term in addition to the brightness and smoothness terms in the energy functional.

In order to test the effects of the divergence term, both the Horn-Schunck and proposed methods were subject to the same smoothness coefficient. When tested on the noisy analytical data, the proposed method achieved a smoother and less noisy interpolated velocity field.

The proposed method was also applied to the PIV data with different values of smoothness and divergence term coefficients, α and γ, respectively. Results indicate that the tradeoff between minimizing errors in velocity magnitude values and errors in divergence can be managed such that both are decreased below levels observed for standard truncated sinc function-based interpolators as well as pure optical flow-based interpolators. The divergence term coefficient, γ, needs to be large enough to reduce divergence in the interpolated data but not so large as to dominate the energy functional and introduce errors into the final interpolated velocity field. The effect of the iteration number on the divergence and MSE numbers was found to be minimal after 200 iterations. The computational cost of the method was similar to that of the Horn-Schunck based approach.

The method uses a numerical scheme that is well-known and straightforward. It is true that a more recent optical flow computation scheme may lead to performance gains in quality and/or speed-up. Methods presented in [59, 60] have become popular because of their speed, simplicity, and flexibility. Adoptation of recent numerical optimization and implementation techniques will be explored for future research.

The proposed method has potential to improve the interpolation of velocimetry data when it is difficult achieve an out-of-plane resolution close to the in-plane resolution. The results also indicate that the effect of the new divergence term in the optical flow functional can be appreciated better as the distance between the interpolated slice and the neighboring slices increases. It was noted that the proposed method outperforms the tested methods in both divergence and MSE values when the slice distance was increased. When the slice distance is small, the proposed method achieves lower divergence than the other methods while achieving similar MSE values.