Exterior/interior problem for the circular means transform with applications to intravascular imaging

Exterior inverse problem for the circular means transform (CMT) arises in the intravascular photoacoustic imaging (IVPA), in the intravascular ultrasound imaging (IVUS), as well as in radar and sonar. The reduction of the IPVA to the CMT is quite straightforward. As shown in the paper, in IVUS the circular means can be recovered from measurements by solving a certain Volterra integral equation. Thus, a tomography reconstruction in both modalities requires solving the exterior problem for the CMT. Numerical solution of this problem usually is not attempted due to the presence of"invisible"wavefronts, which results in severe instability of the reconstruction. The novel inversion algorithm proposed in this paper yields a stable partial reconstruction: it reproduces the"visible"part of the image and blurs the"invisible"part. If the image contains little or no invisible wavefronts (as frequently happens in the IVPA and IVUS) the reconstruction is quantitatively accurate. The presented numerical simulations demonstrate the feasibility of tomography-like reconstruction in these modalities.


Introduction
Mathematical models of various imaging modalities are based on the circular means transform (CMT), which maps a compactly supported function of two variables to its integrals along a twodimensional family of circles. Some examples of such imaging techniques include thermo-and photo-acoustic tomography, ultrasound reflection tomography, sonar and radar imaging. Under certain reasonable assumptions, the problem of image reconstruction in all these procedures can be reduced to the inversion of the CMT, i.e. the reconstruction of the unknown image function from its integrals along circles (e.g. see [9,12,13,5,14,20], as well as Section 1 here). The choice of the family of integration circles (and, in particular, their location with respect to the support of the image function) depends on the setup of the imaging device.
Typically, the centers of integration circles correspond to the locations of transducers recording the physical signals (acoustic or electromagnetic waves) coming from the imaging object. As a result, these centers are limited to some discrete locations distributed along a curve, which we call a data acquisition curve. There are abundant results about the inversion of CMT for the so-called interior problem, arising when the data acquisition curve surrounds the support of the image function (e.g. see [10] and the references there). However only limited results are available for the exterior problem, where the support of the image function lies outside of that curve, or the interior/exterior problem with the support partially inside and partially outside of that curve (e.g. see [2] and the references there).
Some of the applications for the interior/exterior problem are radar and sonar imaging ( [13,5]). However, our present work is mostly motivated by the recent advances in intravascular ultrasound imaging (IVUS) and intravascular photoacoustic imaging (IVPA) [7,15,17,16,19]. The measuring device in these modalities is a thin catheter with a set of embedded ultrasound transducers. The catheter is inserted in a blood vessel and records ultrasound waves coming from the surrounding tissue. In IVUS, the transducer first generates an outgoing ultrasound impulse and then switches to the listening mode and records the wave reflected by the tissue. In IVPA the acoustic waves are created by the photoacoustic effect: a short pulse of an infrared laser heats up the tissue that, in turn, expands and generates an outgoing acoustic wave.
To the best of our knowledge, in IVUS and IVPA a tomographic reconstruction currently is not attempted. Instead, a crude image is reconstructed by superimposing the incoming signal with the directional information. One of the goals of this paper is to show that, by utilizing the equipment and measuring techniques similar to those currently used in IVUS and IVPA, one can reconstruct qualitatively (and sometimes even quantitatively) correct images of certain physical properties of biological tissue. As discussed below, in IVPA the image reconstruction problem reduces to the inversion of the CMT by considerations well known in the standard photoacoustic tomography. The situation in IVUS is more complicated: since the reflected wave depends on the speed of sound one tries to recover, the inverse problem is non-linear (unlike that of IVPA). Since in soft tissues variations of the speed of sound are not very large, one can try to linearize the problem using the Born approximation. In Section 2 we use this approach and show that in this case the circular integrals of the speed of sound are related to the measurements by a Volterra integral equation (in time) whose kernel is given by a certain integral expression. We analyze properties of the kernel and show that this integral equation can be stably solved (by means of successive iterations or by linear algebra techniques) to recover the circular means. This reduces the reconstruction problem, again, to the inversion of the CMT.
The exterior inverse problem for the CMT has received so far very little attention in the literature. The main culprit here is the well-known ill-posedness of this problem. In order for a material interface to be "visible" under the CMT, the normal to the interface should intersect the surface supporting the transducers (in our case, the catheter). It is clear that in a generic exterior problem not all interfaces are visible, and therefore accurate and stable reconstruction of material properties is not possible. However, in the intravascular imaging the walls of blood of vessels mostly run in the visible directions, roughly parallel to the surface of the catheter (see, for example, [17] for good images of cross-sections of aorta). As our numerical simulations show (see Section 3), when invisible interfaces are absent, a properly regularized algorithm can stably reconstruct a quantitatively correct image. If the invisible interfaces are present, they will appear blurred and the corresponding part of the image will not be reconstructed correctly. These numerical experiments suggest the feasibility of tomography-like reconstruction in IVPA and IVUS, which can be obtained by the techniques proposed here.
The rest of the paper is organized as follows: In Section 1 we discuss an exact inversion formula for the CMT in the exterior and interior/exterior problem. We then present an algorithm based on that formula for the approximate reconstruction of a function from its circular means. In Section 2 we show that, under Born approximation, the inverse problem for IVUS can be reduced to the inversion of the CMT by solving a Volterra integral equation with a continuous kernel. Section 3 presents results of numerical simulations demonstrating the work of the reconstructions techniques proposed in the previous sections.
1 Exterior/interior problem for the circular means transform 1.1 Formulation of the inverse problem in IVPA It is not difficult to reduce the inverse problem arising in IVPA to the inversion of the CMT. In fact, a very similar problem arises in the photoacoustic tomography with integrating line detectors, and the derivation presented below can be found in the corresponding literature (see, for example [4]).
Let us assume that the acoustic wave is measured by infinitely long transducers placed on the surface of a cylindrical catheter of radius R 0 and parallel to the (vertical) axis of the catheter. We will assume for simplicity that the speed of sound within the catheter coincides with the constant speed of sound in the surrounding tissue. (The constant speed approximation is frequently made in the problems of photo-and thermoacoustic tomography in soft tissue.) Without loss of generality this speed of sound then can be assumed to be 1. The electromagnetic energy of the incoming laser pulse is absorbed by the tissue, which leads to thermoelastic expansion and generation of an acoustic wave. An excess pressure p(t, x) solves the initial value problem for the wave equation in where initial pressure F (x) is determined by the properties of the tissue; this function carries important biological information and is the object of our interest. It is well known that the integrals of p(t, x) along any selected direction solve the 2D wave equation. We, in particular, will consider the integrals u(t, x) in the direction parallel to transducers These integrals satisfy the following 2D initial value problem (IVP) where A transducer passing through the point z = (z 1 , z 2 ) ∈ R 2 (or, correspondingly, through the point (z 1, z 2 , 0) in R 3 ) measures the time series u(t, z) for t ∈ (0, ∞). Since the transducers cover the surface of the catheter, such measurements are done for all z lying on the circle of radius R 0 centered at the origin in R 2 . Ideally, we would like to reconstruct F (x) from the measurements u(t, z). However, since the transducers are infinitely long and integrate the data along straight lines, there is not enough data to reconstruct F (x). Instead, our goal is to reconstruct f (x). Solution of the IVP (1) can be represented with the help of the Green's function G 2D (t, x) of the two-dimensional wave equation in R 2 (describing the outgoing wave) where H(t, s) is the Heaviside function equal to 1 for t > s and equal to 0 otherwise. The well-known Kirchhoff formula [18] yields where the term in the brackets represents the circular integrals g(z, r) of f (x): Taking into account equations (2) and (3) one obtains if t > r, otherwise u(t, z) = 0. Equation (4) is one of the versions of the well-known, explicitly invertible Abel transform (see [4] or [8] Section 4.3); for every transducer (for every z with |z| = R 0 ) function g(z, r) can be reconstructed by the following formula: This solves the problem of finding circular integrals g(z, r) (or, equivalently, circular means g(z,r) 2π ) from the measurements u(t, z).
Next, we develop an algorithm for the approximate reconstruction of a function from its circular means (or circular integrals) in 2D. The centers of the integration circles lie on the circle S 0 of radius R 0 centered at the origin O, as shown in Figure 1. The function f (x) is supported within a larger  Figure 1: Acquisition geometry and "invisible" directions concentric circle of radius R 1 > R 0 . It will be (approximately) reconstructed from the circular integrals g(z, r) with centers z = z(R 0 , ϕ) of the integration circles sampling the whole circle S 0 , i.e. ϕ ∈ [0, 2π]; the radii r of the integration circles cover the interval [0, R 0 + R 1 ]. In other words, we consider an exterior/interior problem for the circular means Radon transform. The exterior problem for this transform is known to be ill-posed. In particular, some wavefronts in the exterior of the circle S 0 are not "visible", i.e. there are no integration circles passing through the corresponding points orthogonally to the wavefront directions. The "invisible" wavefronts cannot be stably reconstructed (e.g. see [20] and the references there). An example of such a wavefront passing through a point p is shown in Figure 1. Here, the shaded area shows the directions of invisible singularities at point P . In particular, if the function has a jump discontinuity across interface I shown in the figure, then it can not be stably recovered from the circular means.
Since an arbitrary function has all possible wavefronts, including the invisible ones, the general exterior (or interior/exterior) problem can only be solved approximately. The algorithm we present below stably reconstructs theoretically visible wavefronts and blurs the invisible ones.

Image reconstruction from the circular means
We start by computing the Hankel transformĝ(z, λ) of g(z, r)/r: The next step is to substitute into (6) the following expression representing the well-known addition theorem (Thm. 2.10 in [6]) for J 0 : x = (r cos θ, r sin θ), z = R 0 (cos ϕ, sin ϕ).
Since the series in (7) converges absolutely and uniformly we obtain Let us expand bothĝ(z(R 0 , ϕ), λ) and f (x(r, θ)) in the Fourier series with respect to the angular variables: Then, by substituting (8) into the second equation in (10) one obtainŝ Let us formally divide equation (11) by J |l| (λR 0 ) (a discussion of this step follows below). We obtain the following equation where The right hand side of (12) is the self-invertible Hankel transform; thus f l (r) can be computed as By combining the equations (5), (10), (13), (14) and (9) one formally reconstructs f (x). This solution was first proposed by Norton [14] for the interior problem, i.e. for the case when the reconstructed function is supported in the interior of S 0 . An important issue not addressed by Norton is the treatment of the zeros of the Bessel functions J |l| (λR 0 ) in the denominator of equation (13). Since the Hankel transform F l of f l (given by (12)) is a bounded function of λ, the ratio (13) remains bounded for all values of λ, implying that the exactly computedĝ l (λ) vanishes at all zeros λ l,m , (m = 0, 1, 2, ...) of the Bessel function J l (λR 0 ). Therefore, all the singularities in (13) are removable, for example, by application of the L'Hospital's rule. This is not a very practical algorithmic solution since an approximation toĝ l (λ) computed from the noisy data will not, in general, have zeros at λ l,m . Even if one continues to formally use the L'Hospital's rule at λ l,m , computing the fractionĝ l (λ)/J |l| (λR 0 ) at values of λ close to λ l,m is still an ill-posed problem. A better solution was proposed in [9] where the Fourier-Bessel series was used in a way that avoids divisions by zeros. A more elegant solution [3,12] (applicable for the interior problem only) is to replace the Bessel function J 0 (λr) in (5) and in (7) by the Hankel function H (1) 0 (λr). This leads to a ratio similar to (13), but with Bessel functions in the denominator replaced by Hankel functions. Since the latter functions do not have zeroes for real values of λ, a straightforward discretization of (14) leads to a stable algorithm.
In the present case of the exterior or exterior/interior problem the methods of [3,9,12] are no longer applicable. In order to avoid divisions by small values of J |l| (λR 0 ) in (13) we will utilize two techniques. The first of them consists in using complex values of λ and deforming the integration contour in (14) to avoid all the zeros of J |l| (λR 0 ) except the first one at λ = 0. In the next section we show that such a deformation does not change the computed value of f l .
In order to deal with the zero at λ = 0 and with the small values of J |l| (λR 0 ) in the neighborhood of this zero we will have to restrict the values of λ and l for which F l (λ) are computed, and to replace by zero the missing values. After such a replacement the algorithm will no longer be theoretically exact, but the computations will be stable. This technique is described in Section 1.4 .

Deforming the contour
In order to avoid division by zeros in (13) we will work with complex values of λ lying on the curve C consisting of the segments [0, ia] and infinite line [ia, ia + M ] with M real, going to ∞. The values ofĝ(z, λ) are computed by formula (5) for all λ ∈ C. Then f l (r) is computed using instead of (14) the following formula Since all zeros of Bessel functions J |l| (t) are real the denominator in the formula (13) (defining F l (λ)) vanishes only at λ = 0 (except the case l = 0 when the denominator does not vanish at λ = 0). Numerical integration in the neighborhood of λ = 0 requires regularization, as described in the next section.
We need to prove that formulas (14) and (15) are indeed equivalent. First, we observe the following well known property: Proposition 1 For each l ∈ Z, function F l (λ) given by equation (12) is an entire function.
Proof. This fact follows from the well-known relation between the Hankel transforms (12) and the entire 2D Fourier transform of f (x). I.e., up to a constant factor, F l (λ) is the restriction of the entire 2D Fourier transformf l (ξ 1 , ξ 2 ) of a finitely supported function f l (r) exp(ilθ) to the complex plane ξ 1 = −λ ∈ C, ξ 2 = 0. Now, in order to prove that the integral over the positive real values of λ can be replaced by the integral over C it is enough to show that Let us first recall the asymptotic behavior of the Bessel function J ν (z) for large values of argument z (see [1], formula 9.2.1) It follows that for λ lying on the segment Now, in order to prove (16) it is enough to show that F l (λ) decays sufficiently fast in the limit b → ∞. This condition is given by the following Proof. In addition to the asymptotic decay for large values of the argument mentioned above, Bessel function J n (z) is bounded for any integer index n and any z ∈ D. This can be easily deduced, for example, from the well-known representation of the Bessel function ([1], formula 9.1.21) This fact can be used to derive an estimate on F l (λ). Taking into account finite support of f l (r), from the definition of F l (λ) we obtain Since both f l (r) and J |l| (λr) are bounded, the first integral in the right hand side decays as O(b −1 ): The last integral in (18) over the interval [1/ √ b, R 1 ] can be estimated by using the large argument asymptotics (17), since |λr| ≥ √ b on this interval: Next, by adding another small term (uniformly decaying as O 1 b ), the integral in (20) can be extended to the whole interval [0, R 1 ]: Suppose that λ = b + iα, with α ∈ [0, a]. For a fixed α the last integral in (21) can be re-written as a combination of a sine and cosine Fourier transforms of a finitely supported integrable functions cosh(αr) √ rf l (r) and sinh(αr) √ rf l (r): The well-known Riemann-Lebesgue lemma then implies that for a fixed α We cannot claim yet that this decay is uniform with respect to Im λ. However, (22) holds, in particular for α = 0 and α = a, which correspond to the lower and upper boundaries ∂D lower and ∂D upper of the strip D. Now, by combining estimates (18)-(22), one obtains for every fixed α ∈ [0, a]. It follows that the real and imaginary parts of the function F l (λ) are harmonic functions within the strip D, vanishing at infinity. The values of these functions on ∂D lower and ∂D upper decrease as o b −1/2 when b → ∞. Now the uniform decay within D can be established by applying the proof presented in the Appendix of [11], resulting in the estimate Now vanishing of the integral in (16) is guaranteed by combining (24) with (17) and we thus proved the following Theorem 3 Integration over the positive real axis in (14) is equivalent to integration over contour C in (15).

Regularization
The contour deformation described in the previous section eliminates the division by small values in the denominator of (13) for all values of λ, except λ = 0. Unlike other zeros, at λ = 0 Bessel function J ν (λ) has a root of multiplicity |ν| (we consider only integer ν). In addition J ν (t) remains small (much less than one) in the interval t ∈ [0, ν(1 − ε))] (with the value of ε depending on the desired notion of "smallness"). A more quantitative description of this decay is given by the first term of the Debye asymptotics ([1], formula 9.3.7): It follows that in any direction x = kν, with |k| < 1, values of J ν (x) decay exponentially as x → ∞. Such behavior of Bessel functions results in extremely small values of F l (λ) andĝ l (λ) for certain values of λ. Since f l (r) are finitely supported in the interval [0, R 1 ], equation (12) shows that F l (λ) become very small for values of |l| ≥ λR 1 . For simplicity we neglected here the narrow transition zone. We also notice that this effect is preserved when λ is slightly shifted into the imaginary direction (|l| ≫ 1 ≥ | Im λ|). In other words, the vanishing values lie within the cone |l| ≥ R 1 Re λ.
We notice, parenthetically, that this effect is closely related to the well-known "bow-tie" shape of the Fourier spectrum of projections in classical X-ray tomography. Functions F l (λ) can be understood as the Fourier coefficients obtained by expanding the Fourier transformf (ξ) of f (x) represented in polar coordinates in the 1D Fourier series in the angular coordinate. It follows that if f (x) is supported within a circle of radius R 1 then there is little energy in the cone |l| ≥ R 1 λ; these values do not need to be computed and can be set to zero without much effect on the reconstructed image.
However, functionsĝ l (λ) given by (11) have an additional Bessel factor J |l| (λR 0 ) in front of the integral. According to the asymptotic behavior of the Bessel functions, the exact values of g l (λ) will become very small in the larger cone |l| ≥ R 0 Re λ. Since we need to compute (13) from the approximately known values ofĝ l (λ), the division within the cone |l| ≥ R 0 Re λ is an ill-posed operation and should be avoided. The regularization step of the algorithm consists in replacing values of F l (λ) by zeros within the cone |l| ≥ R 0 Re λ. While values of F l (λ) within the smaller region of |l| ≥ R 1 Re λ are very small and can be set to zero without noticeably affecting the image, the additional loss of values in the region R 0 Re λ ≤ |l| ≤ R 1 Re λ will cause the disappearance of certain material interfaces (or wavefronts).
It is interesting to investigate which wavefronts will remain in the image and which will be smoothed out. Notice that the reconstructed image is obtained by combining equations (9) and (15), in other words Setting to zero certain values of F l (λ) is equivalent to removing from the reconstructed image components in the form J |l| (λr)e ilθ for certain values of λ (for simplicity we will ignore in the following crude analysis the small imaginary part of λ).
Let us analyze local behavior of a function J |l| (λr)e ilθ in a small neighborhood N of point p as shown in Figure 1. (Due to the rotational symmetry of the acquisition scheme it is enough to understand this behavior only for the points located on the vertical axis). Within N we have x 1 ≈ Rθ, so that the angular exponent e ilθ can be approximated by exp(ilx 1 /R), where x = (x 1 , x 2 ). The function u(x(r, θ)) ≡ J |l| (λr)e ilθ solves the Helmholtz equation ∆u + λ 2 u = 0.
Locally the dependence of this function on x 1 can be approximated by exp(ilx 1 /R). Then u(x) can be approximated by the formula with some constants c 1 and c 2 , and with the vertical frequency γ satisfying the condition The two plane waves given by (25) propagate in the directions given by the wave vectors v + and v − : v ± = (l/R, ±γ) = λ(cos θ, sin(±θ)), where the angle θ between v + and the horizontal axis is given by the condition cos θ = l λR .
Now, if components J |l| (λr)e ilθ with λR 0 ≤ |l| are filtered out by the regularization step of the algorithm, the boundary of the cone of excluded directions are given by the condition However, this condition coincides with the visibility condition discussed in the beginning of Section 1. It follows that our regularization technique removes from the image plane waves propagating in the "invisible" directions but keeps all the others.

Inverse problem for IVUS and the CMT
While the inverse problem of IVPA is reduced to the inverse problems for the CMT in a rather straightforward way (see Section 1.1), the situation in IVUS is more complicated. The important difference between these modalities lies in the duration of the excitation. In IVPA the acoustic wave is excited by a short laser pulse; it's duration is much shorter than the time needed for an acoustic wave to cover a distance compared to the desired resolution of the image. Thus, the forward problem can be modelled by the wave equation without sources (the source term is absorbed in the initial condition at t = 0). In the IVUS, on the other hand, the reflected acoustic wave that carries the desired information is generated by the incoming excitation wave initiated by the transducer; both waves propagate with the same speed. In addition, while in the IVPA the speed of sound can be assumed constant and known, in the IVUS the variations in the speed of sound is the information one wants to reconstruct. This makes the inverse problem non-linear; we will linearize it by using the Born approximation. In this section we show that under the latter approximation the circular means of the speed of sound (centered at each transducer) can be reconstructed from the measurements by solving a certain Volterra integral equation. After the circular means are found, the exterior problem for the CMT can be solved by the methods presented in the previous section.
We will assume that each transducer is infinitely long and the measurements are made by each transducer sequentially and independently from the others. A transducer initiates an outgoing excitation wave u exc (t, x) and then switches into the recording mode. The transducers are presumed to be infinitely thin, so that after the initial pule is generated, they do not perturb the propagation of acoustic waves. Under these assumptions function u exc (t, x) satisfies the free-space wave equation We will assume that the speed of sound c(x) is close to a constant and this constant (by choosing proper physical units) can be made equal to unity. In other words Now (26) can be re-written in the following form: Let us also consider the solution u 0 (t, x) of the wave equation in a homogeneous media (excited by the same transducer). It satisfies the homogeneous wave equation Now the difference w(t, x) = u exc (t, x) − u 0 (t, x) satisfies the equation The Born approximation results from neglecting the last term in the above equation (this is a second order term with respect to m(x) since ∆w has the same order as m(x)). Taking into account (29) we obtain Our goal is to reconstruct from the measurements some information about m(x). Since the transducers we model do not have any resolution in the vertical direction, we will only be able to partially reconstruct the integrals m(x 1 , x 2 ) of m(x) in x 3 : The first step towards this goal is to reconstruct from the measurements obtained by one transducer (without loss of generality assumed to be located along the x 3 axis) the circular averages M (r) of m(x 1 , x 2 ) defined as follows M (r) = 2π 0 m(r cos θ, r sin θ)dθ.
If the wave u 0 (t, x) is excited by an infinitely thin transducer lying on the Ox 3 axis, u 0 (t, x) should be invariant with respect to x 3 , and invariant with the respect to rotations about the axis Ox 3 . We will represent u 0 with the help of the Green's function G 2D (t, x 1 , x 2 ) of the two-dimensional wave equation in R 2 satisfying the radiation condition at infinity and given by equations (2) and (3): where r(x) = x 2 1 + x 2 2 , and ϕ(τ ) is a C ∞ function in R finitely supported within the interval [−a, 0]. The function ϕ(τ ) is a delta-approximating function describing the initial pressure on the surface of the transducer. It is introduced in order to avoid some technical difficulties; further in this section we will pass to the limit a → 0 and ϕ(τ ) → δ(τ ). We notice that for such a choice of ϕ(τ ), the function u 0 (t, x) still solves the wave equation in the whole space R 2 for all t ≥ 0, and this solution is invariant with respect to x 3 (it depends only on r(x) and t).
Consider the free space Green's function G 3D (t, x) of the 3D wave equation satisfying the radiation condition at infinity Using G 3D (t, x), solution of equation (30) can be re-written in the form Under an additional assumption that m(x) is finitely supported in space in x 3 (or that it decreases at infinity sufficiently fast), one can consider integralsw(t, The transducer measuresw(s, 0, 0), i.e. the integral of w(s, x) over the Ox 3 axis: where we interchanged the order of integrations and made use of the fact that G 2D (t, −y 1 , −y 2 ) = G 2D (t, y 1 , y 2 ) = R G 3D (t, (y 1 , y 2 , y 3 ))dy 3 .
By using (31) equation (33) can be further simplified as follows By utilizing formula (2) and by integrating (34) in polar coordinates,w(s, 0, 0) can be expressed in terms of the circular averages M (r) as follows The substitution of (32) into (35) results in the following integro-differential equation relating circular averages M (r) with the measurementsw(s, 0, 0): Below we show that by a proper choice of ϕ(τ ) the above equation can be reduced to the Volterra integral equation of the second kind. Let us consider the anti-derivative Ψ (t, r) of Φ (t, r) in t: Then, by taking into account the finite support of ϕ(τ ), the expression in the brackets in (36) can be transformed as follows: Now (36) can be re-written in the following form Let us consider the case of ϕ(τ ) equal to the Dirac's delta function δ(τ ). To this end introduce a family of delta-approximating functions ϕ α (τ ) and the corresponding family of measurements w α (s, 0, 0) so that lim α→0 ϕ α (τ ) = δ(τ ).
Then, taking the limit α → 0 yields where the second derivative of K(r, s) should be understood in the sense of distributions. In fact, it will be more convenient for us to work with the anti-derivative of the data W δ (s) ≡ w δ (s, 0, 0)ds. The latter function is related to the averages M (r) by the equation where the derivative is, again, understood in the sense of distributions. Let us investigate the behavior of K(r, s). By combining (3) and (37) we obtain We observe that if s ≤ 2r, the numerator of the integrand identically vanishes and K(r, s) = 0. Otherwise, if s > 2r, By substitution q = t − s/2 this can be further simplified to Integration by parts, followed by substitution η = q/(s/2 − r) further yields We notice that for r > 0 the expression in the denominator of the last integrand is bounded away from zero: Thus, as s approaches 2r, the second term in (39) vanishes as O (s/2 − r) 2 . Moreover, it is easy to check that the derivative (in s) of this term also vanishes as s approaches 2r. Therefore, as a function of s, K(r, s) has a jump at s = 2r which can be represented using the Heaviside function in the form π 2r H(s, 2r) with the remaining part continuous through the point s = 2r: H(s, 2r) + K 1 (r, s), 2r H(s, 2r), s > 2r. Now, the derivative of K(r, s) has form where the second term is a bounded (although discontinuous at s = 2r) function. The second term can be extended by continuity to the point s = 2r. Substituting this expression in (38) yields This is a Volterra integral equation of the second kind, with a continuous kernel. Equations of this type are well-posed and have unique solutions. They are also easy to solve numerically. For example, the method of successive approximations will converge unconditionally; it consists in computing the next approximation M (k+1) from the previous one M (k) by the formula Alternatively, our analysis of the kernel ∂ ∂s K(r, s) suggests that, when properly discretized, equation (38) will reduce to a well-posed system of linear equations, with a triangular matrix. We chose this alternative to conduct numerical simulations described in the next section.

Numerical realization and simulations
We present below the results of numerical simulations illustrating the work of the algorithms proposed in the previous sections. The first two series of simulations test the performance of the CMT inversion technique developed in Section 1.

Inversion of the CMT
We first apply the CMT inversion algorithm to a standard interior problem. This problem is wellposed and a good quantitatively correct reconstruction is expected. As a phantom, we used a set of functions whose gray-scale image is shown in Figure 2(a). These functions are mostly constant (equal 1) within their support, however the transition from 1 to 0 is smooth. This smoothness is especially important in the next section where the wave equation was solved by finite difference methods whose accuracy would be severely compromised if the test phantom were discontinuous. The simulated transducers are passing through the circle of radius 0.5 perpendicular to the plane of the Figure; their locations are shown by the dotted line in Figure 2(a). Since the phantom is completely surrounded by the transducers, this is an interior problem. There were 256 simulated transducers in this experiment; for each of them 401 circular means where computed, with the radii ranging from 0 to 2. Figure 2(b) shows reconstruction of the function from the circular means by the method presented in Section 1. As expected, the reconstruction is quite accurate. In order to test the performance of our method in the exterior/interior problem we added to the phantom an additional function with the intention to crudely model aorta walls, as seen, for example, in [17]. The material interfaces modelled by this function are "visible". The reconstruction from the circular means is shown in Figure 2(d). One can see that the interior part of the phantom is reconstructed as well as before. In spite of the ill-posedness of the exterior problem and the inexact nature of the regularized algorithm, the exterior part of the phantom is also reconstructed quantitatively correct, although a careful reader will notice some variations in the brightness of the "aorta walls". It is interesting to see what effect on the reconstruction will have the presence of "invisible" interfaces. To this end we added to the phantom several circular inclusions (see Figure 3(a)). Figure 3(b) demonstrates the reconstruction from the accurate circular means. In Figure 3(c) before the reconstruction a simulated noise was added to the circular means; the intensity of the noise was 5% of the "signal" in L 2 norm. For a fair comparison images in Figure 3 are drawn using the same gray scale. Images in Figures 3(b) and 3(c) show that, as expected, the "invisible" material interfaces are significantly blurred, and the details (circular inclusions) are not reconstructed correctly. Nevertheless, the visible parts of the corresponding boundaries are captured, and the artifacts are relatively well localized, so that "aorta walls" are still clearly seen.

Simulation of the IVUS
Our next simulation tests the feasibility of tomography reconstruction in IVUS. The geometry of the simulation is similar to that of the previous section. This time, however, for every transducer location using the time-stepping finite-difference technique we solved two wave equations. One of them corresponds to the uniform speed of sound; it models u 0 (t, x) from Section 2. Since the initial condition does not depend on x 3 , this is actually a 2D problem. The other equation we solve is equation (26) with a variable speed of sound, again not depending on x 3 .The simulated transducers "measured" the difference w(t, x) of these two solutions.
The variable speed of sound in the first of our simulations was equal to 1 plus 1% perturbation modelled by the phantom from the previous section. The perturbation is shown in Figure 4(a). The circular means of the perturbation were reconstructed using the method from Section 2, and the perturbation was then reconstructed from the circular means by the method of Section 1. The result is shown in Figure 4(b). Although additional artifacts can be noticed in the image (as compared to Figure 3(b)), the "aorta walls" are still clearly seen and the interior part of the image is reconstructed well. The additional artifacts are caused by the imprecise nature of the method (the Born approximation is used) and by the errors introduced by finite differences when modelling the forward problem.
In order to better understand the effect of the Born approximation we repeated the last simulation with new variable speed of sound equal 1 plus a 5% perturbation given by the same function. The result is shown in Figure 4(c) (with the gray scale modified by the factor of 5). One can see additional artifacts in the central part of the image; they are clearly caused by the increased error in the Born approximation, since everything else remains the same. Nevertheless the rest of the image is reconstructed quite well, and the parts of the phantom with visible material interfaces are reconstructed with a reasonable accuracy.

Concluding remarks
We have considered in the present paper the inverse problems that arise in IVPA and IVUS, and have shown that under several simplifying assumptions these problems can be reduced to the solution of the exterior problem for the CMT. The latter is, in general, severely ill-posed. However, the regularized algorithm proposed in Section 1 yields stable reconstruction of the visible material interfaces while it blurs the invisible ones. If the invisible interfaces are absent or almost absent in the image, the algorithm reconstructs a quantitatively correct image. Such a situation (absence of invisible interfaces) is not uncommon in practical application of IVPA and IVUS, judging from the images we found in the literature.
In order to carry out our analysis we have made several simplifying assumptions. In particular, the transducers were assumed infinitely long and infinitely thin, the body of the catheter was assumed to be made of material with the same speed of sound as that in blood. We also assumed that the speed of sound in the soft tissues constituting and surrounding the vessels is close to the speed of sound in blood. The latter assumption is standard in problems of thermoacoustic tomography; it allowed us to use the constant speed approximation in IVPA and the Born approximation in IVUS. The model of acoustically transparent and infinitely long catheter is more crude; it was used mostly to simplify the analysis and to develop practically useful reconstruction algorithms. It is not clear whether it is feasible or practical to manufacture the catheter whose speed of sound matches that of blood or water. The effect of the finite length of the transducers also requires further investigation.
However, our purpose was not to develop a perfect reconstruction algorithm, but rather to demonstrate the feasibility of tomography-like reconstruction in IVPA and IVUS. More sophisticated techniques, that would take into account the details we neglected or over-simplified in the present, first approach to the problem, will be the subject of the future work.