Automatic skull-template alignment without a guidance image.

Oscar Bates
Dept. Earth Science and Engineering
Imperial College London
o.bates18@imperial.ac.uk
Carlos Cueto
Dept. Earth Science and Engineering
Imperial College London
Ciaran Coleman
Sonalis Imaging Limited
Cameron A. B. Smith
Dept. Bioengineering
Imperial College London
Lluis Guasch
Dept. Earth Science and Engineering
Imperial College London
Oscar Calderon Agudo
Dept. Earth Science and Engineering
Imperial College London


Abstract

Transcranial ultrasound must overcome the significant challenge of the human skull, limiting both imaging and therapeutic applications. While high-fidelity numerical simulations can compensate for skull-induced distortions, they require precise skull templates (typically derived from Computed Tomography, CT) and spatial alignment to the patient’s anatomy. Current alignment relies on concurrent Magnetic Resonance Imaging (MRI) for registration, introducing financial, logistical, and throughput barriers. To overcome these challenges, we present Manifold Optimisation for Full-Waveform Inversion (MOFI), a method to register skull templates without using a guidance image. Our method aligns the skull template by minimising the difference between simulated and observed radio-frequency acoustic data. We demonstrate that MOFI accurately recovers the position of skull templates in silico and in vitro, offering an alternative to MRI guidance-based registration. These results indicate that MOFI has the potential to be a practical alternative to MRI-guided approaches, reducing the barriers to wider clinical adoption of transcranial ultrasound.

2

Transcranial ultrasound encompasses numerous therapeutic and imaging applications. Currently, the major neuroimaging modalities, Magnetic Resonance Imaging (MRI) and X-ray Computed Tomography (CT), are large, expensive, immobile, high-power instruments that are deployable only in specialized facilities. Pulse-Echo Ultrasound is a safe, fast, and portable but transcranial B-mode imaging struggles with limited resolution [1]. Therapeutic ultrasound has few direct analogues in clinical use and encompasses neuromodulation [2], drug delivery [3], and ablation [4].

Ultrasound neurotechnologies are hindered by the high impedance and heterogeneous structure of the skull, which abberates, scatters, and attenuates the ultrasound wavefield [5]. These effects substantially degrade both therapeutic targeting [6] and transcranial imaging [7].

Recent work combining wide-aperture transcranial Ultrasound Tomography (t-UST) systems with Full-Waveform Inversion (FWI) has shown promise for recovering intracranial structure [7][9]. FWI formulates image reconstruction as a non-linear optimisation problem, estimating the acoustic parameter field that best matches measured radio-frequency acoustic data. However, the large impedance contrast between bone and soft tissue continues to pose severe challenges for conventional gradient-based optimisation. FWI is high dimensional, often involving tens to hundreds of millions of elements and is prone convergence issues—such as cycle skipping—that impair the reconstruction [10].

Various strategies have been proposed to address these limitations, including alternative loss functions [11], [12] and preconditioning the image vector via neural networks [13] or finite-element grids [14], though these have yet to be validated on experimental UST data. In therapeutic ultrasound, a long-standing solution is to estimate the skull’s acoustic properties from patient-specific CT data and numerically model skull-induced aberration [15], enabling corrections to be applied. One analogous approach in t-UST is to use prior information from the skull is provided to simplify the optimisation problem [7], [8].

However, an important practical challenge is the alignment of skull templates to patient anatomy. Existing approaches typically rely on MRI guidance images to register CT-derived skull models [16], which undermines the clinical utility of t-UST and complicates deployment. Guidance-free registration remains largely unexplored, with only one prior study using reflected ultrasound time-of-flight data and constrained optimisation, restricting rotations and translations to ±3° and ±3 mm [17].

We present Manifold Optimisation for Full-Waveform Inversion (MOFI), a framework for guidance-free skull-template alignment utilising both reflected and transmitted acoustic data employing unconstrained optimisation. MOFI identifies a smooth, low-dimensional manifold that embeds key properties of the FWI image vector, enabling optimisation over physically meaningful transformations while preserving high-fidelity wave propagation modelling. By operating in this reduced space, MOFI effectively preconditions the FWI gradient and improves robustness to large initial misalignments.

1 Method↩︎

In standard FWI, the gradient is computed for each pixel of the discretised sound speed grid using an adjoint state method. The geometry of this gradient is inherently tied to the discretisation of the sound speed grid, where the dimensionality of the update is determined by the number of pixels in the image. However, this geometry is not necessarily optimal for FWI updates, as the sound speed discretisation is primarily chosen to ensure the stability of numerical simulations rather than to reflect the geometric properties of the underlying object. We hypothesise that a more geometrically informed choice of update parameters can significantly improve the update direction, enabling access to solutions that are unattainable through standard FWI. This principle underpins MOFI, which leverages a low-dimensional, geometrically meaningful parameterisation to enhance the efficiency and effectiveness of the inversion process.

1.1 Full-Waveform Inversion↩︎

We consider a standard time-domain acoustic FWI setup. The unknown is the sound-speed distribution \(\mathbf{c}(\mathbf{x})\), \(\mathbf{x} \in \Omega \subset \mathbb{R}^n\). Given appropriate initial and boundary conditions, including a source \(s(\mathbf{x}, t)\), the acoustic wavefield \(\mathbf{u}_p(\mathbf{x}, t)\) satisfies \[\Box(\mathbf{c}) \mathbf{u}_p(\mathbf{x}, t) = s(\mathbf{x}, t),\] where \[\Box(\mathbf{c}) = \frac{1}{\mathbf{c}(\mathbf{x})^2} \frac{\partial^2}{\partial t^2} - \nabla^2,\]

For a point source at position \(\mathbf{x}_s\) with time signal \(\varsigma(t)\), we write \[s(\mathbf{x}, t) = \delta(\mathbf{x} - \mathbf{x}_s)\,\varsigma(t).\] The predicted data are obtained by sampling the wavefield at receiver positions \(\mathbf{x}_r\): \[\mathbf{d}_p(\mathbf{x}_r, t; \mathbf{c}, \mathbf{x}_s, \varsigma) = \delta(\mathbf{x} - \mathbf{x}_r) \, \mathbf{u}_p(\mathbf{x}, t; \mathbf{c}, \mathbf{x}_s, \varsigma).\] Let \(L(\mathbf{c})\) denote the collection of all predicted data for a given model \(\mathbf{c}\) (stacking all receivers, times, and sources into a single data vector), and let \(\mathbf{d}\) be the corresponding observed data. We use the usual least-squares data misfit \[f(\mathbf{c}) = \frac{1}{2} \| \mathbf{d} - L(\mathbf{c}) \|_2^2.\]

The goal of FWI is to find a model \(\mathbf{c}\) that minimises \(f(\mathbf{c})\). This is typically done by iterative gradient-based optimisation. The gradient of \(f\) with respect to \(\mathbf{c}\) is obtained using the adjoint-state method. The gradient is given by: \[\frac{\partial f}{\partial \mathbf{c}} = - \mathbf{u}_p(\mathbf{c})^T \, \frac{\partial \Box(\mathbf{c})}{\partial \mathbf{c}} \, \mathbf{u}_{\text{adj}}(\mathbf{c}),\] where the adjoint wavefield \(\mathbf{u}_{\text{adj}}(\mathbf{c})\) is the solution of the adjoint problem: \[\Box(\mathbf{c})^T \mathbf{u}_{\text{adj}}(\mathbf{c}) = L(\mathbf{c}) - \mathbf{d}.\]

In practice, each iteration consists of:

  1. Forward solve: Compute \(\mathbf{u}_p\) and the predicted data \(L(\mathbf{c})\).

  2. Adjoint solve: Back-propagate the data residual \(L(\mathbf{c}) - \mathbf{d}\) to obtain an adjoint wavefield \(\mathbf{u}_{\text{adj}}\).

  3. Gradient computation: Form the zero-lag of the cross-correlation between the forward and adjoint wavefields to obtain the gradient vector \(\partial f / \partial \mathbf{c}\) at each grid point.

The FWI gradient indicates how the sound-speed model should be updated to reduce the misfit. In conventional FWI, we use this gradient directly to perform a pixel-wise update of \(\mathbf{c}\).

1.2 Manifold Optimisation for FWI↩︎

The standard pixel-wise parameterisation in FWI is inefficient when the mismatch between predicted and observed data arises from geometric deformations, such as global translations or rotations. To address this, MOFI introduces a low-dimensional parameterisation that captures the geometry of the deformation. In this work, we focus on representing the model update through rigid-body motion—namely, translations and rotations of the image — though MOFI can also be extended to non-rigid deformations and other parametrisations. By constraining updates to a geometrically meaningful manifold, MOFI improves the efficiency, accuracy and robustness of the inversion process.

We start from a reference sound-speed image \(\mathbf{c}(\mathbf{x}_{\text{orig}})\) defined on the regular grid \(\mathbf{x}_{\text{orig}}\). Rather than updating \(\mathbf{c}\) directly, we introduce a small set of motion parameters: \[\boldsymbol{\phi} = \{\theta, \delta_1, \delta_2\},\]

which describe a rotation by angle \(\theta\) around a chosen centre and a translation by distance \((\delta_1, \delta_2)\).

These parameters define a rigid transformation of coordinates using the Special Euclidean Matrix Group \(SE(2)\): \[\begin{align} \mathbf{x}_{\text{orig}}' &= T_{\boldsymbol{\phi}}(\mathbf{x}_{\text{orig}}) \\ \begin{pmatrix} x_{\text{orig},1}' \\ x_{\text{orig},2}' \\ 1 \end{pmatrix} &= \begin{pmatrix} \cos(\theta) & -\sin(\theta) & \delta_1 \\ \sin(\theta) & \cos(\theta) & \delta_2 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x_{\text{orig},1} \\ x_{\text{orig},2} \\ 1 \end{pmatrix}. \end{align}\]

Here, \(T_{\boldsymbol{\phi}}\) is the forward transformation matrix that maps the original coordinates \(\mathbf{x}_{\text{orig}}\) to the transformed coordinates \(\mathbf{x}_{\text{orig}}'\).

The aim is to find pixel values for the transformed image. The transformed coordinates \(\mathbf{x}_{\text{orig}}'\) are not usually integers, so \(\mathbf{c}_{\boldsymbol{\phi}}(\mathbf{x}_{\text{orig}}')\) is not defined on a regular pixelated grid. Instead, a pixelated coordinate grid \(\mathbf{x}_{\text{trans}}\) is defined for a “new” image \(\mathbf{c}(\mathbf{x}_{\text{trans}})\). Then, each pixel of the new image is associated with a location \(\mathbf{x}_{\text{trans}}'\) on the original image \[\mathbf{x}_{\text{trans}}' = T_{\boldsymbol{\phi}}^{-1}(\mathbf{x}_{\text{trans}}).\] where the inverse transformation \(T_{\boldsymbol{\phi}}^{-1}\) maps the new coordinate system back to the original coordinate system. The transformed sound-speed image is then parameterised in terms of the original image, given by: \[\mathbf{c}_{\boldsymbol{\phi}}(\mathbf{x}_{\text{trans}}') = \mathbf{c}\bigl(T_{\boldsymbol{\phi}}^{-1}(\mathbf{x}_{\text{trans}})\bigr).\]

Intuitively, \(\mathbf{c}_{\boldsymbol{\phi}}(\mathbf{x}_{\text{trans}}')\) represents the original sound-speed image \(\mathbf{c}(\mathbf{x}_{\text{orig}})\) after it has been rotated by an angle \(\theta\) and translated by \((\delta_1, \delta_2)\). Since the coordinates \(\mathbf{x}_{\text{trans}}' = T_{\boldsymbol{\phi}}^{-1}(\mathbf{x}_{\text{trans}})\) typically do not align with the original grid points, an interpolation algorithm is used to estimate the sound-speed values in the transformed image.

By parametrising the sound-speed image as a function of \(\boldsymbol{\phi}\), we shift the optimisation process from updating individual pixels of \(\mathbf{c}\) to adjusting the geometric transformation parameters. Consequently, the misfit function \(f\) depends on \(\boldsymbol{\phi}\) solely through the transformed model \[f(\boldsymbol{\phi}) = f(\mathbf{c}_{\boldsymbol{\phi}}).\] enabling a more efficient search direction and geometrically meaningful optimisation process.

Using the chain rule, the gradient of \(f\) with respect to the motion parameters is: \[\frac{\partial f}{\partial \boldsymbol{\phi}} = \left( \frac{\partial \mathbf{c}_{\boldsymbol{\phi}}}{\partial \boldsymbol{\phi}} \right)^T \frac{\partial f}{\partial \mathbf{c}_{\boldsymbol{\phi}}}.\] Here, \(\partial f / \partial \mathbf{c}_{\boldsymbol{\phi}}\) is the usual FWI gradient vector with respect to the sound speed (computed by the adjoint-state method), and \(\partial \mathbf{c}_{\boldsymbol{\phi}} / \partial \boldsymbol{\phi}\) is a Jacobian that describes how a change in the motion parameters \(\boldsymbol{\phi}\) changes the sound speed \(\mathbf{c}\) at each pixel.

This formula shows that we can obtain the gradient with respect to the motion parameters \(\boldsymbol{\phi}\) by re-weighting and combining the values of the standard FWI gradient.

1.3 Numerical Implementation↩︎

The implementation is analogous to rigid registration. Each pixel in the sound-speed image is assigned a coordinate, which is a 2D vector representing its position relative to a chosen centre of rotation. Given the current parameters \(\boldsymbol{\phi} = \{\theta, \delta_1, \delta_2\}\), we construct the corresponding rigid transform \(T_{\boldsymbol{\phi}}\):

  1. The transformed sound-speed field \(\mathbf{c}_{\boldsymbol{\phi}}(\mathbf{x}_{\text{trans}}')\) is computed by mapping each grid point \(\mathbf{x}_{\text{trans}}\) back to the original domain via \(\mathbf{x}_{\text{trans}}' = T_{\boldsymbol{\phi}}^{-1}(\mathbf{x}_{\text{trans}})\).

  2. Bilinear interpolation is used to evaluate \(\mathbf{c}\) at coordinates \(\mathbf{x}_{\text{trans}}'\) to form \(\mathbf{c}_{\boldsymbol{\phi}}\).

  3. Solve the forward problem to obtain \(L(\mathbf{c}_{\boldsymbol{\phi}})\).

  4. Compute the data residual \(L(\mathbf{c}_{\boldsymbol{\phi}}) - \mathbf{d}\).

  5. Solve the adjoint problem to form the standard FWI gradient vector \(\partial f / \partial \mathbf{c}_{\boldsymbol{\phi}}\).

  6. Obtain the gradient with respect to the motion parameters by applying \(\frac{\partial f}{\partial \boldsymbol{\phi}} = \left( \frac{\partial \mathbf{c}_{\boldsymbol{\phi}}}{\partial \boldsymbol{\phi}} \right)^T \frac{\partial f}{\partial \mathbf{c}_{\boldsymbol{\phi}}}\).

In practice, we do not provide an explicit expression for the Jacobian \(\partial \mathbf{c}_{\boldsymbol{\phi}} / \partial \boldsymbol{\phi}\). Instead, we use reverse-mode automatic differentiation (Pytorch [18]) to directly compute the product of this Jacobian transpose with the vector \(\partial f / \partial \mathbf{c}_{\boldsymbol{\phi}}\). The FWI gradient with respect to \(\mathbf{c}\) is obtained using the adjoint-state method implemented in Stride [19].

The resulting gradient \(\partial f / \partial \boldsymbol{\phi}\) gives the steepest descent direction in the space of motion parameters. We perform a gradient-based update of \(\boldsymbol{\phi}\), with gradient normalisation and line search to improve stability near the minimum. This allows MOFI to capture the dominant geometric mismatch (translations and rotations) by using a more efficient and accurate search direction, in contrast to the pixel-wise update. Once this coarse alignment is achieved, the standard FWI gradient can, if desired, be used in the usual way to refine smaller-scale features.

2 Results↩︎

2.1 In-silico↩︎

Figure 1: (a) The ground truth model in the “correct” position and (b) the skull template in the “shifted” position. (c) The correct and shifted data – respectively the black and red traces – show a large offset. (d) The FWI gradient captures the leftward shift, but the pixel-wise details are unsatisfactory. (d) By parametrising the image using coordinates, (f) MOFI (Manifold Optimisation for FWI) can represent the update in terms of translation and rotation of the image, providing an update direction (white arrows) that captures the shift.

Figure 1 (a) shows an in-silico phantom representing the sound speed of the head. This phantom is constructed from the MIDA model [20] and the IT’IS database of tissue properties [21]. The shape of the grid is (560, 450) and isotropic spacing of the grid is \(0.5 \, \mathrm{mm}\).

Figure 1 provides an initial validation of our methodology. A large offset between (a) the “correct” and (b) the “shifted” p-wave velocity models causes (c) a large offset between the corresponding data (black and red traces). Figure 1 (d) shows that the standard FWI gradient is sensitive to this shift, but applying this gradient pixel-wise leads to inefficient updates. MOFI replaces pixel-wise updates with a small set of motion parameters (translations and rotations of the image), which provides a more efficient update direction. Figure 1 (e) is an initial validation of our hypothesis—the white arrows show the direction of the MOFI update, which is the reverse of the shift induced in Figure 1 (b).

Figure 2: In-silico wavefield data are simulated from (a) the ground truth sound speed. (b) The skull template at the initial position produces (c) a failed FWI reconstruction. (d) MOFI identifies the correct translation (\mathrm{\delta1, \delta2}) and rotation (\mathrm{\theta}). (d) The MOFI aligned skull template produces (f) a successful FWI reconstruction.

Figure 2 demonstrates that MOFI is able to accurately locate the skull template when the ground truth sound speed undergoes a substantial shift in position and orientation (Fig 2 (a)). Relative to the the original image, the rotation is \(10^\mathrm{o}\)—relative to the vertical—about the centre of the image. The translation is \((14, -5.85) \, \mathrm{mm}\).

Figure 2 (b) shows a starting model where the skull template has the correct sound speed and shape, but its position (i.e. translation and rotation) is incorrect. Figure 2 (c) shows that FWI cannot reconstruct the brain from this starting model. Except for the position of the skull, this example is an inverse crime—the modelling and reconstruction physics are identical, the exact source wavelet and transducer positions are known, and the data are not cycle-skipped as the inversion starts from 100kHz.

Figure 2 (d) shows the three parameters \(\boldsymbol{\phi}\) of SE(2) during the optimisation process. The translation starts at \((0,0) \, \mathrm{mm}\) and the rotation starts at \(0^\mathrm{o}\). In total there are 500 iterations. The last 200 iterations oscillate around a stationary point, which give a translation \((13.9 \pm 0.3, -6.8 \pm 0.3) \, \mathrm{mm}\) and rotation \(9.8 \pm 0.3^\mathrm{o}\). With respect to the known displacement, the horizontal translation is incorrect by 1mm, or two pixels, the other values are accurate to less than a pixel. The stationary point is reached after approximately 200 iterations, which is computationally feasible. Furthermore, because the group has a limited number of parameters, the problem is—in a sense—overdetermined, which means that a small batch size (10% of the emitters) can be used. Figure 2 (e) shows the skull template at the optimal position, which is determined by the optimisation on SE(2). The FWI reconstruction in Figure 2 (f) confirms that the starting model is now correctly located.

2.2 In vitro↩︎

Figure 3: image.

FWI reconstructions of lab data are substantially more challenging. First, the impulse response of the transducer is unknown. This means that the spatial response of the emitter and sensitivity of the receiver are difficult to determine with precision. Second, the physics of the problem are unknown. To reduce computation time, we solve the acoustic wave equation with constant density and zero attenuation. In reality, materials have variable density, attenuate wave energy, and may support elastic or anisotropic physics. Finally, the properties of the starting model are unknown. Different sources in the literature provide conflicting estimates of the skull sound speed, which suggests that these estimates have wide errors. Furthermore, the shape of the skull can only be determined imprecisely from CT or MRI data (typical resolution of \(0.5 \, \mathrm{mm}\)), and its position is only known approximately despite careful placement.

Figure 3 demonstrates that MOFI can optimise the position of the in vitro skull template. Figure 3 (column 1) shows the parameter values \(\boldsymbol{\phi}=(\delta_1, \delta_2, \theta)\) for iterations 1, 50, 150, and 210. These correspond to the sound speed images in Figure 3 (column 2), which show the template moving smoothly from the initial to the final position.

Figures 3 (columns 3-5) provides a detailed view of the data misfit during the optimisation process. The observed acoustic data are in black, the predicted acosutic data are in red.

Broadly, transmitted energy in the predicted data (panel 2, red traces) is smoothly shifted to adjacent receivers, but the arrival time remains unchanged. In contrast, reflected energy (panel 1&3, red traces, second arrival) is shifted to earlier arrival times. The first arrival in panels 1&3 travels directly through water and remains unchanged throughout the optimisation. In detail:

  • Figure 3, Iteration 1: The predicted and observed data exhibit significant misalignment across all panels.

  • Figure 3, Iteration 50: In panels 1 and 3, the reflected predicted data (second arrivals) shift to earlier times, moving toward the observed data. In panel 1, a complete alignment in phase and amplitude is achieved, while the transmitted data in panel 2 remain largely unchanged.

  • Figure 3, Iteration 150: In panel 3, the predicted reflected data overshoot the observed data, while in panel 1 (on the opposite side of the array) the arrival time continues to steadily decrease. In panel 2, the timing of the transmitted data remains stable, though energy is redistributed to adjacent receivers.

  • Figure 3, Iteration 210: The data in all three panels align closely: transmitted energy is correctly redistributed, and reflected energy is appropriately time-shifted. Minor misalignments persist, particularly for reflected arrivals in panel 1, which are cycle-skipped.

Figure 4: (a) MOFI identifies a translation \boldsymbol{\delta}=(7.2 \pm 0.2, 0.7 \pm 0.2) \, \mathrm{mm} and a rotation \mathrm{\theta} = 10.71 \pm 0.07^\mathrm{o}. (b) The skull template at the initial position produces (c) a failed FWI reconstruction. (d) The MOFI aligned skull template provides (e) a successful FWI reconstruction.

Figure 4 demonstrates that MOFI itentifies an accurate position and orientation for the skull template. Figure 4 (a) shows that the optimal translation is \((7.2 \pm 0.2, 0.7 \pm 0.2) \, \mathrm{mm}\) and the optimal rotation is \(10.71 \pm 0.07^\mathrm{o}\), measured from the final 300 iterations. The stability of the translation and rotation values between iteration 200 and iteration 500 confirms that the optimisation process has converged. Figure 4 (b) shows the skull template in its initial position. Previous work [8] shows that this template has approximately the correct sound speed and shape, but none of the values are known precisely. Figure 4 (c) confirms the incorrect placement because the FWI reconstruction has failed. Figure 3 (d) shows the skull template in the optimal position, identified at iteration 500. Figure 4 (e) shows an accurate FWI reconstruction, which confirms that the starting model is now correctly positioned.

3 Discussion↩︎

Our approach delivers strong results whenever an approximate sound-speed template is available and requires alignment. Ultrasound neurotechnologies rely on accurate skull template alignment, as misalignment compromises the precision of numerical refocusing corrections. Rigid MOFI addresses this challenge by enabling automatic skull template alignment with sub-pixel accuracy. This innovation reduces the need for manual intervention in FWI neuroimaging and provides a robust alternative — or complement — to MRI-guided image registration in transcranial focused ultrasound therapy.

The effectiveness of MOFI stems from its use of the Special Euclidean Group SE(2) to constrain the full-waveform inversion (FWI) optimisation process. SE(2) is part of a broader class of mathematical structures known as Lie Groups, which combine algebraic group properties with the smooth geometry of manifolds. These groups are particularly suited to modelling continuous transformations—such as rotations and translations—while preserving geometric symmetries like angles and distances. As a result, Lie Groups offer a natural framework for optimisation problems. Their utility is well-documented in fields such as image registration and robotics. For instance, diffeomorphic non-rigid image registration, which leverages the Diffeomorphism group, significantly outperforms free-form deformation methods when handling large deformations [22]. In robotics, manifold optimisation using Lie Groups enhances the accuracy of trajectory mapping [23][25].

Another perspective emerges from convex optimisation, where computational complexity scales polynomially with the number of dimensions [26]. Although full-waveform inversion (FWI) is non-convex, the convex case implies that the complexity of non-convex optimisation also increases with dimensionality. Conventional FWI optimisation algorithms operate on a discretised sound-speed grid, where the grid is designed to ensure the stability of the numerical simulator rather than optimisation efficiency. In Figures 2 and 4, the discretised sound-speed grids consist of 252,000 and 672,000 pixels, respectively — each pixel representing an individual dimension in the optimisation space.

Thus, by constraining the optimisation to just three dimensions—rotation and translation—MOFI circumvents the challenges posed by high-dimensional discretised grids, such as those in Figures 2 and 5, where each pixel introduces an additional optimisation variable. Dimensionality reduction by manifold optimisation is not unique to FWI [14], [27] but, through the use of SE(2), the Lie Group structure of MOFI ensures a smooth optimisation manifold that retains algebraic group properties. This contrasts with level set methods [27], which lack a group structure and may struggle to enforce specific deformations. Furthermore, unlike methods that invert for finite-element grid control points [14], MOFI guarantees invertibility, preventing grid points from self-intersecting and preserving the image’s topology. As a result, MOFI delivers a robust and efficient optimisation process, readily adaptable to any discretisation (e.g. finite-difference or finite-element ) while maintaining stability for large deformations.

Looking ahead, we aim to explore more expressive Lie Groups, such as the Diffeomorphism Group, and investigate how concepts from Lie Groups could enhance Deep Learning approaches. The design of Deep Learning models often results in non-smooth representations — exemplified by the challenges in training Generative Adversarial Networks. Approaches that regularise inverse problems but do not explicitly enforce smoothness, such as the Deep Image Prior [28] or Score-based Generative Models [29], may therefore benefit from incorporating this principle. Similarly, Normalising Flows, which employ smooth, invertible transformations, have already demonstrated their utility in FWI [13], [30].

It is well documented that FWI struggles to converge when the predicted and observed data arrive out-of-phase [10]. MOFI is much more resilient than the typical requirement that the phase difference is less than \(90^{\textrm{o}}\), but it will eventually struggle if this distance becomes too large. Therefore, MOFI also benefits from the multi-scale method [31]. Similarly, the current loss function uses the L2-norm, but alternative formulations from the FWI literature could further expand the basin of attraction. In terms of optimisation parameters, both fixed-step and explicit line-search methods were evaluated, with explicit line-search providing greater stability at the optimum. Gradient normalisation was found to be essential, and using the norm of the first gradient proved most effective, whereas normalising by the maximum of the current gradient led to instability. Future work should also focus on extending this approach to 3D in-vivo datasets to assess its clinical efficacy.

4 Conclusion↩︎

Our findings demonstrate that a skull sound-speed template can be accurately and automatically located using only acoustic data. This eliminates the need for manual alignment of the “starting model” in ultrasound neuroimaging and potentially enhances the image-to-image alignment used during ultrasound therapy. We greatly improve the convergence properties of FWI by leveraging the Special Euclidean Lie Group, which reduces the dimensionality of the FWI optimisation problem while preserving a smooth optimisation manifold. This strategy has been successful in other fields, particularly non-rigid image registration, and suggests that more expressive Lie Groups could entirely eliminate the need for patient-specific prior information in wave-based image-reconstruction algorithms. Together, these findings support the use of MOFI as a principled pathway toward practical, guidance-free transcranial Ultrasound Tomography.

5 Acknowledgements↩︎

This work was funded by EPSRC grant number EP/X033651/1. Oscar Calderon-Agudo was funded by MRC, Future Leader Fellowship (FLF), MR/V024086/1. We would like to thank the team at Sonalis Imaging Limited for their support and insightful discussions.

References↩︎

[1]
H. Estrada, T. Deffieux, J. Robin, M. Tanter, and D. Razansky, “Imaging the brain by traversing the skull with light and sound,” Nature Biomedical Engineering, vol. 9, no. 10, pp. 1574–1590, Oct. 2025, doi: 10.1038/s41551-025-01433-5.
[2]
Y. Tufail et al., “Transcranial Pulsed Ultrasound Stimulates Intact Brain Circuits,” Neuron, vol. 66, no. 5, pp. 681–694, Jun. 2010, doi: 10.1016/j.neuron.2010.05.008.
[3]
K. Hynynen, N. McDannold, N. Vykhodtseva, and F. A. Jolesz, “Noninvasive MR Imaging–guided Focal Opening of the Blood-Brain Barrier in Rabbits,” Radiology, vol. 220, no. 3, pp. 640–646, Sep. 2001, doi: 10.1148/radiol.2202001804.
[4]
W. J. Elias et al., _eprint: https://www.nejm.org/doi/pdf/10.1056/NEJMoa1300962“A Pilot Study of Focused Ultrasound Thalamotomy for Essential Tremor,” New England Journal of Medicine, vol. 369, no. 7, pp. 640–648, Aug. 2013, doi: 10.1056/NEJMoa1300962.
[5]
S. Pichardo, V. W. Sin, and K. Hynynen, “Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls,” Physics in Medicine & Biology, vol. 56, no. 1, p. 219, Dec. 2010, doi: 10.1088/0031-9155/56/1/014.
[6]
J. Robertson, E. Martin, B. Cox, and B. E. Treeby, “Sensitivity of simulated transcranial ultrasound fields to acoustic medium property maps,” Physics in Medicine & Biology, vol. 62, no. 7, p. 2559, Mar. 2017, doi: 10.1088/1361-6560/aa5e98.
[7]
L. Guasch, O. Calderon Agudo, M.-X. Tang, P. Nachev, and M. Warner, “Full-waveform inversion imaging of the human brain,” npj Digital Medicine, vol. 3, no. 1, pp. 1–12, Mar. 2020, doi: 10.1038/s41746-020-0240-8.
[8]
T. C. Robins et al., “Dual-Probe Transcranial Full-Waveform Inversion: A Brain Phantom Feasibility Study,” Ultrasound in Medicine & Biology, vol. 49, no. 10, pp. 2302–2315, Oct. 2023, doi: 10.1016/j.ultrasmedbio.2023.06.001.
[9]
T. Mitcham et al., “Transcranial ultrasound tomography for brain imaging: Ex vivo results and potential for stroke imaging,” Medical Physics, vol. 52, no. 10, p. e18090, Oct. 2025, doi: 10.1002/mp.18090.
[10]
J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” GEOPHYSICS, vol. 74, no. 6, pp. WCC1–WCC26, Nov. 2009, doi: 10.1190/1.3238367.
[11]
T. van Leeuwen and F. J. Herrmann, “Mitigating local minima in full-waveform inversion by expanding the search space,” Geophysical Journal International, vol. 195, no. 1, pp. 661–667, Oct. 2013, doi: 10.1093/gji/ggt258.
[12]
M. Warner and L. Guasch, “Adaptive waveform inversion: Theory,” GEOPHYSICS, vol. 81, no. 6, pp. R429–R445, Nov. 2016, doi: 10.1190/geo2015-0387.1.
[13]
R. Orozco et al., ADS Bibcode: 2024JOSS....9.6554OInvertibleNetworks.jl: A Julia package for scalable normalizing flows,” The Journal of Open Source Software, vol. 9, p. 6554, Jul. 2024, doi: 10.21105/joss.06554.
[14]
P. Marty, C. Boehm, and A. Fichtner, “Shape optimization for transcranial ultrasound computed tomography,” in Medical Imaging 2023: Ultrasonic Imaging and Tomography, Apr. 2023, vol. 12470, pp. 77–88, doi: 10.1117/12.2654328.
[15]
J.-F. Aubry, M. Tanter, M. Pernot, J.-L. Thomas, and M. Fink, “Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans,” The Journal of the Acoustical Society of America, vol. 113, no. 1, pp. 84–93, Jan. 2003, doi: 10.1121/1.1529663.
[16]
N. McDannold, G. T. Clement, P. Black, F. Jolesz, and K. Hynynen, “Transcranial Magnetic Resonance ImagingGuided Focused Ultrasound Surgery of Brain Tumors: Initial Findings in 3 Patients,” Neurosurgery, vol. 66, no. 2, p. 323, Feb. 2010, doi: 10.1227/01.NEU.0000360379.95800.2F.
[17]
M. A. O’Reilly, R. M. Jones, G. Birman, and K. Hynynen, _eprint: https://aapm.onlinelibrary.wiley.com/doi/pdf/10.1118/1.4960362“Registration of human skull computed tomography data to an ultrasound treatment space using a sparse high frequency ultrasound hemispherical array,” Medical Physics, vol. 43, no. 9, pp. 5063–5071, 2016, doi: 10.1118/1.4960362.
[18]
A. Paszke et al., PyTorch: An Imperative Style, High-Performance Deep Learning Library,” in Advances in Neural Information Processing Systems, 2019, vol. 32, Accessed: Aug. 26, 2022. [Online].
[19]
C. Cueto et al., “Stride: A flexible software platform for high-performance ultrasound computed tomography,” Computer Methods and Programs in Biomedicine, vol. 221, p. 106855, Jun. 2022, doi: 10.1016/j.cmpb.2022.106855.
[20]
M. I. Iacono et al., MIDA: A Multimodal Imaging-Based Detailed Anatomical Model of the Human Head and Neck,” PLOS ONE, vol. 10, no. 4, p. e0124126, Apr. 2015, doi: 10.1371/journal.pone.0124126.
[21]
P. Hasgall, E. Neufeld, M. Gosselin, A. Klingenbock, and N. Kuster, ITIS Database for thermal and electromagnetic parameters of biological tissues.” IT’IS Foundation, 2015, doi: 10.13099/VIP21000-04-0.
[22]
M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes, “Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms,” International Journal of Computer Vision, vol. 61, no. 2, pp. 139–157, Feb. 2005, doi: 10.1023/B:VISI.0000043755.93987.aa.
[23]
R. W. Brockett, “Robotic manipulators and the product of exponentials formula,” in Mathematical Theory of Networks and Systems, 1984, pp. 120–129, doi: 10.1007/BFb0031048.
[24]
B. K. P. Horn, “Closed-form solution of absolute orientation using unit quaternions,” JOSA A, vol. 4, no. 4, pp. 629–642, Apr. 1987, doi: 10.1364/JOSAA.4.000629.
[25]
J. Solà, J. Deray, and D. Atchuthan, arXiv:1812.01537 [cs]“A micro Lie theory for state estimation in robotics.” arXiv, Dec. 2021, doi: 10.48550/arXiv.1812.01537.
[26]
S. P. Boyd and L. Vandenberghe, Google-Books-ID: mYm0bLd3fcoCConvex Optimization. Cambridge University Press, 2004.
[27]
A. Kadu, T. van Leeuwen, and W. A. Mulder, “Salt Reconstruction in Full-Waveform Inversion With a Parametric Level-Set Method,” IEEE Transactions on Computational Imaging, vol. 3, no. 2, pp. 305–315, Jun. 2017, doi: 10.1109/TCI.2016.2640761.
[28]
D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep Image Prior,” 2018, pp. 9446–9454, Accessed: May 30, 2025. [Online].
[29]
Y. Song, L. Shen, L. Xing, and S. Ermon, arXiv:2111.08005 [eess]“Solving Inverse Problems in Medical Imaging with Score-Based Generative Models.” arXiv, Jun. 2022, doi: 10.48550/arXiv.2111.08005.
[30]
C. Sun, A. Malcolm, R. Kumar, and W. Mao, “Enabling uncertainty quantification in a standard full-waveform inversion method using normalizing flows,” GEOPHYSICS, vol. 89, no. 5, pp. R493–R507, Sep. 2024, doi: 10.1190/geo2023-0755.1.
[31]
C. Bunks, F. M. Saleck, S. Zaleski, and G. Chavent, “Multiscale seismic waveform inversion,” Geophysics, vol. 60, no. 5, pp. 1457–1473, Oct. 1995, doi: 10.1190/1.1443880.