Eikonal Equation Research Articles

Feasibility study of clinical target volume definition for soft-tissue sarcoma using muscle fiber orientations derived from diffusion tensor imaging

Objective. Soft-tissue sarcoma spreads preferentially along muscle fibers. We explore the utility of deriving muscle fiber orientations from diffusion tensor MRI (DT-MRI) for defining the boundary of the clinical target volume (CTV) in muscle tissue. Approach. We recruited eight healthy volunteers to acquire MR images of the left and right thigh. The imaging session consisted of (a) two MRI spin-echo-based scans, T1- and T2-weighted; (b) a diffusion weighted (DW) spin-echo-based scan using an echo planar acquisition with fat suppression. The thigh muscles were auto-segmented using the convolutional neural network. DT-MRI data were used as a geometry encoding input to solve the anisotropic Eikonal equation with the Hamiltonian Fast-Marching method. The isosurfaces of the solution modeled the CTV boundary. Main results. The auto-segmented muscles of the thigh agreed with manually delineated with the Dice score ranging from 0.8 to 0.94 for different muscles. To validate our method of deriving muscle fiber orientations, we compared anisotropy of the isosurfaces across muscles with different anatomical orientations within a thigh, between muscles in the left and right thighs of each subject, and between different subjects. The fiber orientations were identified reproducibly across all comparisons. We identified two controlling parameters, the distance from the gross tumor volume to the isosurface and the eigenvalues ratio, to tailor the proposed CTV to the satisfaction of the clinician. Significance. Our feasibility study with healthy volunteers shows the promise of using muscle fiber orientations derived from DW MRI data for automated generation of anisotropic CTV boundary in soft tissue sarcoma. Our contribution is significant as it serves as a proof of principle for combining DT-MRI information with tumor spread modeling, in contrast to using moderately informative 2D CT planes for the CTV delineation. Such improvements will positively impact the cancer centers with a small volume of sarcoma patients.

Boundary Estimation from Point Clouds: Algorithms, Guarantees and Applications

We investigate identifying the boundary of a domain from sample points in the domain. We introduce new estimators for the normal vector to the boundary, distance of a point to the boundary, and a test for whether a point lies within a boundary strip. The estimators can be efficiently computed and are more accurate than the ones present in the literature. We provide rigorous error estimates for the estimators. Furthermore we use the detected boundary points to solve boundary-value problems for PDE on point clouds. We prove error estimates for the Laplace and eikonal equations on point clouds. Finally we provide a range of numerical experiments illustrating the performance of our boundary estimators, applications to PDE on point clouds, and tests on image data sets.

Simulation of non-transmural ablation lines that effectively block electrical signal propagation in the heart

Atrial fibrillation (AF), the most common cardiac rhythm disorder, increases patient mortality and social burdens. To achieve a permanent cure for AF, catheter-based radiofrequency ablation destroys the ability of local cardiac tissue to conduct electrical signals, thus terminating the arrhythmia. It is generally assumed that only lesions that penetrate the whole thickness of the heart wall can block the electrical signals. However, there was no in-depth theoretical research on the relationship between ablation depths and electrical signal blocking. We use the Barkley model and the Luo-Rudy model to simulate the spiral activity during AF episodes. We find that, when the ablation depth exceeds a certain critical value, although there is still a tiny slit, the electrical signals cannot pass through, and this phenomenon can be explained by the eikonal equation. Moreover, the slit has a filtering effect. When the height of the slit is within a certain interval, the signal period detected on the other side of the ablation line is almost twice that of the wave source, which is caused by the slow conduction at the slit. In addition, we find that when the ablation line is not deep enough to completely block the passage of the spiral wave, increasing the width of the ablation line can achieve blocking. We simulate two different obstacles caused by the ablation under the no-flux boundary condition constructed by the phase-field method and under the condition of fixed potential, and get different results. The electrotonic leak current may be the main reason for this difference. These results can help clinicians understand the blocking phenomenon in ablation procedures and develop more effective ablation strategies.

Geometry of the Singular Set of Hypersurfaces and the Eikonal Equation

Geometry of the Singular Set of Hypersurfaces and the Eikonal Equation

A Continuous Finite Element Method with Homotopy Vanishing Viscosity for Solving the Static Eikonal Equation

A Continuous Finite Element Method with Homotopy Vanishing Viscosity for Solving the Static Eikonal Equation

Erratum to “Improving the accuracy of the expanded anisotropic eikonal equation at larger offsets using Levin’s T-transformation”

Erratum to “Improving the accuracy of the expanded anisotropic eikonal equation at larger offsets using Levin’s T-transformation”

Surface wave tomography using active sources in engineering applications: A 3D experiment at a test site of known conditions

In this study a specific workflow is proposed for the application of Surface wave tomography in near surface applications using active sources. The proposed methodology can be implemented in 2D or 3D acquisition schemes, resulting in high resolution phase and group surface wave velocity maps that can be used for Vs estimation in the survey area. Our methodology starts with the implementation of a pre-processing procedure in which a series of common steps for the quality control and optimization of the data is applied. Then, a reliable reference phase velocity dispersion curve can be estimated internally from our data. This reference curve is used for the automatic assignment of phase and group delays between nearby stations in the time domain. Afterwards, high precision frequency-dependent phase and group velocity maps can be produced through inversion with the Eikonal equation. The proposed methodology is evaluated by application in a 3D acquisition example in a controlled test of known conditions. The resulting frequency-dependent phase and group velocity maps delineate efficiently the expected lateral and depth variations in the survey area. The comparison of the results with available geotechnical information and other geophysical methods confirms the method's validity, demonstrating the effectiveness of the technique and its value in near surface applications as an attractive tool for the calculation of quasi-3D shear wave velocity for a survey area.

Neighbor-Gradient Single-Pass Method for Solving Anisotropic Eikonal Equation

We develop a single pass method for approximating the solution to an anisotropic eikonal equation related to the anisotropic min-time optimal trajectory problem. Ordered Upwind Method (OUM) solves this equation, which is a single-pass method with an asymptotic complexity. OUM uses the search along the accepted front (SAAF) to update the value at considered nodes. Our technique, which we refer to as "Neighbor-Gradient Single-Pass Method", uses the minimizer of the Hamiltonian, in which the gradient is substituted with neighbor gradient information, to avoid SAAF. Our technique is considered in the context of control-theoretic problem. We begin by discussing SAAF of OUM. We then prove some properties of the value function and its gradient, which provide the key motivation for constructing our method. Based on these discussions, we present a new single-pass method, which is fast since it does not require SAAF. We test this method with several anisotropic eikonal equations to observe that it works well while significantly reducing the computational cost.

Geodesics in Jet Space

The space $J^k$ of $k$-jets of a real function of one real variable $x$ admits the structure of Carnot group type. As such, $J^k$ admits a submetry (\sR submersion) onto the Euclidean plane. Horizontal lifts of Euclidean lines (which are the left-translates of horizontal one-parameter subgroups) are thus globally minimizing geodesics on $J^k$. All $J^k$-geodesics, minimizing or not, are constructed from degree $k$ polynomials in $x$ according to Anzaldo-Meneses and Monroy-Per\'ez, reviewed here. The constant polynomials correspond to the horizontal lifts of lines. Which other polynomials yield globally minimizers and what do these minimizers look like? We give a partial answer. Our methods include constructing an intermediate three-dimensional "magnetic" sub-Riemannian space lying between the jet space and the plane, solving a Hamilton-Jacobi (eikonal) equations on this space, and analyzing period asymptotics associated to period degenerations arising from two-parameter families of these polynomials. Along the way, we conjecture the independence of the cut time of any geodesic on jet space from the starting location on that geodesic.

Explicit complex-valued solutions of the 2D eikonal equation

We present a method to obtain explicit solutions of the complex eikonal equation in the plane. This equation arises in the approximation of the Helmholtz equation by the WKBJ or Evanescent Wave Tracking methods. We obtain the complex-valued solutions (called eikonals) as parametrizations in a complex variable. We consider both the cases of constant and non-constant index of refraction. In both cases, the relevant parametrizations depend on some holomorphic function. In the case of a non-constant index of refraction, the parametrization also depends on some extra exponential complex-valued function and on a quasi-conformal homeomorphism. This is due to the use of the theory of pseudo-analytic functions and the related similarity principle. The parametrizations give information about the formation of caustics and the light and shadow regions for the relevant eikonals.

Solving the eikonal equation for compressional and shear waves in anisotropic media using peridynamic differential operator

SUMMARY Traveltimes of compressional (P) and shear (S) waves have proven essential in many earthquake and exploration seismology applications. An accurate and efficient traveltime computation for P and S waves is crucial for the success of these applications. However, solving the eikonal equation with a complex phase velocity field in anisotropic media is challenging. The eikonal equation is a first-order nonlinear hyperbolic partial differential equation. It represents the high-frequency asymptotic approximation of the wave equation. The fast marching and sweeping methods are commonly used due to their efficiency in numerically solving the eikonal equation. However, these methods suffer from numerical inaccuracy in anisotropic media with sharp heterogeneity, irregular surface topography and complex phase velocity fields. This study presents a new method for the solution of the eikonal equation by employing the peridynamic differential operator (PDDO). The PDDO provides the non-local form of the eikonal equation by introducing an internal length parameter (horizon) and a weight function with directional non-locality. The operator is immune to discontinuities in the form of sharp changes in field or model variables and invokes the direction of traveltime in a consistent manner. The weight function controls the strength of association among points within the horizon. Solutions are constructed in a consistent manner without upwind assumptions through simple discretization. The robustness of this approach is established by considering different types of eikonal equations on complex velocity models in anisotropic media. The examples demonstrate its unconditional numerical stability and results compare well with the reference solutions.

Eikonal surface wave tomography with smoothing splines—application to Southern California

SUMMARY The densification of both permanent and temporary seismic networks has raised new interest in surface wave eikonal tomography from which phase velocity maps can be obtained without resolving a tomographic inverse problem. However, eikonal tomography requires to reconstruct traveltime surfaces from a discrete number of measurements obtained at the station locations, which can be challenging. We present a new method to reconstruct these traveltime surfaces with smoothing splines discretized in a regular 2-D Cartesian grid. We impose Neumann boundary conditions so that the phase gradients on the edges of the grid are equal to the apparent slownesses of the average plane wave along the normal direction measured by beamforming. Using the eikonal equation, phase velocity maps are then derived from the norm of the gradient of the interpolated traveltime maps. The method is applied to Rayleigh waves recorded by the Southern California Seismic Network to derive phase velocity surfaces. Robust, stable and finely resolved phase velocity maps at 25 and 33 s period are obtained after averaging the phase velocity maps derived from the analysis of a selection of recent large (Mw ≥ 6.5) teleseismic events. The phase velocity map at 25 s mainly constrains the thickness of the Southern California crust, with results that are in excellent agreement with previous tomographic studies.

On the solution of hyperbolic equations using the peridynamic differential operator

Numerical solution of hyperbolic differential equations, such as the advection equation, poses challenges. Classically, this issue has been addressed by using a scheme known as the upwind scheme. It simply invokes more points from the upwind side of the flow stream when calculating derivatives. This study presents a generalized upwind scheme, referred to as directional nonlocality, for the numerical solution of linear and nonlinear hyperbolic Partial Differential Equations (PDEs) using the peridynamic differential operator (PDDO). The PDDO provides the nonlocal form of the differential equations by introducing an internal length parameter (horizon) and a weight function. The weight function controls the degree of interaction among the points within the horizon. A modification to the weight function, i.e., upwinded-weight function, accounts for directional nonlocality along which information travels. This modification results in a stable PDDO discretization of hyperbolic PDEs. Solutions are constructed in a consistent manner without special treatments through simple discretization. The capability of this approach is demonstrated by considering time dependent linear and nonlinear hyperbolic equations as well as the time invariant Eikonal equation. Numerical stability is ensured for the linear advection equation and the PD solutions compare well with the analytical/reference solutions.

Numerical Solutions to Wave Propagation and Heat Transfer Non-Linear PDEs by Using a Meshless Method

The applications of the Eikonal and stationary heat transfer equations in broad fields of science and engineering are the motivation to present an implementation, not only valid for structured domains but also for completely irregular domains, of the meshless Generalized Finite Difference Method (GFDM). In this paper, the fully non-linear Eikonal equation and the stationary heat transfer equation with variable thermal conductivity and source term are solved in 2D. The explicit formulae for derivatives are developed and applied to the equations in order to obtain the numerical schemes to be used. Moreover, the numerical values that approximate the functions for the considered domain are obtained. Numerous examples for both equations on irregular 2D domains are exposed to underline the effectiveness and practicality of the method.

Upwind, No More: Flexible Traveltime Solutions Using Physics-Informed Neural Networks

The eikonal equation plays an important role across multidisciplinary branches of science and engineering. In geophysics, the eikonal equation, and its characteristics, are used in addressing two fundamental questions pertaining to seismic waves: what paths do the seismic waves take (its spreading)? and how long do they take? There have been numerous attempts to solve the eikonal equation, which can be broadly categorized as finite-difference and physics informed neural network (PINN) based approaches. While the former has been developed and optimized over the years, it still inherits some numerical inaccuracies and also the cost scales exponentially with the velocity model size. More importantly, it requires upwind calculations to satisfy the viscosity solution. PINNs, on the other hand, have shown great promise due to several features allowing for higher accuracy and scalability than conventional approaches. In this paper, we demonstrate another unique feature of PINN solutions, specifically its flexibility resulting from the global nature of its NN functional optimization, allowing for functional gradients referred to as automatic differentiation. This feature allows us to overcome the inability of conventional methods to handle large areas of missing information (gap) in the velocity model. We find empirically that the PINNs interpolation-extrapolation inherent capability enables us to circumvent a scenario when traveltime modelling is performed on velocity models containing gaps. Such a capability is crucial when performing traveltime modelling using the global tomographic Earth velocity model.

Consistency analysis of bilevel data-driven learning in inverse problems

One fundamental problem when solving inverse problems is how to find regularization parameters. This article considers solving this problem using data-driven bilevel optimization, i.e. we consider the adaptive learning of the regularization parameter from data by means of optimization. This approach can be interpreted as solving an empirical risk minimization problem, and we analyze its performance in the large data sample size limit for general nonlinear problems. We demonstrate how to implement our framework on linear inverse problems, where we can further show the inverse accuracy does not depend on the ambient space dimension. To reduce the associated computational cost, online numerical schemes are derived using the stochastic gradient descent method. We prove convergence of these numerical schemes under suitable assumptions on the forward problem. Numerical experiments are presented illustrating the theoretical results and demonstrating the applicability and efficiency of the proposed approaches for various linear and nonlinear inverse problems, including Darcy flow, the eikonal equation, and an image denoising example.

A linear finite-difference scheme for approximating randers distances on cartesian grids

Randers distances are an asymmetric generalization of Riemannian distances, and arise in optimal control problems subject to a drift term, among other applications. We show that Randers eikonal equation can be approximated by a logarithmic transformation of an anisotropic second order linear equation, generalizing Varadhan’s formula for Riemannian manifolds. Based on this observation, we establish the convergence of a numerical method for computing Randers distances, from point sources or from a domain’s boundary, on Cartesian grids of dimension 2 and 3, which is consistent at order 2/3, and uses tools from low-dimensional algorithmic geometry for best efficiency. We also propose a numerical method for optimal transport problems whose cost is a Randers distance, exploiting the linear structure of our discretization and generalizing previous works in the Riemannian case. Numerical experiments illustrate our results.

Learning cardiac activation maps from 12-lead ECG with multi-fidelity Bayesian optimization on manifolds

We propose a method for identifying an ectopic activation in the heart non-invasively. Ectopic activity in the heart can trigger deadly arrhythmias. The localization of the ectopic foci or earliest activation sites (EASs) is therefore a critical information for cardiologists in deciding the optimal treatment. In this work, we formulate the identification problem as a global optimization problem, by minimizing the mismatch between the ECG predicted by a cardiac model, when paced at a given EAS, and the observed ECG during the ectopic activity. Our cardiac model amounts at solving an anisotropic eikonal equation for cardiac activation and the forward bidomain model in the torso with the lead field approach for computing the ECG. We build a Gaussian process surrogate model of the loss function on the heart surface to perform Bayesian optimization. In this procedure, we iteratively evaluate the loss function following the lower confidence bound criterion, which combines exploring the surface with exploitation of the minimum region. We also extend this framework to incorporate multiple levels of fidelity of the model. We show that our procedure converges to the minimum only after 11.7 ± 10.4 iterations (20 independent runs) for the single-fidelity case and 3.5 ± 1.7 iterations for the multi-fidelity case. We envision that this tool could be applied in real time in a clinical setting to identify potentially dangerous EASs.

Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique

First-arrival traveltime tomography has been widely used for upper crustal velocity modeling, but it usually suffers from the problem of complex surface topography. To overcome this problem, we have developed a new topography-dependent eikonal tomography scheme that combines a developed accurate and efficient traveltime modeling method and introduces a flexible and robust adjoint inversion scheme in the presence of irregular topography. A surface-flattening scheme is used to handle the irregular surface, where the real model is discretized by curvilinear grids and the irregular free surface is mathematically flattened through the transformation from Cartesian to curvilinear coordinates. Based on this parameterization, the forward traveltime modeling is conducted by a monotone fast-sweeping method that discretizes the factored topography-dependent eikonal equation with a point-source condition. This algorithm can circumvent the source-singularity problem and decrease the numerical error in the vicinity of a point source in the curvilinear system. Then, the gradient-based inversion is used to minimize the misfit function, which is achieved by a matrix-free adjoint-state method without cumbersome ray tracing and explicit estimation of the Fréchet derivative matrix in the curvilinear coordinate system. The new tomographic scheme is evaluated through numerical examples with different seismic structures with complex topography, and then applied to a wide-angle profile acquired in the northeastern Tibetan Plateau. The results validate the effectiveness and efficiency of our tomography scheme in constructing shallow crustal velocity models with irregular topography.

Non-recursive transport of intensity phase retrieval with the transport of phase.

The transport of intensity equation (TIE) is a non-interferometric phase retrieval method that originates from the imaginary part of the Helmholtz equation and is equivalent to the law of conservation of energy. From the real part of the Helmholtz equation, the transport of phase equation (TPE), which represents the Eikonal equation in the presence of diffraction, can be derived. The amplitude and phase for an arbitrary optical field should satisfy these coupled equations simultaneously during propagation. In this work, the coupling between the TIE and TPE is exploited to improve the phase retrieval solutions from the TIE. Specifically, a non-recursive fast Fourier transform (FFT)-based phase retrieval method using both the TIE and TPE is demonstrated. Based on the FFT-based TIE solution, a correction factor calculated by the TPE is introduced to improve the phase retrieval results.