Iterative Image Correction Scheme for Optoacoustic Tomography

An iterative algorithm to enhance the quality of optoacoustic images is proposed based on the Banach fixed point theorem. Numerical simulation of the propagation of ultrasound waves and image reconstruction for media with characteristics close to soft biological tissues has been performed. The problem of eliminating distortions and artifacts in optoacoustic images has been modeled for four iterative schemes. The efficiency of the approach has been demonstrated even for a few iterations.


INTRODUCTION
Noninvasive nondisturbing diagnostic tools such as magnetic resonance and X-ray computed tomography are currently widely used in scientific research and clinical practice. The world scientific community appreciated the great importance of these tools by awarding their creators Nobel Prizes in physiology and medicine. However, the hard fields and radiation used in these tools are potentially hazardous to human health, which severely limits their application. Optical coherence tomography and optical diffusion tomography have considerable potential in terms of safety and obtaining reliable diagnostically meaningful information. However, due to the strong absorption and scattering of light waves in biological tissues, purely optical methods can only be used to efficiently visualize structures only in subsurface domains. Therefore, one of the most promising and safest trends in medical diagnostics is the approach that uses the optoacoustic (OA) effect.
Optoacoustic imaging uses the phenomenon of optical radiation absorption by inhomogeneities in a biotissue, which results in its heating and subsequent thermal expansion. The thermal expansion of optical absorbers leads to generation of ultrasound waves, which are detected by an array located on the surface of a studied object. Since the propagation and formation characteristics of OA signals are determined by the thermophysical, acoustic, and optical properties of a medium, the relationship between these characteris-tics makes it possible to use the detected ultrasound signals to quantify the properties of a medium by solving the inverse problem of optoacoustics.
In addition to biomedical safety, an important advantage of OA imaging is that this approach combines the capabilities of high-contrast optical absorption and deep ultrasound penetration. In this case, the quality of an image reconstructed from the ultrasound signals detected on the tissue surface is mainly determined by an optoacoustic reconstruction algorithm. Distortions and artifacts in a reconstructed image can occur due to both the features of the reconstruction method, and obstacles and noises of various natures inevitable in practical studies. The OA images typically contain an intense non-stationary background with artifacts similar to a pattern of useful signals, the signal-to-noise ratio is usually low, and the digital image is of low quality, and has a few digitizing levels, specklelike pattern, and fuzzy boundaries. Therefore, the problem of eliminating distortions and artifacts in the reconstructed OA images is very relevant for the efficient use of this method in clinical practice and scientific research [1].
The aim of this work is to develop and test in silico algorithms for eliminating distortions and artifacts from reconstructed two-(2D) and three-dimensional (3D) OA images. First, the direct and inverse OA problems are formulated together with a description of the proposed iterative algorithm for correcting reconstructed OA images. The next section describes in more detail the algorithm, numerical experiment, and the results of simulating the reconstruction of the 2D and 3D OA images. In the Conclusions, the results and prospects of the developed algorithm are presented.

ITERATIVE OPTOACOUSTIC RECONSTRUCTION ALGORITHM
Optoacoustic imaging is based on the effect of thermoelasticity, when optical radiation absorption by inhomogeneities of a medium results in conversion of optical energy into thermal. At a moderate density of the released energy, no phase transformations occur in the absorption region, and thermal optical absorption results in the generation of ultrasound waves.
The direct problem of OA tomography is to determine the pressure field p(r, t) using the known distribution of heat sources H(r, t) excited by a short light (laser) pulse (r is the spatial coordinate of a point, t is the time) (Fig. 1).
The spatiotemporal dependence p(r, t) being sought in an acoustically homogeneous infinite medium is determined by the following equation [2]: (1) where c is the ultrasound velocity, β is the isobaric expansion coefficient, and c p the specific heat capacity at constant pressure when the heat source H(r, t) is presented as the product of the spatial distribution of the absorbed energy and the time dependence of the The initial acoustic pressure occurred due to the absorption of pulsed radiation by optical inhomogeneities at the time instant t = 0 can be expressed in the form p 0 (r) = ΓQ(r), where Γ is the dimensionless Gruneisen parameter, which characterizes the efficiency of the OA transformation of absorbed light into sound.
In this case, the solution of the direct problem is [2] (2) where V' is the volume containing the distributed OA sources. The inverse problem of optoacoustics is to reconstruct the initial acoustic pressure p 0 (r, t = 0) = p 0 (r) using the pressure signals p s (r, t) detected on surface S of volume V' [3].
The algorithms for reconstructing the heat source distribution in a medium can be divided into several classes: the Fourier-domain algorithm, the timereversal algorithm, and the back-projection method [4,5].
The back-projection procedure [6] is the classical method of the OA reconstruction. It can be implemented either in the space-time domain or in the Fourier domain for several detection modes in plane [6], cylindrical [6], or spherical [7] geometries. In this case, the solutions are valid only for a closed perfect

Ultrasound detectors
Ultrasound waves Optical inomogeneity surface (each surface point serves as a detector). In addition, it is usually considered that the target object is located in an infinitely homogeneous medium without dispersion with a constant ultrasound velocity, absorption coefficient and density.
In practice, these conditions are not fulfilled, which leads to distortions of the reconstructed images. The features of these distortions are different when using different reconstruction methods, different locations of detectors, and different geometries of a reconstructed object.
To compensate distortions and overcome the strong dependence of image quality on these factors, we developed an iterative distortion reduction scheme for OA images based on the Banach fixed point theorem. Let where x is the unknown input image, and f(⋅) is the mapping (operator) that transforms input image x into result y of the solution of the inverse problem of optoacoustics. In this case, the operator f(⋅) is the result of the subsequent application of the operators f 1 (⋅) and is the operator that determines the solution of the direct problem of optoacoustics (Eq. (2)), and f 2 (⋅) is an operator that determines the solution of the inverse problem. According to the Banach fixed point theorem it is known that if f is a contraction mapping of the set F ⊂ (M, ρ) into itself, (M, ρ) is a complete metric space, and F is a closed set, then there exists, and moreover exactly one, fixed point x* ∈ F of mapping f (if the point is fixed, then f(x*) = x*).
In this case, contraction mapping means the following mapping f: In accordance with this definition, for the reconstruction problem F = R m×n is a digital image consisting of m × n pixels (for the 2D case), and ρ(x, y) = ||x-y|| is the Euclidean distance. It can be shown that F ⊂ (M,ρ) is a complete metric space, for which the Banach fixed point theorem is valid [8]. In this case, the fixed point x* ∈ F can be found by the method of successive approximations: x n = f(x n-1 ), where the zero approximation x 0 ∈ F is an arbitrary point of the metric space.
Thus, a refined solution to the inverse problem of optoacoustics can be obtained by forming a sequence of images {I q } such that I q = R(I q-1 ), where R(⋅) is an operator defined by the mappings f 1 (⋅) and f 2 (⋅) and transforming the original pressure distribution p 0 (r) to its zero approximation . In other words, the solution of the inverse problem of optoacoustics p 0 (r) is taken as the zero approximation, and occurring distortions are removed during an iterative process.
In the proposed approach, two algorithms of successive approximations were tested: the algorithm based on the Picard theorem [9] and the algorithm proposed in [10] (formulas (4) and (5), respectively): To quantify the quality of image restoration, there were used such criteria as the signal-to-noise ratio: the relative error (here, is the qth iteration of the algorithm), and the structure similarity index (SSIM), which characterizes the closeness of images X and Y in terms of brightness, contrast, and pattern. The structure similarity index is considered as an unofficial standard for estimating image quality, when a reference exists. In general, the value of SSIM(X,Y) is determined by the formula where x, y are the coordinates of a pixel, is the brightness functional, is the contrast functional, is the pattern functional, and μ x , μ y , σ x , σ y , and σ xy are the local mean values, standard deviations, and the cross covariance for corresponding images, respectively. The SSIM index takes values from -1 to 1. The value 1 is obtained if the compared images are equal [11].
Note that the implementation of the proposed iterative scheme is possible in two modes. The first mode is implemented using formulas (4) and (5). In the second mode, the sequence of images {I q } is formed by an iterative process I q = I(I q-1 ), where I(⋅) is an operator defined by the same mappings f 1 (⋅) and f 2 (⋅), but acting in a different order, i.e. f(⋅) = f 1 (f 2 (⋅)). In this case, the input image is not the pressure distribution reconstructed by solving the inverse problem, but the Here, we use formulas similar to (4) and (5), but for the operators I 1 (⋅) and I 2 (⋅), respectively.
The next section presents the results of testing both iterative schemes.

TESTING OF AN ALGORITHM
The k-Wave software code operating in the MAT-LAB toolbox was used for numerical simulation of the problem of propagation of acoustic waves. It allows on to model systems with acoustic sources and detectors of arbitrary shapes and sizes. In this case, the numerical model is based on the transition to the k-space. The spatial gradients in this space are calculated using a fast Fourier transform scheme. The time gradients are calculated using an adjusted k-space difference scheme [12].
A numerical model was set close in its characteristics to soft biological tissues: a homogeneous medium with the density ρ 0 = 1020 kg/m 3 and ultrasound velocity c 0 = 1510 m/s. The problem was solved for 2D and 3D cases.
In the 2D case, the numerical phantoms (objects for the OA reconstruction) were a disk and a twodimensional model of the vascular tree; in the 3D case, it was a three-dimensional numerical model of the aorta with an aneurysm. The physical sizes of the 2D and 3D samples were 4.6 × 4.6 mm and 10.3 × 10.3 × 5.3 mm, respectively. In the 2D case, detectors were located linearly on the upper surface of a rectangular sample; in the 3D case, detectors were distributed over the upper plane of a parallelepiped. The back-projection algorithm with Fourier transform was used to reconstruct the given objects. Figure 2a shows the OA source p 0 (x, y) in the form of a circular disk and the result of its reconstruction implemented using the k-Wave algorithm is shown in Fig. 2b. The shape of the reconstructed image is mapped quite accurately. However, the image is aggravated with arc artifacts touching the reconstructed image; the intensity of the reconstructed image is much lower than that of the original sample, and its edges are diffused. This is especially noticeable in  Fig. 2d, where the given OA source p 0 (x, y) is shown in an isometric projection (the source intensity is plotted along the applicate axis). The reconstruction result is shown in Fig. 2d in the foreground.
The two middle images and in Fig. 2d represent the fourth iteration of the scheme I q = R 2 (I q-1 ) and the 4-th iteration of the scheme I q = I 2 (I q-1 ), respectively.
The results of the implementation of the algorithm can be seen in more detail in Fig. 3. It shows the linear profiles of the reconstructed images and together with the original image p 0 (x, y). Figure 3a shows zero approximation of circular disc p 0 . Arrows mark the location and orientation of the vertical and horizontal sections of the reconstructed image. The resulting linear profiles of original object p 0 (bold black line) and its reconstructions are shown in Figs. 3b and 3c. Clearly, the implemented algorithm resulted in that the edges of reconstructed and corrected image become sharper, its intensity practically coincides with the original, and errors and reconstruction artifacts present in zero approximation are eliminated, which leads to a considerable increase in the signal-to-noise ratio. Similar results were obtained for more complex 2D and 3D objects. Figure 4 shows the results of the implemented algorithm for a 2D phantom of a vascular tree. It should be emphasized here that the proposed iterative correction algorithm efficiently operates for both horizontal and vertical linear structures, which are especially poorly reconstructed with a limited view and linear arrangement of detectors only on the upper surface of the sample (see, for example, vertical segments of branches in Fig. 4b).
A similar result is obtained in 3D simulation, when a 3D numerical model of the aorta with an aneurysm was used as the original object for reconstruction (Fig. 5).
The original image of the aorta is shown in Fig. 5a, the result of the k-Wave reconstruction ( in our notations) is shown in Fig. 5b, and corrected modification of the reconstructed image is shown in Fig. 5c. More detailed differences between the reconstructed images are shown in Fig. 5d, which presents the results of reconstruction in the central x-z section of the simulated 3D vascular tree (line notations are the same as in Fig. 3). Similar to the 2D case, the algorithm considerably improves the quality of the reconstructed object and much more accurately reproduces its boundaries, edges, and intensity.
Quantitative estimates of the reconstruction quality for the described reconstructed objects in terms of the SSIM index are summarized in the Table 1. Clearly, the proposed algorithm considerably  The rate of convergence of the proposed iterative schemes can be estimated from the results shown in Fig. 6. As it was mentioned above, such characteristics as the SSIM index and relative error E were used as criteria of the efficiency of the image correction. Figure 5 shows that with an increase in iteration number q, both proposed iterative schemes I q = R(I q-1 ) and I q = I(I q-1 ) provide a steady enhancement in image quality both in terms of SSIM and E. In both cases, faster conver-  gence is provided by the algorithms I q = R 2 (I q-1 ) and I q = I 2 (I q-1 ). Figure 6 shows the dependences for different iterative schemes when reconstructing a plane circular disk. Similar patterns were observed when reconstructing other model objects.

CONCLUSIONS
The main objective of this work was to develop and study a numerical algorithm to correct artifacts and distortions resulting from the reconstruction of OA images. The problem was to develop an algorithm capable of compensating for the features of a particular reconstruction method. The proposed iterative scheme of the correction of OA images is based on the Banach fixed point theorem.
To study the efficiency of the proposed algorithm, a numerical model of an OA experiment was constructed that simulates a biological tissue with an embedded object to be reconstructed. The linear (2D case) or planar (3D case) detecting acoustic arrays located on the surface of the studied samples were considered.
The k-Wave software code operating in the MATLAB toolbox was used to solve the inverse problem of the reconstruction of the original OA sources. This software code makes it possible to simulate a medium for propagation of acoustic waves using such characteristics as the density and ultrasound velocity. The solution of the inversed problem of optoacoustics obtained using the algorithms of the k-Wave Matlab toolbox was corrected using an iterative algorithm based on the Banach fixed point theorem. Four iterative schemes of the correction of a reconstructed image were proposed.
The reconstruction quality was determined using both the quantitative and visual estimates of the obtained results. The efficiency of the iterative algorithm of enhancing the reconstruction quality was quantitatively estimated using the structure similarity index SSIM and relative reconstruction error E. It was shown that this algorithm makes it possible to considerably improve the image quality as compared with the original OA reconstruction.
The results obtained in this work may be important from the point of view of the prospects for their further practical application to solve biomedical imaging problems. The final conclusions about the efficiency of the approach proposed in this work can be drawn after its testing on arrays of practical experimental data.
The work was supported by the Volkswagen Foundation project "Modeling, Analysis and Approximation Theory towards Applications in Tomography and Inverse Problems".