Documentation	Build Status

PhotoAcoustic operators (forward and adjoint) with a high level linear algebra interface. This package is an extension of JUDI for photoacoustic simulations and introduces the use of initial value problems.

This package is based on JUDI so please read that documentation to understand how to setup basic simulation parameters such as model size and receiver locations. Once a simulation has been designed in pure JUDI, the photoacoustic simulation is defined by using the judiPhoto operator and a judiInitialState as a source. Given F, a judiModeling operator describing the acoustic simulation and recGeometry a Geometry object describing the receiver setup we can setup and run a photoacoustic simulation with a point source as follows:

# Forward Photoacoustic operator
A = judiPhoto(F, recGeometry;)

# Photoacoustic source (n contains spatial dimensions)
init_dist = zeros(Float32, n)
init_dist[div(n[1],2), div(n[2],2)] = 1
p = judiInitialState(init_dist);

# Forward simulation
d_sim = A*p

A complete runnable script with the above example is here

In order to solve photoacoustic inverse problems in a variational framework, we also need the adjoint photoacoustic operator A^{\top}. In this package, we derive and implement the adjoint. The operator can be accessed with simple linear algebra notation:

# Forward simulation
d_sim = A*p

# Adjoint simulation (adjoint(A) also works)
p_adj = A'*d_sim 

The adjoint operator can be used to find sensitivities of data with respect to the initial source. This allows performing gradient descent and also iterative methods to solve a least squares variational framework. Furthermore, we also derive the adjoint simulation that defines the sensitivity with respect to the speed of sound model used in the underlying acoustic simulation. This enables independent or joint optimization of the initial value p and also acoustic parameters m. The sensitivity with respect to the acoustic model is related to the linearization of wave equation around a point model0 so we can implement this in the familiar JUDI way and apply the adjoint Jacobian on the data residual dD:

# Forward simulation
F0 = judiPhoto(model0, recGeometry;)
J = judiJacobian(F0, p)
dm = J'*dD

A promising new topic is scientific machine learning that marries together physics based operators with learned operators. We exemplify this with the deep image prior a technique that requires chaining together gradients of the forward operator (in our case related to the solution of a partial differential equation) and gradients of learned layers.

loss, grad = Flux.withgradient(Flux.params(unet)) do
    norm(A*unet(z) - y).^2
end
Flux.Optimise.update!(opt, ps, grad)

This package was written by Mathias Louboutin and Rafael Orozco from the Seismic Laboratory for Imaging and Modeling (SLIM) at the Georgia Institute of Technology.