Optoacoustic Model-Based Inversion Using Anisotropic Adaptive Total-Variation Regularization

TL;DR

This paper introduces an adaptive anisotropic total-variation regularization framework for optoacoustic image reconstruction, improving boundary preservation and contrast in noisy or incomplete data scenarios.

Contribution

The novel adaptive anisotropic total-variation regularization scheme enhances boundary preservation and contrast over traditional methods in optoacoustic tomography.

Findings

Outperforms total-variation-$L_1$ scheme in simulations and experiments.

Better preserves complex boundaries in blood-vessel images.

Enhances image contrast in noisy and incomplete data conditions.

Abstract

In optoacoustic tomography, image reconstruction is often performed with incomplete or noisy data, leading to reconstruction errors. Significant improvement in reconstruction accuracy may be achieved in such cases by using nonlinear regularization schemes, such as total-variation minimization and $L_1$-based sparsity-preserving schemes. In this paper, we introduce a new framework for optoacoustic image reconstruction based on adaptive anisotropic total-variation regularization, which is more capable of preserving complex boundaries than conventional total-variation regularization. The new scheme is demonstrated in numerical simulations on blood-vessel images \textcolor{black} {as well as on experimental data} and is shown to be more capable than the total-variation-$L_1$ scheme in enhancing image contrast.

Equations

$$\frac{{\partial }^{2}p(\mathbf{r},t)}{\partial t^2}-c^2{\nabla }^{2}p(\mathbf{r},t)=\Gamma H_r(\mathbf{r})\frac{\partial H_t(t)}{\partial t}$$

$$p_{\delta }(\mathbf{r},t)=\frac{\Gamma }{4\pi c}\frac{\partial }{\partial t}\frac{\delta (|\mathbf{r}-{\mathbf{r}}^{\prime }|-ct)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$p(\mathbf{r},t)=\frac{\Gamma }{4\pi c}\frac{\partial }{\partial t}{\int }_{|\mathbf{r}-{\mathbf{r}}^{\prime }|=ct}\frac{H_r({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$\mathbf{p}=\mathbf{Mu}$$

$${\mathbf{u}}^{\ast }=\arg \min ||\mathbf{p}-\mathbf{Mu}|{|}_2^2,$$

$${\mathbf{u}}^{\ast }=({\mathbf{M}}^{\mathbf{T}}\mathbf{M}{)}^{-1}{\mathbf{M}}^T\mathbf{p}.$$

$${\mathbf{u}}^{\ast }=\arg \min \Vert \mathbf{p}-\mathbf{Mu}{\Vert }_2^2+\lambda \Vert \mathbf{Lu}{\Vert }_2^2$$

$${\mathbf{u}}^{\ast }=\arg \min ||\mathbf{p}-\mathbf{Mu}|{|}_2^2+\mu \Vert \mathbf{\Phi u}{\Vert }_1,$$

$$\Vert \mathbf{u}{\Vert }_{TV}=\sum _n\sqrt{|u_{x,y}^n-u_{x-1,y}^n{|}^2+|u_{x,y}^n-u_{x,y-1}^n{|}^2},$$

$${\mathbf{u}}^{\ast }=\arg \min ||\mathbf{p}-\mathbf{Mu}|{|}_2^2+\alpha \Vert \mathbf{u}{\Vert }_{TV}$$

$${\mathbf{u}}^{\ast }=\arg \min ||\mathbf{p}-\mathbf{Mu}|{|}_2^2+\mu \Vert \mathbf{\Phi u}{\Vert }_1+\alpha \Vert \mathbf{u}{\Vert }_{TV}.$$

$${\mathcal{J}}_{A^{2}TV}(\mathbf{u})={\int }_{\Omega }\Vert A(\mathbf{x})\nabla \mathbf{u}(\mathbf{x}){\Vert }_{2}d\mathbf{x}={\int }_{\Omega }\Vert {\nabla }_A\mathbf{u}(\mathbf{x}){\Vert }_{2}d\mathbf{x}$$

$${\mathcal{J}}_{\rho }(\nabla {\mathbf{u}}_{0;\sigma })={\kappa }_{\rho }\ast \left(\nabla {\mathbf{u}}_{0;\sigma }\otimes \nabla {\mathbf{u}}_{0;\sigma }\right),$$

$${\mathcal{J}}_{\rho }(\mathbf{x})=\left(\begin{array}{cc}{\kappa }_{\rho }\ast u_{x;\sigma }^2(\mathbf{x})& {\kappa }_{\rho }\ast u_{x;\sigma }(\mathbf{x})u_{y;\sigma }(\mathbf{x})\\ {\kappa }_{\rho }\ast u_{x;\sigma }(\mathbf{x})u_{y;\sigma }(\mathbf{x})& {\kappa }_{\rho }\ast u_{y;\sigma }^2(\mathbf{x})\end{array}\right)$$

$${\mathcal{J}}_{\rho }(\mathbf{x})={\kappa }_{\rho }\ast \left(\begin{array}{cccc}{u}_{x;\sigma }^2(\mathbf{x})& u_{x;\sigma }(\mathbf{x})u_{y;\sigma }(\mathbf{x})& u_{x;\sigma }(\mathbf{x})u_{z;\sigma }(\mathbf{x})& \\ u_{x;\sigma }(\mathbf{x})u_{y;\sigma }(\mathbf{x})& u_{y;\sigma }^2(\mathbf{x})& u_{y;\sigma }(\mathbf{x})u_{z;\sigma }(\mathbf{x})& \\ u_{x;\sigma }(\mathbf{x})u_{z;\sigma }(\mathbf{x})& u_{y;\sigma }(\mathbf{x})u_{z;\sigma }(\mathbf{x})& u_{z;\sigma }^2(\mathbf{x})& \end{array}\right).$$

$${\mathcal{J}}_{\rho }=VD{V}^{-1}$$

$$D=\left(\begin{array}{cc}{\mu }_1& 0\\ 0& {\mu }_2\end{array}\right);V=(v_1|v_2)$$

$$A=V\tilde{D}{V}^{-1},$$

$$\tilde{D}=\left(\begin{array}{cc}c({\mu }_1/{\mu }_{1,avg};k)& 0\\ 0& 1\end{array}\right).$$

$$c(s;k)=\left\{\begin{array}[]{rcl}1&\mbox{,}&s\leq 0\\ 1-exp(-\frac{c_{m}}{(\frac{s}{k})^{m}})&\mbox{,}&s>0\end{array}\right.$$

$$\tilde{D}=\left(\begin{array}{ccc}c({\mu }_1/{\mu }_{1,avg};k)& 0& 0\\ 0& c({\mu }_2/{\mu }_{1,avg};k)& 0\\ 0& 0& 1\end{array}\right).$$

$${\mathbf{u}}^{\ast }=\arg \underset{\mathbf{u}}{\min }{\mathcal{J}}_{A^{2}TV}(\mathbf{u})+\frac{\lambda }{2}\Vert \mathbf{Mu}-\mathbf{p}{\Vert }_2^2,$$

$$\text{MAD}=\frac{1}{N}\Vert {\mathbf{u}}_{\text{orig}}-{\mathbf{u}}^{\ast }{\Vert }_1,$$

$$\begin{gathered} {\mathbf{u}}^{\ast }=\arg \underset{\mathbf{u}\in X}{\min }\mathcal{E}(\mathbf{u})= \\ \arg \underset{\mathbf{u}\in X}{\min }{\mathcal{J}}_{A^{2}TV}(\mathbf{u})+\frac{\lambda }{2}\Vert \mathbf{Mu}-\mathbf{p}{\Vert }_2^2. \end{gathered}$$

$$\arg \underset{\mathbf{u}\in X}{\min }F(K\mathbf{u})+G(\mathbf{u}).$$

$$\begin{array}[]{cc}F(\mathbf{y})=\|\mathbf{y}\|_{1},G(\mathbf{u})=\frac{\lambda}{2}\|\mathbf{M}\mathbf{u}-\mathbf{p}\|_{2}^{2}\\ K=\nabla_{A}=\textrm{grad}_{A},K^{*}=\nabla_{A}^{T}=\textrm{div}_{A}\end{array}.$$

$$\arg \underset{\mathbf{u}\in X}{\min }\arg \underset{\mathbf{z}\in Y}{\max }<K\mathbf{u},\mathbf{z}>+G(\mathbf{u})-F^{\ast }(\mathbf{z}).$$

$$\arg \underset{\mathbf{u}\in X}{\min }\arg \underset{\mathbf{z}\in Y}{\max }<{\nabla }_A\mathbf{u},\mathbf{z}>+\frac{\lambda }{2}\Vert \mathbf{Mu}-\mathbf{p}{\Vert }_2^2-{\delta }_P(\mathbf{z}),$$

$$F^{\ast }(\mathbf{z})={\delta }_P(\mathbf{z}),$$

$${\delta }_P(\mathbf{z})=\{\begin{array}{ll}0,& \mathbf{z}\in P\\ \infty ,& \mathbf{z}\notin P\end{array},P=\{\mathbf{z}\in Y|\Vert \mathbf{z}{\Vert }_{\infty }\le 1\}.$$

Full text

††journal: Photoacousticsmyfootnotemyfootnotefootnotetext: These authors contributed equally to this work

1 Introduction

Optoacoustic tomography (OAT) is a hybrid imaging modality capable of visualizing optically absorbing structures with ultrasound resolution at tissue depths in which light is fully diffused [1, 2, 3, 4, 5]. The excitation in OAT is most often performed via high-energy short pulses whose absorption in the tissue leads to the generation of acoustic sources via the process of thermal expansion. The acoustic signals from the sources are measured over a surface that partially or fully surrounds the imaged object and used to form an image of the acoustic sources, which generally represents the local energy absorption in the tissue [6]. Since hemoglobin is one of the strongest absorbing tissue constituents, optoacoustic images often depict blood vessels and blood-rich organs, such as the kidneys [7] and heart [8].

Numerous algorithms exist for reconstructing optoacoustic images from measured tomographic data [6]. In several imaging geometries, analytical formulae exist that may be applied directly on the measured data in either time [9, 10] or frequency domain [11] to yield an exact reconstruction. The popularity of analytical formulae may be attributed to the simplicity of their implementation and low computational burden [12] However, analytical algorithms are not exact for arbitrary detection surfaces or detector geometries and lack the possibility of regularizing the inversion in the case of noisy or incomplete data. In such cases, it is often preferable to use model-based algorithms, in which the relation between the image and measured data is represented by a matrix whose inversion is required to reconstruct the image.

In the last decade, numerous regularization approaches have been demonstrated for model-based image reconstruction. The most basic approach is based on energy minimization and includes techniques such as Tikhonov regularization [13] and truncated singular-value decomposition [14]. In these techniques, a cost function on the image or components thereof is used to avoid divergence of the solution in the case of missing data, generally without making any assumptions on the nature of the solution. More advanced approaches to regularization exploit the specific properties of the reconstructed image. Since natural images may be sparsely represented when transformed to an alternative basis, e.g. the wavelet basis, using nonlinear cost functions that promote sparsity in such bases may be used for denoising and image reconstruction from missing data [15, 16, 17, 18, 19, 20, 21]. In images of blood vessels, in which the boundaries of the imaged structures may be of higher importance than the texture in the image, total-variation (TV) minimization has been shown to enhance image contrast and reduce artifacts [22, 23, 24]

In TV regularization, the cost function regularizer is the $L_{1}$ norm of the image gradient, which is generally lower for images with sharp, yet very localized, variations than for images in which small variations occur across the entire image. Therefore TV regularization enhances boundaries and reduces texture, where over-regularization may lead to almost piecewise constant images, which are often referred to as cartoon-like. While TV regularization is capable of accentuating the boundaries of imaged objects, it does not treat all boundaries the same. In particular, boundaries with short lengths will lead to a lower TV cost function than boundaries with long lengths. Thus, complex, non-convex boundaries may be rounded by TV regularization to the closest convex form. Recently, Wang et al. have shown that if the directionality of the TV functional is adapted to the image features, TV regularization may be applied for optoacoustic reconstruction of objects with non-convex boundaries without distorting the boundaries [25].

In this paper, we demonstrate a new regularization framework for model-based optoacoustic image reconstruction that overcomes the limitations of TV regularization and is compatible with objects with complex non-convex boundaries. In our scheme, an adaptive anisotropic total variation (A2TV) cost function is used, in which the cost function is determined by the specific geometry of the imaged objects [26]. In particular, the A2TV cost function wishes to minimize the total variation in directions that are orthogonal to the boundary of the object. In contrast to [25], where the boundaries were calculated using geometrical considerations limited to 2D images, the A2TV framework developed in this work is based on eigenvalue decomposition of the image structure tensor, which may be applied in higher dimensions. The proposed formalism in the current study is based on a recent work by part of the authors which is concerned with nonlinear spectral analysis of the A2TV functional [26]. The work of Ref. [26] examines shapes which are perfectly preserved under A2TV regularization. Earlier works concerning TV have shown that only convex rounded shapes of low curvature are preserved [27]. For A2TV, however, it is shown in [26] that a parameter controlling the local extent of directionality is directly related to the degree of convexity (in the sense of [28]) and to the curvature magnitude of structures which are preserved. Thus, with appropriate parameters, long vessels of complex-nonconvex structure can be better regularized, keeping the original structure intact.

The performance of A2TV regularization was tested numerically in this work in numerical simulations for complex images of blood vessels and experimentally on 2D phantoms. The simulations were performed for the cases of noisy data and missing data and compared to unregularized reconstructions as well as to TV-$L_{1}$ reconstructions [15]. In both the numerical and experimental reconstructions, A2TV significantly increased the image contrast and was more capable than TV-$L_{1}$ in preserving non-convex structures when strong regularization was performed. In the experimental reconstructions, A2TV achieved a higher level of contrast enhancement of weak structures than the one achieved by TV-$L_{1}$.

The rest of the paper is organized as follows: in Section 2 we give the theoretical background for OAT image reconstruction. Section 3 introduces the framework of A2TV and the A2TV algorithm for OAT image reconstruction developed in this work. The simulation results are given in Section 4, while the experimental ones are given in Section 5. We conclude the paper with a Discussion in Section 6.

2 Optoacoustic image reconstruction

2.1 The forward problem

The acoustic waves in OAT are commonly described by a pressure field $p(\mathbf{r},t)$ that fulfills the following wave equation[29]:

[TABLE]

where $c$ is the speed of sound in the medium, $t$ is time, $\mathbf{r}=(x,y,z)$ denotes position in 3D space, $p(\mathbf{r},t)$ is the generated pressure, $\Gamma$ is the Grüneisen parameter, and $H_{r}(\mathbf{r})H_{t}(t)$ is the energy per unit volume and unit time. The spatial distribution function of energy deposited in the imaged object, $H_{r}(\mathbf{r})$, is referred to in the rest of the paper as the optoacoustic image.

The analysis of (1) for an optoacoustic point source at $\mathbf{r}^{\prime}$, i.e. $H_{r}(\mathbf{r})=\delta(\mathbf{r-r^{\prime}})$ , may be performed in either time or frequency domain. In the time domain, a short-pulse excitation $H_{t}(t)=\delta(t)$ is used. In this case, the solution to Eq. (1) is given by [30]

[TABLE]

For a general image $H_{r}(\mathbf{r})$, the solution for $p(\mathbf{r},t)$ may be obtained by convolving $H_{r}(\mathbf{r})$ with the expression in Eq. (2), which yields:

[TABLE]

Since the relation between the measured pressure signals, or projections, and the image is linear, it may be represented by a matrix relation in its discrete form:

[TABLE]

where $\mathbf{p}$ and $\mathbf{u}$ are vector representations of the acoustic signals and originating image respectively, and $\mathbf{M}$ is the model matrix that represents the operations in Eq. (3). In our work, the images is given on a two dimensional grid, and the measured pressure signals are also two dimensional, where one dimension represents the projection number, and the other time. An illustration of the image grid and projection and their respective mapping to the vectors $\mathbf{u}$ and $\mathbf{p}$, is shown in Fig. 1. Briefly, the vector $\mathbf{u}$ is divided into sub-vectors, each representing the image values for a given column, whereas the vector $\mathbf{p}$ is divided to sub-vectors, each of which represents the time-domain pressure signal for a given location of the acoustic detector.

The ith column of the matrix $\mathbf{M}$ represents the set of acoustic signals generated by a pixel corresponding to the location of the ith entry in the vector $\mathbf{u}$. Accordingly, to calculate the matrix $\mathbf{M}$ one needs to define time domain signal for a given detector location expected for a discrete pixel. Since the operations in Eq. (3) relate to continuous, rather than discrete images, one first needs to define the continuous representation of a single pixel, and then calculate its respective time-domain signals. For example, in [31] it was assumed that the image value was constant over each of the square pixels, leading to an image $H_{r}(\mathbf{r})$ that is piece-wise constant. While simple to implement, a simple piece-wise uniform model for $H_{r}(\mathbf{r})$ includes discontinuities that lead to significant numerical errors owing to the derivative operation in Eq. (3). In the current work, we use the model of [32], in which the image $H_{r}(\mathbf{r})$ represented by a linear interpolation between its grid points.

2.2 The inverse problem

While several approaches to OAT image reconstruction exist, we focus herein on image reconstruction within the discrete model-based framework described in the previous sub-section, which involves inverting the matrix relation in Eq. (4) to recover $\mathbf{u}$ from $\mathbf{p}$. The most basic method to invert Eq. (4) is based on solving the following optimization problem:

[TABLE]

where $\mathbf{u}^{*}$ is the solution and $||\cdot||_{2}$ is the $L_{2}$ norm. A unique solution to Eq. (5) exists, which is given by the Moore–Penrose inverse:

[TABLE]

Alternatively, Eq. (5) may be solved via iterative optimization algorithms. In particular, since the matrix $\mathbf{M}$ is sparse, efficient inversion may be achieved by the LSQR algorithm [33]

In many cases, the measured projection data $\mathbf{p}$ is insufficient to achieve a high-quality reconstruction of $\mathbf{u}$ that accurately depicts its morphology. For example, when the density or coverage of the projections is too small, the matrix $\mathbf{M}$ may become ill-conditioned, leading to significant, possibly divergent, image artifact. In other cases, $\mathbf{M}$ may be well conditioned, but the measurement data may be too noisy to accurately recover $\mathbf{u}$. In both these cases, regularization may be used to improve image quality by incorporating previous knowledge on the properties image in the inversion process.

One of the simplest forms of regularization is Tikhonov regularization, in which an additional cost function is added [13]:

[TABLE]

where $\mathbf{L}$ is weighting matrix. In the simplest form of Tikhonov regularization, $\mathbf{L}$ is equal to the identity matrix, i.e. $\mathbf{L}=\mathbf{I}$, thus putting a penalty on the energy of the image. The value of the regularization parameter $\lambda>0$ controls the tradeoff between fidelity and smoothness, where over-regularization may lead to the smearing of edges and texture in the image.

An alternative to the energy-minimizing cost function of Tikhonov are sparsity-maximizing cost functions. Since natural images may be sparsely represented in an alternative basis, e.g. the wavelet basis, a cost function that promotes sparsity may reduce reconstruction errors. Denoting the transformation matrix by $\mathbf{\Phi}$, one wishes that $\mathbf{\Phi u}$ be sparse, i.e. that most of its entries be approximately zero. In practice, sparsity is often enforced by using the $L_{1}$ norm because of its compatibility of optimization algorithms [16, 17, 18, 19, 20, 21]. Accordingly, the inversion is performed by solving the following optimization problem:

[TABLE]

where $\mu>0$ is the regularization parameter, controlling the tradeoff between sparsity and signal fidelity. When over regularization is performed, compression artifacts may appear in the reconstructed image.

Sparsity may be enforced not only on alternative representations of the image, but also on image variations. The discrete TV cost function approximates the $l_{1}$ norm on the image gradient, and is given by

[TABLE]

where $u_{x,y}^{n}$ is the $n$th entry of the vector $\mathbf{u}$, and the subtraction of ”1” to $x$ or $y$ in the subscript corresponds to an entry of $\mathbf{u}$ that is respectively shifted in relation to $u_{x,y}^{n}$ by one pixel in the $x$ or $y$ direction of the 2D image. The inversion using TV is thus given by [22, 23, 24]

[TABLE]

where $\alpha>0$ is the regularization parameter. In the case of TV minimization, over-regularization may lead to cartoon-like images and rounding of complex boundaries into convex shapes.

In some cases, it is beneficial to promote both sparsity of the image in an alternative basis and TV minimization. In such cases, the optimization problem is as follows [15]:

[TABLE]

We will refer to the regularization described in Eq. (11), as TV-$L_{1}$ regularization.

3 Adaptive Anisotropic Total Variation

3.1 The Functional

We would like to use a regularizer that is adapted to the image in such a way that it regularizes more along edges (level-lines of the image) and less across edges (in the direction of the gradient). This idea has been introduced for nonlinear scale-space flows by Weickert [34] in the anisotropic diffusion formulation. However there is no known functional associated with anisotropic diffusion, and it is therefore not trivial to include an anisotropic-diffusion operation in our inverse-problem formulation. A more recent study of Grasmair et al. [35] uses a similar adaptive scheme within a TV-type formulation. In the study of [26], a comprehensive theoretical and numerical analysis was performed for adaptive-anisotropic TV (A2TV). It was shown that stable structures can be non-convex and in addition can have high curvature on the boundaries. Illustrations of the stable sets characterize the TV and A2TV regularizers are shown in Fig. 2. The degree of anisotropy directly controls the degree of allowed nonconvexity and the upper bound on the curvature. We adopt the formulation of [26], in which the mathematical underpinnings of A2TV are described in detail.

Let the A2TV functional be defined by

[TABLE]

where $\mathbf{x}=(x,y)\in\Omega$ denotes position in 2D space, $A(\mathbf{x})\in\mathbb{R}^{2\times 2}$ is a tensor or a spatially adaptive matrix and $\nabla_{A}=A(\mathbf{x})\nabla$ is an “adaptive gradient”. This functional is convex and can be optimized by standard convex solvers given the tensor $A(\mathbf{x})$. We now turn to the issue of how this tensor is constructed to allow a good regularization of vessel-like structures.

3.1.1 Constructing the tensor $A(\mathbf{x})$

We assume to have a rough approximation of the image to be reconstructed. This can be done, for instance, by having an initial least-square non-regularized approximation or a standard TV reconstruction, as in Eq. (10). We thus have an initial estimation $\mathbf{u}_{0}$. The tensor $A(\mathbf{x})$ determines the principle and the secondary directions of the regularization at each point and their magnitude. The construction of $A(\mathbf{x})$ is performed using $\mathbf{u}_{0}$ according to the following principles: In regions in which $\mathbf{u}_{0}$ is relatively flat, i.e. $\nabla\mathbf{u}_{0}$ almost vanishes, $A(\mathbf{x})$ should resemble the identity matrix and have no preferred direction, thus leading to the conventional TV regularizer. In regions with dominant edges, $A(\mathbf{x})$ should capture the principle axes of the edge. See Fig. 3 for an illustration of $A(\mathbf{x})$.

Mathematically, the tensor $A(\mathbf{x})$ is defined by adapting the eigenvalues of a smoothed structure tensor of a smoothed image $\mathbf{u}_{0;\sigma}$ (with a Gaussian kernel of standard deviation $\sigma$), defined by,

[TABLE]

where $\kappa_{\rho}$ is a Gaussian kernel with a standard deviation of $\rho$, $*$ denotes an element-wise convolution and $\otimes$ denotes an outer product. The Gaussian kernel’s standard deviations $\sigma$ and $\rho$ are chosen in coordination with the noise level of the image and the smallest object resolution within the image. High standard deviation values might diminish the effect of small objects on the structure tensor while low values might lead to the unwanted scenario in which reconstruction errors have a significant effect on the structure tensor. The explicit expressions for $\mathcal{J}_{\rho}$ in the 2D and 3D cases are respectively given by

[TABLE]

and

[TABLE]

The structure tensor matrix has eigenvectors corresponding to the direction of the gradient and tangent at each point $\mathbf{x}$; and eigenvalues corresponding to the magnitude of each direction. In order to preserve structure, we should change the relation between those eigenvalues so that for flat-like areas in the image we will smooth the image in an isotropic way, while for edge-like areas, we will perform more smoothing in the tangent direction rather the gradient one. For 2D, we begin by looking in the eigen-decomposition of the structure-tensor,

[TABLE]

Where $V$ is a matrix whose columns, $v_{1},v_{2}\in\mathbb{R}^{2}$, are the eigenvectors of $\mathcal{J}_{\rho}$, and $D$ is a diagonal matrix with eigenvalues in the diagonal,

[TABLE]

Assuming $\mu_{1}\geq\mu_{2}$. The spatially adaptive matrix $A(\mathbf{x})$, used in Eq. (12), is constructed from a modification of the eigenvalue matrix, denoted by $\tilde{D}$, as follows [36]:

[TABLE]

where $\tilde{D}$ is given by

[TABLE]

In Eq. (19), $\mu_{1,avg}$ is the average value of $\mu_{1}$ across the image and $c(\cdot;\cdot)$ is a function of two parameters defined as follows:

[TABLE]

where the chosen values for the parameter are the ones recommended in Ref. [34]: $c_{m}=3.31488$ and $m=4$. The parameter $k\leq 1$ determines which regions in the image will be regularized anisotropically and is chosen based on the desired level of anisotropy in the reconstructed image, as shown in Sections 4 and 5. In pixels in which $s\ll k$, i.e. $\mu_{1}/\mu_{1,avg}\ll k$, we will obtain $c\approx 1$, and $\tilde{D}$ will be reduced to the unitary matrix. Accordingly, for regions in which the image gradient is sufficiently small from the average image gradient, as regulated by $k$, the A2TV functional (Eq. (12)) is reduced to the standard, isotropic TV functional (Eq. (9)). In the rest of the image, where $c$ is sufficiently smaller than 1, regularization is performed more strongly in the direction of the eigenvector $v_{2}$, i.e. the direction in which the image gradient is smaller, thus enhancing the anisotropy in those regions

It is worth noting that in the 3D case, assuming $\mu_{1}\geq\mu_{2}\geq\mu_{3}$, the only modification to the analysis above is that $\tilde{D}$ accepts the following form:

[TABLE]

The structure of $\tilde{D}$ in Eq. (21) will be highly anisotropic for tube-like structures and will enforce variation-reducing regularization along the structure length (eigenvector $v_{3}$), while maintaining low regularization in the cross-section plane (eigenvectors $v_{1}$ and $v_{2}$).

3.2 Reconstruction based on A2TV

The reconstruction based on the A2TV minimizes the following functional:

[TABLE]

where $\mathbf{M}$ is the model matrix and $\mathbf{p}$ is the acoustic pressure wave. The solution of which is done by the modified Chambolle-Pock projection algorithm [37] described in Appendix A.

In the process of minimization, the tensor $A$ is initialized as the identity matrix for all $\mathbf{x}\in\mathbb{R}^{2}$, which reduces the A2TV energy to the TV one as it performs the diffusion isotropically. The tensor $A$ is then updated according to the initial solution $u_{0}$. This is repeated until numerical convergence is reached.

We note that while the energy of the A2TV is convex for a fixed tensor $A(\mathbf{x})$, it is not convex when $A(\mathbf{x})$ is adaptive and depends on the imaged object. Thereby, we do not have a mathematical proof of convergence. Nonetheless, it has been shown both in our work [26] and in Refs. [36, 35] that heuristically, both the image $u$ and the tensor $A(\mathbf{x})$ converge.

4 Numerical Simulations

In this section, we demonstrate the performance of A2TV-based inversion for the circular detection geometry illustrated in Fig. 4. The simulations were performed on a 2D vascular image of a mouse retina, obtained via confocal microscopy. The vasculature image, shown in Fig. 5a, was represented over a square grid with a size of 256$\times$256 pixels with pixel size of 0.1$\times$0.1 mm. The projections were simulated over a 270-degree semi-circle with a radius of 4 cm that surrounded the object, in accordance with conventional OAT systems [4]. A magnification of four square regions of the image is shown in Fig. 5b.

The reconstructions were performed using the conventional $L_{2}$-based regularization-free approach (Eq. (5)) performed via LSQR, TV-$L_{1}$ regularization (Eq. (11)), and the proposed A2TV regularization (Eq. (22)) method. Since the scaling of the model matrix $\mathbf{M}$ (Eq. (4)) depends on the exact implementation of its construction [32], we normalized $\mathbf{M}$ by $\frac{1}{160}\sqrt{\|\mathbf{M}\|_{\infty}\|\mathbf{M}\|_{1}}$ to assure that the regularization parameters are independent of the scaling of $\mathbf{M}$. To assess the quality of the reconstructions, we used the mean absolute distance (MAD) given by the following equation:

[TABLE]

where $N$ is the number of pixels in the image.

Two cases were tested: In the first case, zero-mean Gaussian noise with a standard deviation of 0.6 times the maximum value of $\mathbf{p}$ was added to the projection. The number of projections was chosen to be 256, corresponding to the geometry found in state-of-the-art optoacoustic systems [38] and sufficient for the accurate reconstruction of the tested image in the noiseless case. In the second case, the number of projections was reduced to 32, which is half the number of projections used in low-end optoacoustic systems characterized by reduced lateral resolution [38]. Accordingly, 32 projections are insufficient for producing detailed optoacoustic images using conventional reconstruction techniques. In all the examples, the number of iterations was chosen to be sufficiently high to achieve convergence.

Fig. 6 and Fig. 7 respectively show the reconstructions obtained using TV-$L_{1}$ and A2TV, performed with 3000 iterations, for the cases of additive Gaussian noise. In both figures, 9 reconstructions are shown, corresponding to a scan in the regularization parameters. In the case of TV-$L_{1}$, $\mu$ and $\alpha$ represent the strength of the $L_{1}$ and TV regularization terms Eq. (11), whereas in the case of A2TV, $\lambda$ represents the strength of the fidelity term in Eq. (22) with respect to the A2TV term and $k$ determines the strength of the anisotropy, where lower values of $k$ correspond to higher anisotropy. In the A2TV reconstructions, the standard deviations of the smoothing kernels (Eq. 14) were $\sigma=1.5$ pixels and $\rho=3$ pixels. For all the reconstructions, the MAD values appear on the top-right corner of the image. In Fig. 8, we show a comparison between the regularization-free LSQR reconstruction (Fig. 8a) and the TV-$L_{1}$ (Fig. 8b) and A2TV (Fig. 8c) of Fig. 6d and Fig. 7e, respectively, which correspond to the regularization parameters that achieved the lowest MAD values. The middle panel of the figure (Fig. 8d-f) shows a magnification of 4 patches taken from the reconstructions, whereas the bottom panel (Fig. 8g) presents a 1D slice taken over the yellow line in Fig. 8b and Fig. 8c. While both TV-$L_{1}$ and A2TV significantly improved the reconstruction quality, in A2TV more of the noise-induced texture between the blood vessels could be removed without damaging the structure of the blood vessels, thus leading to a lower MAD.

The reconstructions for the case of 32 projections are presented in a similar manner to the case of noisy data. Fig. 9, obtained with 1000 iterations, and Fig. 10, obtained with 1500 iterations, respectively show the TV-$L_{1}$ and A2TV reconstructions for a scan regularization parameters, whereas Fig. 11 shows a comparison between the LSQR and TV-$L_{1}$ and A2TV reconstructions that achieved the lowest MAD (Fig. 9d and Fig. 10e, respectively). All the A2TV reconstructions were obtained with $\sigma=1.5$ pixels and $\rho=1$ pixel. Fig. 11 shows that both TV-$L_{1}$ and A2TV eliminated the streak artifacts that appeared with in the LSQR reconstruction, where the lowest MAD was achieved by the TV-$L_{1}$ reconstruction. Since both regularization methods eliminated the streak artifacts, the lower MAD achieved by TV-$L_{1}$ is a result of its ability to better preserve texture within the blood vessels, whereas in the A2TV much of the blood-vessel texture was lost. Indeed, when examining the 1D slice in Fig. 11g, it is easy to see that the TV-$L_{1}$ reconstruction captures the variations within the blood vessels better, whereas in the A2TV reconstruction, these variations are smoothed.

In both the cases studied, A2TV exhibited a higher capability than TV-$L_{1}$ to perform regularization without harming non-convex structures. Even in the case in which TV-$L_{1}$ achieved a lower MAD, the higher ability of A2TV to preserve the fine details of the blood vessel morphology can be observed when comparing Fig. 11e to Fig. 11f. Additionally, in all the examples, one can observe that when TV-$L_{1}$ was performed with a high level of TV regularization (bottom row in Fig. 6 and Fig. 9) significant smearing of the blood vessels was observed. In contrast, in the A2TV reconstructions, the smearing owing to over-regularization (low values of $\lambda$) could be diminished by increasing the anisotropy in the regularization, i.e. reducing the value of $k$. For the lowest values of $k$, higher levels of regularization (low values of $\lambda$) created undesirable anisotropic vessel-like texture in the reconstructed images.

5 Experimental Results

To further validate the suitability of A2TV regularization for OAT image reconstruction, we tested its performance on experimental data. The optoacoustic setup comprised an optical parametric oscillator (OPO), which produced nanosecond optical pulses with an energy 30 mJ at a repetition rate of 100 Hz and at a wavelength of $\lambda=680$ nm. (SpitLight DPSS 100 OPO, InnoLas Laser GmbH, Krailling, Germany). The OPO pulses were delivered to the imaged object using a fiber bundle (CeramOptec GmbH, Bonn, Germany). Ultrasound detection was performed by a 256-element annular array (Imasonic SAS, Voray sur l’Ognon, France) with a radius of 4 cm, and an angular coverage of 270 degrees, comparable to the geometry shown in Fig. 4. The ultrasound detectors were cylindrically focused to a plane, approximating a 2D imaging scenario.

The imaged object was a transparent agar phantom which contained four intersecting hairs. A photo of the phantom is shown in Fig.12a, where Fig.12b shows 4 magnified parts of the phantom. The phantom preparation involved mixing 1.3% (by weight) agar powder (Sigma-Aldrich, St. Louis, MO) in boiling water and pouring the solution in a cylindrical mold until solidification. To assure that all four hair strands lie approximately in the same plane, we first prepared a clear cylindrical agar phantom, on which the hairs were placed; additional agar solution was then poured on the structure to seal the hairs.

The optoacoustic image was reconstructed from the measured data using TV-$L_{1}$ regularization with 3000 iterations and A2TV regularization with 6000 iterations, $\sigma=1.5$ pixels, and $\rho=1$ pixel. Figs. 13 and 14 respectively show the images obtained via TV-$L_{1}$ and A2TV reconstructions for a set of regularization parameters. As in the previous section, over-regularization in the TV-$L_{1}$ led to significant loss of structure, whereas for the A2TV the ability to capture the image morphology under strong regularization (low $\lambda$) was improved when the anisotropy was increased via low values of $k$.

Fig.15 compares between the unregularized LSQR reconstruction (Fig. 15a) and the TV-$L_{1}$ (Fig. 15b) and A2TV (Fig. 15c) reconstructions, respectively taken from Figs. 13e and 14e. The second row in the figure (Figs. 15d-15f) shows a magnification of 4 patches from the three reconstructions of the top row (Figs. 15a-15c), whereas the bottom panel (Fig.15g) shows a 1D slice from the vertical yellow line in Figs. 15b and 15c. To allow an easy comparison, the 1D slices were normalized by their maximum values. We note that the negative values in the reconstruction are a common result of the limited detection bandwidth of the ultrasound detector [6]. The figure clearly shows that the A2TV reconstruction obtained the highest image quality, in particular for the weak hair structure at the image bottom. In particular, the 1D slice shows that the bottom hairs that appear around the position of 1 cm achieve a peak-to-peak signal over 4-times higher in the A2TV reconstruction than in the TV-$L_{1}$ reconstruction.

6 Discussion

In this work, a novel regularization framework was developed for OAT image reconstruction. The new framework is based on an A2TV cost function, which represents a generalization of the conventional TV functional, that is compatible with objects that possess complex, nonconvex boundaries. In contrast to TV, the A2TV cost function is an adaptive functional whose form depends on the characteristics of the image. When using A2TV, one first roughly defines the boundaries of the objects in the image. Then, one uses these boundaries to determine the directions in which the gradients are applied on the image.

Similar to TV, A2TV is most appropriate as a regularizer in optoacoustic image reconstruction when the images are comprised of objects with well-defined boundaries. However, A2TV is useful also when these boundaries are incompatible with TV regularization due to their complexity. A common category of optoacoustic images that fits the above description is blood-vessel images. Since blood is a major source of contrast in optoacoustic imaging, OAT systems often produce images that are dominated by a complex structure of interwoven blood vessels. In particular, high-resolution images of the micro-vasculature are characterized by a complex network of arterioles, venues, and capillaries with extremely complex, nonconvex boundaries.

In our current implementation, A2TV required setting 4 parameters: $\sigma$, $\rho$, $\lambda$, and $k$. The first two parameters, $\sigma$ and $\rho$, determined the image smoothing used in calculating the image gradients. While smoothing reduces the noise, thus limiting the effect of reconstruction errors on the detection of the principle axes, only moderate smoothing may be used without the risk of merging the boundaries of different objects in the image. Therefore, in all our examples, the smoothing was performed with Gaussian kernels whose standard deviations, $\sigma$ and $\rho$, were 3 pixels or less. While the choice of $\sigma$ and $\rho$ depended on the level of noise and artifacts in the regularization-free reconstruction, the parameters $\lambda$ and $k$ were determined by the amount of regularization desired in the reconstructed image, where $k$ determined the amount of anisotropy and $\lambda$ determined the strength of the regularizer. In our examples, performing over-regularization ($\lambda=0.0001$) led to image smearing in the isotropic case ($k=1$) and vessel-like artifacts in the case of high anisotropy ($k=0.01$). While such artifacts are undesirable, it is worth noting that their presence did not obscure the underlying image morphology, whereas over-regularization in TV-$L_{1}$ led to loss of image details.

We compared the performance of the A2TV algorithm TV-$L_{1}$ algorithm for both numerically simulated data and experimental data. In the numerical simulations, an image of blood vessels was reconstructed for the cases additive noise and sparse sampling of the projection data. In the experimental reconstructions, the imaged object was four intersecting hair strands whose structure emulated the morphology of blood vessels. In both the numerical and experimental examples, A2TV demonstrated a higher ability to preserve the blood-vessel morphology for high regularization parameters. Nonetheless, in the numerical example in which the reconstructions were performed with a low number of projections, TV-$L_{1}$ regularization achieved a lower MAD owing to the loss of texture in the A2TV reconstruction. In the experimental example, A2TV led to a considerable improvement in the contrast in the weak structures of the image in comparison to the TV-$L_{1}$ reconstruction.

The reconstruction performance demonstrated in this work suggests that A2TV may be a useful tool for improving the ability of optoacoustic systems to perform vasculature imaging, which is a major application in the field. Since the texture within the blood vessels is affected by the random distribution of the red-blood cells within the blood vessels, its elimination by A2TV may be considered as an acceptable price for better visualization of the blood-vessel morphology. In deep-tissue OAT systems, sub-millimeter vasculature imaging has been suggested as a potential diagnostic tool and has been demonstrated in the human extremities [39, 40] and breast [41]. We note that while some OAT systems can also produce images of the tissue bulk, characterized by low frequencies and representative of the density of the microvasculature and fluence map, such systems require transducers capable of detecting ultrasound frequencies considerably below 1 MHz [6]. In high-resolution OAT systems, operating at frequencies above 1 MHz and capable of reaching resolutions better than 100 $\mu$m [42, 43], only signals from blood vessels may be detected. Finally, when performing optoacoustic imaging at resolutions better than 10 $\mu$m, e.g. using raster-scan optoacoustic mesoscopy (RSOM) [44], the image is dominated by the microvasculare and generally lacks any bulk component associated with blood.

We note that the formalism of A2TV, in which the structure tensor matrix analysis is performed via eigenvalue decomposition Eq. (13) enables its adaptation to higher image dimensions. TV regularization has been recently performed in 4D optoacoustic reconstruction that included 3 spatial dimensions and time [45]. Since the variation of the pixels in time is generally different than the one space, it may be expected that A2TV can further improve image fidelity in such cases.

7 Acknowledgement

We are thankful to Anna M. Randi for providing information on the image used in Section 4.

Appendix A Solving the generalized ROF model

A2TV based reconstruction amounts to solving the following generalized ROF model,

[TABLE]

Where $\mathbf{M}$ is the model matrix, $\mathbf{p}$ is the acoustic pressure wave and $\mathbf{u}^{*}$ is the reconstructed image. The algorithm used for inversion is Chambolle-Pock, which solves the following primal form [37]:

[TABLE]

Where, in our case,

[TABLE]

Let X and Y be two finite-dimensional real vector spaces and $K$, a linear operator $K:X\rightarrow Y$ and its hermitian adjoint $K^{*}:Y\rightarrow X$. The appropriate saddle point problem is,

[TABLE]

Where $\mathbf{z}\in Y$ is the dual variable. In our case,

[TABLE]

where

[TABLE]

and

[TABLE]

But, because of the appearance of the model matrix, we shall adjust the saddle point problem to

[TABLE]

using

[TABLE]

Here, the original saddle point problem is transformed into,

[TABLE]

Where

[TABLE]

Now, in order to use the algorithm, we shall use the proximal operator of $F^{*}$ and $G$,

[TABLE]

The algorithm is given in Fig. 1, where in our case, we normalize the model matrix $\mathbf{M}$ by $\frac{1}{160}\sqrt{\|\mathbf{M}\|_{\infty}\|\mathbf{M}\|_{1}}$, where

$\sqrt{\|\mathbf{M}\|_{\infty}\|\mathbf{M}\|_{1}}$ is an approximation to the model matrix $M$ Lipchitz constant. $L=L_{\mathbf{M}}+L_{\nabla}=160+8=168$, $\gamma=0.7\lambda$ and $\tau_{0}=0.5$