Download the proposal

Type of internship

Master 2 internship/ engineering school end of study internship

Duration

6-month internship between January and September 2022

Context

Photoacoustic tomography (PAT) is a recent multi-wave image modality that shows great potential in biological and clinical research. This modality combines an optical excitation (laser) of the tissue with an ultrasound detection allowing to generate high-resolution images of optical absorption at large depths. The optical absorption in biological tissues can be due to endogenous molecules such as hemoglobin or melanin, and to injected contrast agents. This absorption results in a local increase of the temperature in the tissue, creating a small pressure increase which in turn generates an acoustic wave. In 3D photoacoustic tomography (PAT), the ultrasound wavefield is measured, as an electrical signal, by several sensors distributed on a surface surrounding the tissue. The reconstruction of the image is then performed through the resolution of an inverse problem from the observed (measured) signals (see Fig.1). The PAT modality is therefore intricate with its numerical reconstruction method and as such is an instance of computational imaging systems that have transformed modern imaging.

Physical modelling
Physically, the PAT can be modelled with two different phenomena:

The wave propagation in the tissue The initial pressure p₀ : R³ → R generated by the optical absorption and to be recovered in the image is mapped to the whole pressure field at all times p : R³ × [0,T] →R by a linear operator H : p₀ → p. This linear operator represents the wave equation: a second order linear partial differential equation (PDE). Once discretized, the matrix H is of size (N³× n_T) × N³ where N³ is the size of the grid discretizing the tissue domain and n_T is the number of time samples. In practice this matrix is huge and cannot be stored.
The measurements In PAT the pressure field p is measured by n sensors outside the tissue. Each sensor is associated to a 2D surface S_i ⊂ R³ that records the pressure. The associated measurement operator M_i maps the pressure p to the signal. The global measurement operator is thus the concatenation of all measurements.

Overall, the physical modeling leads to the following linear inverse problem:

g = Ap₀                                                    (1)

with A = MH and where g are the observed signals and p0 the initial pressure to be recovered. Once discretized, the matrix A is of size (nn_T) × N³. Note that the optimization algorithms implemented to solve (1) require the iterative application of A and its adjoint A∗. Great care must be taken when implementing these matrix-vector products as they are critical for the efficiency of the reconstruction.

Numerical issues
Two issues make the resolution of this inverse problem difficult. First, the computation of H through wave propagation simulations is numerically intensive (3D + time). Second, the adjoint A∗ is often implemented through back-propagation which has a simple expression when the sensors are isotropic and reduced to points. This however, requires to finely discretize the surfaces of each sensor, say with m points (that has to be chosen such that m ∝ N² points), leading to a matrix A with m more nonzero coefficients. This prevents the reconstruction of the image p0 in decent times leading researchers to resort to approximate sensor geometries which degrades the quality of reconstructed images.

References
[1] Rosenthal, Ntziachristos, Razansky, Acoustic Inversion in Optoacoustic Tomography: A Review. Curr. Med. Imaging Rev. 9, 318 (2013)
[2] X. L. Dean-Ben, V. Ntziachristos and D. Razansky, ”Acceleration of Optoacoustic Model-Based Reconstruction Using Angular Image Discretization,” in IEEE Transactions on Medical Imaging, vol. 31, no. 5, pp. 1154-1162 (2012)
[3] E. Candes, L. Demanet, L. Ying, Fast Computation of Fourier Integral Operators, SIAM Journal on` Scientific Computing, 29:6, 2464-2493 (2007)
[4] P. L. Combettes, J.-C. Pesquet, Proximal Splitting Methods in Signal Processing, in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer New York, 185-212 (2011)
[5] N. Komodakis and J. Pesquet, ”Playing with Duality: An overview of recent primal-dual approaches for solving large-scale optimization problems,” in IEEE Signal Processing Magazine, vol.32, no. 6, pp. 31-54 (2015)

Objectives and tasks

In preliminary works, we investigated an original approach, appearing to be new in the PAT community, that allows to fast reconstruct p0 without compromise on its quality. This approach is based on a fine implementation of the matrix-vector products with A and A∗ in the optimization algorithms. With this approach, one iteration of the algorithm now costs O(nN3) operations.

The objectives of the internship are to make this proof of concept a universal method that can be used routinely in laboratories. We have, for instance, identified several needs:

• Although sensors are fully characterized in our method, they have an electrical impulse response that is crucial to model and incorporate in our method.
• The current implementation is made for images in 2D albeit real images are 3D. In order to scale the code in this setting, we need to implement the algorithms on massively parallel architectures such as graphics processing units (GPU). A speed-up of a factor ×60 can be expected.
• The current implementation is based on fast Fourier transforms (FFT) that are extremely well implemented on any scientific computing library making it universal. However, the Fourier transform implicitly assumes periodic boundary conditions which forced us to significantly extend the domain to avoid the periodic waves to pollute the signals. A possible solution would be to implement a Perfectly Matched Layer which might lead to significant speed-ups (‘×30) in 3D.
• We believe that the matrix vector-products in our method implemented in O(nN³) which depends linearly in the number of sensors can be further reduce to O(rN³) with r << n using dimensionality reduction techniques. This would have a dramatic impact on the speed of reconstruction algorithms. A speed-up of a factor ×10 can be expected.
• If time allows, we might also investigate the co-conception of the PAT device. This means to adequately design the sensor locations to optimize the quality of the reconstructed images. This would allow to reduce the number of measurements needed and speed up the acquisition of the signals which is important in research and clinical applications.

Overall, the previous ideas might lead to a speed-up of order 104 which would make the algorithms efficient enough to reconstruct high-quality images without any compromise on their quality. This internship could have a great impact on the PAT community.

The candidate will be trained and could develop skills in optimization, image processing, high performance computing and approximation theory. These competences are actively being in demand in the industry and the academic research.

Profile

Candidate profile: master of science with strong skills in statistics, signal processing, machine learning, or optimization. Languages: Python/Matlab.
The internship might lead to a PhD thesis.

Internship Stipend

about 600 euros per month

Contact

Paul Escande (paul.escande[at]univ-amu.fr), Caroline Chaux (caroline.chaux[at]univ-amu.fr).
Jérôme Gateau (jerome.gateau[at]sorbonne-universite.fr)

Application: send resume, master grades and motivation letter to supervisors; please name documents as LASTNAME Firstname-Resume.*, LASTNAME Firstname-Grades.* or LASTNAME Firstname-Letter.*

Location: Signal and Image team I2M, Institut de Mathématiques de Marseille.