A distributed resistance inverse method for flow obstacle identification from internal velocity measurements

We present a penalization parameter method for obstacle identification in an incompressible fluid flow for a modified version of the Oseen equations. The proposed method consists in adding a high resistance potential to the system such that some subset of its boundary support represents the obstacle. This allows to work in a fixed domain and highly simplify the solution of the inverse problem via some suitable cost functional. Existence of minimizers and first and second order optimality conditions are derived through the differentiability of the solutions of the Oseen equation with respect to the potential. Finally, several numerical experiments using Navier–Stokes flow illustrate the applicability of the method, for the localization of a bi-dimensional cardiac valve from MRI and ultrasound flow type imaging data.


Introduction
The pathway of blood flow through the heart is regulated by four membrane structures or valves, opening and closing by differential blood pressures. Valvular heart diseases affect 20% of the population. Among them, aortic valve stenosis is the most prevalent one in developed countries [22].
Imaging the 3D shape of the valves is a challenging task. For instance, the aortic valve is very thin (0.5 mm), and therefore its shape can be imaged nowadays using only two modalities: computerized tomography (CT) and transesophageal echocardiography (TEE). Since CT is based on x-rays, it is only used in patients that are selected for valvular replacement in order to obtain the aortic root dimensions for sizing the prosthesis. Such CT images are usually obtained when the valve is closed. Obtaining the image at open valve position requires about 5 times larger radiation dose since the complete cardiac cycle has to be imaged. This is equivalent to the annual recommended total radiation dose (see [27]). TEE is a newer technique, but highly invasive: a wire is inserted through the esophagus of the patient involving a cumbersome procedure. TEE is most of the time applied to monitor valve surgeries (see [19]), and is therefore rarely applied in a diagnostic phase.
Magnetic resonance imaging (MRI) allows to image anatomical structures in a non-invasive and non-ionising way. Unfortunately, the aortic valve geometry is not directly visible with MRI, since the valve thickness is smaller than the image voxel size.
However, visual inspection of 3D Flow MRI imaging (4D flow, see [18]) data sets clearly shows that the valvular shape could be retrieved from the flow pattern in the valve surroundings. This fact motivates to formulate, analyze and numerically assess the inverse problem of inferring rigid obstacles from interior velocity measurements, with the ultimate goal of recovering the shape of cardiac valves at the opening position from velocity images.
Available approaches studied and reported in the mathematical literature are limited to the detection of obstacles in viscous fluid flows using boundary stress measurements, which limits the inverse problem to oversimplified shapes, usually of circular nature [3,7,13,14]. Noteworthy, boundary stress measurements would mean in practice to introduce a catheter close to the valve.
The novelty of this work is the identification of flow obstacles by a distributed resistance inverse method. That is, we propose the incorporation of a large resistance term that allows to model the effect of an obstacle in the viscous fluid. This method reduces the obstacle identification to a simpler potential inverse problem. This idea is inspired by [5,12,20], who model cardiac valves using a resistive immerse surface (RIS) given by a Dirac distribution. The problem of using RIS for valve identification would be the need of modifying the discretization of the domain at every iteration of the identification procedure. In contrast, the distributed resistance term that we propose here allows us to work in a fixed domain to solve the valve shape identification problem. This distributed resistance method can also be useful to estimate porosity in porous media flows following Brinkman's law [6,11].
The paper is structured as follows. In section 2, a parameter identification problem is defined for the Oseen equations (as a linearization of Navier-Stokes equations) adding a resistive potential term with the form γu, where u is the velocity field and γ 0 is a function that takes large values in zones where is the obstacle should be located, or 0. In section 3, a proof of the existence of minimizers using appropriate techniques is presented. In sections 4 and 5, first and second order optimality conditions for some suitable cost functions are established, motivating the choice of suitable spaces for the parameter and the feasibility of implementing optimization algorithms to numerically solve this problem. Section 6 presents numerical examples that illustrate the feasibility of the proposed penalty method using Navier-Stokes equations from a 2D reference case representing the blood flow passing through the aortic valve. We reconstruct the position of the valve from global or local velocity field measurements. For this purpose, synthetic images based on MRI and ultrasound type measurements are used. Finally, conclusions and future work are presented in section 7.

Model problem
Consider a non-empty bounded domain Ω ⊆ R N , N ∈ {2, 3}. The Lebesgue measure of Ω is denoted by |Ω|, which extends to lesser dimension spaces. The norm and seminorms for Sobolev spaces W m,p (Ω) is denoted by · m,p,Ω and |·| m,p,Ω , respectively. For p = 2, the norm, seminorms and inner product of the space W m,2 (Ω) = H m (Ω) are denoted by · m,Ω , |·| m,Ω and (·, ·) m,Ω , respectively. Also, · ∞,Ω denotes the norm of L ∞ (Ω). The spaces H m (Ω) and The notation for norms, seminorms and inner products will be extended from W m,p (Ω) or H m (Ω). Consider and a ∈ W 1,∞ (Ω) such that div a = 0 and a · n 0 on Γ N , where n is the outer normal vector. The model problem is defined by The modified version of Oseen equations given by equation (2) models the velocity u and pressure p of a fluid that passes through the control volume Ω. The term γu, with γ 0, represents the distributed resistance that Ω imposes to the fluid. A zero value of γ represents that the fluid follows the Oseen equations. As the γ values increases in a certain area, the magnitude of u tends to decrease to 0. According to Brinkman's law [6,11], 1/γ is an indicator of permeability. The media is completely permeable when γ = 0, while γ tends to +∞ in zones where there are obstacles. Bounded values of γ model porous media, where porosity decreases as γ increases.
In order to simplifying the analysis of this problem, an equivalent formulation with homogeneous Dirichlet boundary condition is proposed.
, the model problem can be written in a new equivalent form.
where (x) − = min {x, 0}. For this case, the existence of a solution of the Oseen equations is verified [10].

Existence of solution
The demonstration of the existence of an optimal solution presented follows the same scheme as [4]. A first step is to introduce some helpful notations.
Γ D (Ω) and p ∈ L 2 0 (Ω) are the unique variational solutions of the problem (4). In what follows, it is denoted H = H 1 Γ D (Ω) × L 2 0 (Ω), which is a Banach space behind the norm

Remark 5.
It is possible to ensure that A is well defined, because (4) has an unique solution for every γ ∈ A (see lemma 5.8 on [23]). Futhermore, there exists a constant C > 0, independent of γ, ν, u D and G, such that The uniformly boundedness of v 1,Ω and p 0,Ω are obtained because γ ∞,Ω M.

Theorem 6. Problem (3) has at least one solution γ * ∈ A with optimal states
Proof. First, J (γ) = J (γ, A (γ)) is bounded below by 0 and Then, {γ n } n∈N is uniformly bounded in H s (Ω). Since A is weakly closed, there exists a subsequence (denoted by {γ n } n∈N ) such that γ n γ * in H s (Ω), with γ * ∈ A. In particular, γ n γ * in L 2 (Ω). By the same way, {v n } n∈N and {p n } n∈N are also uniformly bounded (see remark 5) on H 1 Γ D (Ω) and L 2 0 (Ω), respectively. Applying the Sobolev embedding theorem (see section 6.6 in [16]), there exists a subsequence (also denoted by {v n } n∈N ) such that v n v * weakly in 6) when N 3, or for p 2 when N = 2. In particular, v n → v * in L 4 (Ω) and L 2 (Ω) .
Repeating this argument, there exists a subsequence also denoted by {p n } n∈N such that p n p * in L 2 0 (Ω) .

First order necessary optimality condition
In order to implement a descent method to numerically solve this problem, it is necessary to establish the differentiability of functional J. However, this directly depends on the differentiability of map A. From this result, it is possible to establish necessary optimality conditions. 1 , p 1 can be described by the weak solution of the following problem

Theorem 7. The map A is Fréchet
,ω and C (γ) = α 2 γ 2 s,Ω , an expression for Frechét derivatives of B and C is given by Applying chain rule, it is obtained that In order to reduce this expression, the following definition is introduced similarly to [1].
where χ ω is the indicator function of ω.
Using this definition, it is possible to rewrite J (γ 1 ) depending of the adjoint state. That expression is simpler to analyze, since it depends on the adjoint state, allowing a simple form of a first order optimality condition using a variational inequality. Theorem 9. Let γ, γ 1 ∈ A and s 0. Then, are the states and adjoint states of γ * , respectively.

Proof. First, using integration by parts with the adjoint states [ϕ, ξ] as tests functions, it is obtained that
Then, Later, if γ * ∈ A is optimal for the problem, then proving this theorem.

Second order sufficient optimality condition
The stability results of the optimization algorithms depend, for the most part, on the existence of the second derivative of J. Likewise, it is possible to establish second order sufficient optimality conditions. In first place, it is necessary to introduce new forms of embedding inequalities.

Theorem 10.
There exists a constant q * > 2, dependent only of Ω, such that for each p ∈ [2, q * ] there exists a constant C > 0, dependent only of Ω, M, ν and p such that Then, there exists a constant C p > 0, depending only of p and Ω, such that Futhermore, using theorem 10, So, v 1,p,Ω is uniformly bounded for each p ∈ [2, q * ]. Also, theorem 1 in [21] allows us to be more precise about the asymptotic behavior of q * . Indeed, if the value of γ ∞,Ω increases, q * decreases to 2. Proof. See section 6.6 in [16] Remark 13. Using this embedding result with s N q * and p = 2, lemma 12 can be written by Then, for 2 p In order to establish and prove second order optimality conditions, a first step is to show that the map A is twice Fréchet-differentiable. and Proof. See appendix B Let γ, γ 1 , γ 2 ∈ A. An expression for the Fréchet second derivative of J (γ) on directions γ 1 and γ 2 is given by where, reasoning as in the proof of theorem 9, In consequence, In what follows, a second order optimality condition is proved. For this, a series of technical results are required. Consider r such that 1 2}. Using this, it is possible to obtain the following estimations.
Using theorem 10 twice and triangular inequality, there exist c 1 , Taking C = c 2 G 1,∞,Ω , the estimation is proved.
Proof. First, testing the equation of v h,k with v h,k and applying Friedichs-Poincaré and Hölder inequalities, there exists a constant c 1 > 0 such that , remark 13 points out that there exists c 3 > 0 such that γ k 0,r,Ω c 3 γ k s,Ω . In conclusion, proving this estimation.
Proof. Using Cauchy-Schwartz and Friedrichs-Poincaré inequalities, there exist c 1 , Using previous lemmas, Friedrichs-Poincaré inequality and remark 13, there exist c 3 , c 4 , c 5 , c 6 > 0 such that The proof of the second inequality is similar.
Proof. First, applying Hölder and Friedrichs-Poincaré inequalities, there exist c 1 , c 2 > 0 such that By repeating the previous procedures, since |ϕ h | 1,Ω is uniformly bounded, there exist c 2 , c 3 > 0 such that Proof. Applying triangular inequality, we obtain where every term were bounded with estimations of the form C h ∞,Ω |γ 1 | s,Ω |γ 2 | s,Ω (see lemmas 17 and 18). In conclusion, there exists L > 0 such that Corollary 20. There exists L > 0 such that, for every θ ∈ [0, 1] Finally, a second order sufficient optimality condition is presented and proved.
Theorem 21. Let s N q * and γ * ∈ A such that γ * verifies the first order optimality condition. If there exists δ > 0 such that Then, there exist σ, ε > 0, independent of γ and γ * , such that In consequence, J has a local minimum at γ * .
proving the theorem.

Remark 22.
The result obtained in theorem 21 is conditioned to s N q * . Considering the most pessimistic case, when q * = 2, it is possible to establish that the simplest spaces for γ * , where theorem 21 is valid, are H 1 (Ω) and H 2 (Ω) when N = 2 and N = 3, respectively.

Numerical results
In this section, the previous analysis is complemented by numerical experiments for the new parameter identification problem. Realistic synthetic cases are analyzed in 2D, which represent a longitudinal section of the structure of a cardiac valve in an arbitrary position. However, the numerical results presented below were obtained using the Navier-Stokes equation. So, the numerical problem to be solved is given by This special formulation of (9) is similar to (2) in terms of the resistance term. In the following subsections, the configurations of the reference case is explained, as well as the numerical  solutions of the inverse problems associated with MRI images or Doppler ultrasound. Due to our lack of real images, our experiments are fully synthetic (figures 1-3).

Reference test
First, a reference geometry is defined. This geometry Ω represents the area around the aortic valve. The inflow Γ I has a parabolic profile following Poiseuille's law given by  In the walls of the structure, represented by Γ W , the fluid has no-slip conditions given by u = 0. The valves are modeled on the resistance term γu using the function γ. This function assumes a constant value M 1 in the regions where the valve is and assumes the value 0 where the valve is not. In order to define γ, a parabolic figure is drawn on each side Ω with an approximate thickness of 1 mm. When the valves are fully closed, they are symmetrical with respect to the vertical axis of symmetry. When the valves are open, these parabolas are rotated with respect to a reference system whose origin is at the point where the valve coincides with the walls of the structure given by Γ W . Since blood flow is modeled, the kinematic viscosity is considered equal to ν = 0.035 cm 2 s −1 , the density is assumed as ρ = 1 g cm −3 , and the dimensions d = 2 cm and U = 30 cm s −1 , resulting in a peak Reynolds number on the inlet of The Navier-Stokes equations are discretized using the finite element method (FEM) with Taylor-Hood elements (P 2 for the velocity u and P 1 for the pressure) on an unstructured triangular mesh. The mesh used was generated by domain triangulation with h = 0.05 cm, which corresponds to 9950 elements and 4976 nodes. The solver is implemented using the finite element library FEniCS [2] with the default configuration. To solve the nonlinear problems, a Newton's method was used. The resistance γ is defined as a discontinuous function with discrete values in each element, assuming the value of M = 10 8 if the element intersects the valve or, otherwise, assumes the value 0. We define the set O, that represents the valve inside of Ω, by O = x ∈ Ω|γ(x) = 10 8 .

Numerical solution of the inverse problem
Using the solution of the reference test as reference solution u R , the following version of the model problem is solved numerically using FEniCS and Dolfin-adjoint For this example, we considered a measurement area ω = Ω and the values C = 10 3 , α = 10 −4 and β = 10 −8 . The use of two different weights for the norm and seminorm is consistent with the theoretical analysis of the previous sections, so this problem has a solution. The dolfin-adjoint library [28] allows to implement automatic derivation of the discrete adjoint equations for pde models and implement minimization algorithms from the Python 3 libraries. In particular, the L-BFGS-B algorithm (see section 4.3 in [17]) was used with the following stopping criteria on the step k where (∇J k ) j is the jth component of the projected gradient on the step k. To start the algorithm, γ 0 = 0 was used as the initial solution.
As a way to define a valve reconstruction algorithm sketch, we follow these steps. Numerical results are presented in figure 4. The polyline is drawn in white, which presents a great approximation to the interface between the different values of the reference given by γ.
The optimal γ * has values close to 0 in the areas before and after the valves. On the other hand, in the interior area where there are no valves, the optimal solution takes values close to 0. Likewise, the magnitude and direction of u * is similar to u R , where u * corresponds to the optimal state.  It is necessary to corroborate that this algorithm is able to solve the inverse problem measuring only a part of u R given by the reference velocity. In this case, choosing ω as the area where the valves should be (see figure 5 right), given by the expected result should be similar to that found above. Figure 5 presents the numerical results, the polyline and the references. In the case of the reference solution u R , a pink rectangle was drawn that allows delimiting the measurement area ω. In particular, this problem required more iterations than the case with measurements on Ω, obtaining very similar results.

Measurements of MRI type
Using the reference solutions, it is possible to generate synthetic measurements that represent the behavior of MRI velocity type measurements. A 2D plane is chosen on which the velocity is measured in one specified direction d ∈ R 2 , with |d| = 1. Then, the measurement is given by The MRI type velocity measurement data is represented in a mesh of uniform quadrilaterals of 2 mm × 2 mm. They are obtained by projecting u R to the Q 0 finite element space, given by piecewise constant discontinuous functions on the quadrilateral mesh. In practical applications, the input speed is unknown and must be estimated. Then, a direction d for this MRI can be chosen as that which is an inner normal to Γ I . To determine U, the projection of u R in L 2 (Γ I ) with respect to Ux(2 − x) is calculated by Thus, U approaches in a least squares sense. Figure 6 shows the synthetic MRI generated from the reference solution.
The new problem to solve is given by where the values C = 10 3 , α = 10 −4 and β = 10 −8 were used. It is possible to prove the existence of solution of this problem in the same way as in the proof of the theorem. Figure 7 shows the numerical results and the references. The parameter γ has some jumps that coincide with the vertical lines between the MRI voxels, while the polyline tends to have segments parallel to those vertical lines, but it acceptably approximates the space between the valves. The white noise intensity in the velocity measurements from MRI is proportional to the velocity encoding parameter (also called VENC [18]) of the scan. This parameter is configured with a value greater than the maximum expected velocity, in order to eliminate signal aliasing. Then, the noise in all voxels is proportional to the maximum velocity in the measurement area. In the clinical practice it can be expected that high-quality MRI contains a velocity noise of 10% of the maximum velocity in each voxel [18]. Gaussian noises were added to this MRI with a standard deviation of 5%, 10% and 20% of the maximum value of u R . Figure 8 shows the results of this experiment with a 5% of Gaussian noise, but changing the weights to α = 10 −4 and β = 10 −6 in terms to decrease the effects of noise. The results are similar to the experiment without noise in terms to the tendency of the polyline to approximate the valve shape and draw lines parallel to the voxels.
This approximation seems weaker as noise increases, in the sense that the polyline has a lower quality in its approximation and that the gamma function tends to overfit the data. The following two figures show the result of the experiment with 10% and 20% of Gaussian noise, respectively. The noise is exactly the same than the 5% Gaussian noise case, but increasing the level noise. Table 1 shows the mean square error (MSE) between the reconstructed valve given by the polylines obtained using MRI in figures 7-10, and the polyline obtained in the reference test (see figure 4). To quantify this error, we consider only the points of the polyline in figure 4 that are at a distance less than or equal to 0.5 mm from O. There are minor differences between the valve reconstructions for the cases with a noise level of 5% and 10%. However, the quality of the reconstruction decreases when the level noise increases up to 20%. Figure 11 shows the reconstructed polyline for three experiments with independent noises at same level of noise. Each polyline color represents the final result of this experiment with Figure 8. Real γ and reconstructed valve, optimal γ, optimal u, and synthetic MRI (from left to right). MRI with 5% of noise. 1.5016 × 10 −2 Figure 9. Real γ and reconstructed valve, optimal γ, optimal u, and synthetic MRI (from left to right). MRI with 10% of noise.
a Gaussian noise independent of the others, whose amplitude was adjusted for the noise level of 5%, 10% and 20%. We can see that the rebuilt valves are more irregular as the noise level increases, especially for 20% of level noise, which is consistent with table 1.

Measurements of ultrasound imaging type
It is possible to generate synthetic measurements similar to ultrasound images. From a given focal point, directional velocities are measured at each point in a domain with the form of Figure 10. Real γ and reconstructed valve, optimal γ, optimal u, and synthetic MRI (from left to right). MRI with 20% of noise.   The measurement data is represented in a structured triangular mesh of 2116 nodes and 4050 triangles with h = 0.067 cm. The measurements are obtained by interpolation of u R to the P 1 finite element space defined on the structured mesh ( figure 12). Again, the input speed is unknown and must be estimated from a preliminary measurement. To determine U, the projection of u R in L 2 (Γ I ) with respect to Ux(2 − x) [n · d (x)] is calculated by The new problem to solve is given by subject to −ν u + (∇u) u + ∇p where ω is the subdomain on Ω covered by the ultrasound imaging and the values C = 10 3 , α = 10 −4 and β = 10 −8 were used. This problem has solution, and the proof of the existence of solution is similar to the proof of the theorem. Figure 13 presents the numerical results, the polyline and the references. The results obtained are very similar to those obtained in the reference tests with measurements in the entire domain and in a subdomain.
There are multiple sources of noise in ultrasound images, such as reverberation, ghosting, or fake echo, among others. Therefore, it is not possible to define a simplified structure of noise in this type of images [9,25]. In order to study the effect of the noise, we assume an additive Figure 14. Real γ and reconstructed valve, optimal γ, optimal u, and synthetic ultrasound imaging (from left to right) with 10% of noise. Gaussian noise that is proportional to the maximum of |u R | in the measurement area. From [24], an experiment with a straight-vessel phantom had a 9.6% of standard deviation. Then, is acceptable to simulate an ultrasound imaging with a velocity noise of 10% of the maximum velocity in each node. Gaussian noises were added to this ultrasound imaging with a standard deviation of 10% and 20% of the maximum of |u R |. Figures 14 and 15 show the results of this experiment with a 10% and 20% of Gaussian noise, respectively, but changing the weights to α = 10 −4 and β = 10 −6 in terms to decrease the effects of noise. The results are similar to the experiment without noise in terms to the tendency of the polyline to approximate the valve shape. Table 2 shows the MSE between the reconstructed valve given by the polylines obtained using MRI in figures 13-15, and the polyline obtained in the reference test with measurement in a subdomain (see figure 5). The MSE was quantified by the same way as in table 1. There are minor differences between the valve reconstructions for the cases with a noise level of 10% and 20%.

Conclusions
We have presented a new parameter identification problem for the Oseen equation, contributing to the detection of obstacles in fluid dynamic studies. For this problem, the existence of solution and optimality conditions were determined. The study of optimality conditions is necessary to validate the use of optimization algorithms. In particular, the second order optimality condition allowed to specify the H s space where the parameter γ belongs.
The numerical experiments developed in this work allow us to conjecture that it would be possible to extend this analysis for the Navier-Stokes equation. Indeed, the numerical tests without noise had satisfactory results in terms of rebuilding the simulated valve. Furthermore, experiments with Gaussian noise show that disturbances in the reference image are reflected in bounded changes in the parameter and in the reconstructed valve. The quality of the solutions is worse as noise increases, as expected.
There are several ways to deepen this work. From an analytical perspective, it is necessary to repeat the theoretical analysis for problem (9), in the same way as we worked with problem (1), looking for obtaining similar results. Future work will tackle the uniqueness and stability of the solution of this new inverse problem. Related works on obstacle estimation via shape optimization in the Brinkman problem [29] or source estimation in the Brinkman [26] and transient Navier-Stokes equations [15] could serve as a starting point.
Based on the favorable numerical results, the next step is to perform experiments with 3D domains, including structures and real images (both MRI and ultrasound imaging), whether obtained from phantoms or real patients. The 3D case is of medical interest, since it will contribute to simplify the detection of defects in the function of aortic valves, among other structures. Therefore, designing an algorithm to recreate the aortic valve in 3D is also one of the future improvements.