Resolution analysis of inverting the generalized Radon transform from discrete data in $\mathbb R^3$

TL;DR

This paper analyzes the spatial resolution and accuracy of reconstructing singularities of functions in 3D from discretized generalized Radon transform data, providing explicit singular behavior and conditions for artifact-free recovery.

Contribution

It offers a detailed analysis of the resolution limits and singularity reconstruction accuracy for the GRT in 3D from discrete data, including explicit formulas and conditions for artifact suppression.

Findings

Explicit leading singular behavior of reconstructed functions near discontinuities.

Conditions under which singularities do not produce non-local artifacts.

Numerical experiments confirm theoretical predictions.

Abstract

A number of practically important imaging problems involve inverting the generalized Radon transform (GRT) $\mathcal R$ of a function $f$ in $\mathbb R^3$. On the other hand, not much is known about the spatial resolution of the reconstruction from discretized data. In this paper we study how accurately and with what resolution the singularities of $f$ are reconstructed. The GRT integrates over a fairly general family of surfaces $\mathcal S_y$ in $\mathbb R^3$. Here $y$ is the parameter in the data space, which runs over an open set $\mathcal V\subset\mathbb R^3$. Assume that the data $g(y)=(\mathcal R f)(y)$ are known on a regular grid $y_j$ with step-sizes $O(\epsilon)$ along each axis, and suppose $\mathcal S=\text{singsupp}(f)$ is a piecewise smooth surface. Let $f_\epsilon$ denote the result of reconstruction from the descrete data. We obtain explicitly the leading singular behavior of $f_\epsilon$ in an $O(\epsilon)$-neighborhood of a generic point $x_0\in\mathcal S$, where $f$ has a jump discontinuity. We also prove that under some generic conditions on $\mathcal S$ (which include, e.g. a restriction on the order of tangency of $\mathcal S_y$ and $\mathcal S$), the singularities of $f$ do not lead to non-local artifacts. For both computations, a connection with the uniform distribution theory turns out to be important. Finally, we present a numerical experiment, which demonstrates a good match between the theoretically predicted behavior and actual reconstruction.

Equations

$$\mathcal{C}:=\{(x,y)\in \mathcal{U}\times \mathcal{V}:\Phi (x,y)=0\}.$$

$${\mathcal{T}}_x:=\{y\in \mathcal{V}:\Phi (x,y)=0\},x\in \mathcal{U}.$$

$${\Phi }_x^{\prime }(x,y)\ne 0,{\Phi }_y^{\prime }(x,y)\ne 0,(x,y)\in \mathcal{C};$$

$$\text{det}\left(\frac{{\partial }^2\Phi (x,y)}{\partial x_i\partial y_j}\right)\ne 0,(x,y)\in \mathcal{C},$$

$$\text{det}\left(\begin{array}{cc}{\Phi }_{xy}^{\prime \prime }& ({\Phi }_x^{\prime }{)}^T\\ {\Phi }_y^{\prime }& 0\end{array}\right)\ne 0.$$

$$f(x)=\sum _j{\chi }_{D_j}{f}_j(x),$$

$$g(y)=(\mathcal{R}f)(y):={\int }_{S_y}b(x,y)f(x)dx,y\in \mathcal{V},$$

$$g(ϵj),j\in r+{\mathbb{Z}}^3,$$

$$C:=\{((x,\lambda {\Phi }_x^{\prime }(x,y)),(y,-\lambda {\Phi }_y^{\prime }(x,y)):\Phi (x,y)=0,\lambda \in \mathbb{R}\setminus 0,x\in \mathcal{U},y\in \mathcal{V}\},$$

$$N(x):={\text{II}}_{{\mathcal{S}}_y}(x)-{\text{II}}_{\mathcal{S}}(x),$$

$${\Gamma }_x:=\{y\in \mathcal{V}:x\in {\mathcal{S}}_y,{\mathcal{S}}_y\text{ is tangent to S at some z},z\ne x\}.$$

$$\sum _{j\in {\mathbb{Z}}^3}{j}^m\phi (u-j)=u^m,m\in (0\cup \mathbb{N}{)}^3,|m|\le 2,u\in {\mathbb{R}}^3;$$

$$g_{ϵ}(y):=\sum _{j\in r+{\mathbb{Z}}^3}g(ϵj)\phi \left(\frac{y-ϵj}{ϵ}\right).$$

$$\begin{array}{rl}{f}_{\chi }(x)& =-\frac{1}{4{\pi }^2}{\int }_{S_{+}^2}\frac{\chi (Y(\alpha ,0;x))}{b(x,Y(\alpha ,0;x))}{{\left(\frac{\partial }{\partial t}\right)}^{2}g(Y(\alpha ,t;x))|}_{t=0}d\alpha .\end{array}$$

$$\begin{array}{r}x=x_1{\alpha }_0+x^{\perp },x^{\perp }\in {\alpha }_0^{\perp },\end{array}$$

$$x_{ϵ}:=x_0+ϵ\tilde{x},\tilde{x}\in \tilde{\mathcal{U}},$$

$$\begin{array}{rl}{f}_{\chi }(x_{0^{\pm }})& :=\underset{ϵ\to 0^{\pm }}{\lim }{f}_{\chi }(x_0+ϵ{\alpha }_0),f_0:=\underset{ϵ\to 0^{+}}{\lim }(f(x_0+ϵ{\alpha }_0)-f(x_0-ϵ{\alpha }_0)),\\ \nu & :=\frac{|{\Phi }_x^{\prime }|}{|{\Phi }_y^{\prime }|},{\beta }_0=\frac{{\Phi }_y^{\prime }}{|{\Phi }_y^{\prime }|}.\end{array}$$

$$y=(y_1,y^{\perp })=y_1{\beta }_0+y^{\perp }.$$

$$\hat{f}(u,s):=\int f(x)\delta (u\cdot x-s)dx,u\in {\mathbb{R}}^3\setminus 0,$$

$$\underset{ϵ\to 0}{\lim }{f}_{\chi ϵ}(x_{ϵ})=f_{\chi }(x_{0^{+}})-f_0{\int }_{\nu h}^{\infty }\hat{\phi }({\beta }_0,s)ds,$$

$$\underset{ϵ\to 0}{\lim }{f}_{\chi ϵ}(x_{ϵ})=f_{\chi }(x_{0^{-}})+f_0{\int }_{-\infty }^{\nu h}\hat{\phi }({\beta }_0,s)ds.$$

$$\underset{ϵ\to 0}{\lim }{f}_{\chi ϵ}(x_{ϵ})=f_{\chi }(x_{0^{+}})-f_0{\int }_{\nu h}^{\infty }\hat{\phi }({\beta }_0,s)ds,$$

$$x=z+pn(z),z\in \mathcal{S},n(z)\text{ is normal to S at z},|n(z)|\equiv 1.$$

$${\Phi }_x^{\prime }(z+pn(z),y)-\lambda n(z)=0,\Phi (z+pn(z),y)=0,\Psi (z)=0.$$

$$P_y^{\prime }(y_0)=\frac{1}{|{\Phi }_x^{\prime }|}{\Phi }_y^{\prime }\ne 0;$$

$$\left[\begin{array}{cccc}{\Phi }_{xx}^{\prime \prime }-\lambda n_x^{\prime }& {\Phi }_{xx}^{\prime \prime }{\alpha }_0& -{\alpha }_0& {\Phi }_{xy}^{\prime \prime }\\ {\Phi }_x^{\prime }& {\Phi }_x^{\prime }\cdot {\alpha }_0& 0& {\Phi }_y^{\prime }\\ {\Psi }_x^{\prime }& 0& 0& 0\end{array}\right].$$

$${\Phi }_{x^{\perp }{x}^{\perp }}^{\prime \prime }-\lambda (n^{\perp }{)}_{x^{\perp }}^{\prime }.$$

$${\Phi }_{x^{\perp }{x}^{\perp }}^{\prime \prime }-\lambda (n^{\perp }{)}_{x^{\perp }}^{\prime }=-\lambda N(x_0).$$

$$g(\overline{Y}(z,p))=p_{+}G(z,p)+G_1(z,p)\text{ and }G(z,0)=f_0(z)b(z,\overline{Y}(z,0))\frac{2\pi }{\sqrt{\text{det}N(z)}},$$

$$\begin{array}{rl}& g(y)={\int }_{{\mathcal{S}}_y}f(x)b(x,y)\theta \left(\Psi (x)\right)dx,\end{array}$$

Full text

1Department of Mathematics, University of Central Florida, Orlando, FL 32816.

This work was supported in part by NSF grant DMS-1615124.

1. Introduction

A large number of practically important imaging problems involve inversion of the generalized Radon transform (GRT), i.e. recovering an unknown function $f$ from its integrals over a family of surfaces. The reconstruction may involve finding $f$ itself, or finding $f$ modulo smoother terms. Most of the times, the surfaces are not planes. Below is a list of some of the most common integral transforms with some of the most prominent examples of their use.

1. (1)

Integration over spheres. Applications include ultrasound imaging (or, SONAR) (see [15] and references therein) as well as thermoacoustic and photoacoustic tomography [12, 21]. 2. (2)

Integration over ellipses. This transform arises in linearized seismic imaging with a common offset between the sources and receivers [7]. 3. (3)

Integration over cones arises in Compton camera imaging. Applications are single-scattering optical tomography, Compton camera medical imaging, and homeland security (see [19] for a recent review).

In all of the above cases, one collects a discrete data set and reconstructs $f$ using a numerical algorithm. Frequently, reconstruction is achieved by applying a linear inversion formula (as opposed to a non-linear reconstruction algorithm based on fidelity functional minimization). In all of the above examples it is of fundamental importance to know the resolution of the method as a function of (1) the data sampling rate, and (2) specific implementation of the inversion formula that is used. Despite the significance of this problem, not much is known about the resolution of reconstruction from discrete data. The main reason for this is that the classical sampling theory, which addresses such problems, can be applied only in a few simplest settings of the classical Radon transform (CRT) [14]. The known results are quite scarce, and they are of semi-qualitative nature (see e.g. pp. 784–786 in [6]). Very recently, a more flexible approach to sampling based on semiclassical analysis was proposed in [18]. Let $A$ be a Fourier Integral Operator (FIO). The idea of [18] is to determine how the data $Af$ should be sampled to allow for accurate interpolation of its values on a lattice provided that $f$ is semiclassically bandlimited. If the sampling condition is violated, then reconstruction from the discrete values of $Af$ (i.e., applying a parametrix $A^{-1}$ to the interpolated $Af$) leads to aliasing artifacts, which are also analyzed in [18].

An alternative approach to the analysis of resolution was proposed recently in [10, 11]. The idea is to investigate how accurately and with what resolution the singularities of $f$ are reconstructed. For some of the above problems there is no exact inversion formula, and inversion modulo smoother terms is the most one can hope for. In such cases, spatial resolution of the recovery of singularities is all one needs. Note that in this paper both $f$ and $g=\mathcal{R}f$ are assumed to have singularities in the sense of a conventional, classical wavefront set (see e.g. [8]). In contrast, the main assumption in [18] is that $f$ and, consequently, the data $Af$ have only semiclassical singularities (see e.g. [23]). It is possible to apply the approach of [18] to the analysis of classical singularities, but this would require summing a series over “folded” frequencies in the Fourier domain, which is complicated.

In [10, 11] the author considers the inversion of the CRT of $f$ in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$. The parametrization of the data is standard, i.e. in terms of the affine and angular variables. Suppose the step-sizes along the angular and affine variables are $O(\epsilon)$. Let $f_{\epsilon}$ denote the result of reconstruction from the discrete data. The author picks a point $x_{0}$, where $f$ has a jump singularity, and obtains explicitly the leading singular behavior of $f_{\epsilon}$ in an $O(\epsilon)$-neighborhood of $x_{0}$ as $\epsilon\to 0$. The obtained behavior, which we call edge response, provides the desired resolution of the reconstruction algorithm. It is shown also that convex parts of the singular support of $f$ do not create non-local artifacts. The case when $f$ changes during the scan (so-called, dynamic CT) is considered in the 2D setting as well [10].

In this paper we generalize the approach of [10, 11]. The reconstruction problem is now formulated in terms of the GRT $\mathcal{R}$, which integrates $f$ over a fairly general family of surfaces $\mathcal{S}_{y}$ in $\mathbb{R}^{3}$. Here $\text{supp}(f)\subset\mathcal{U}$, where $\mathcal{U}\subset\mathbb{R}^{3}$ is an open set, and $y$ is the parameter in the data space. For the problem to be well-determined, we assume that $y$ runs over an open set $\mathcal{V}\subset\mathbb{R}^{3}$. As is seen, our setting is fairly general and covers all the problems mentioned above. The GRT in this paper is very close to that considered by Beylkin in [2], only the parametrization of the surfaces $\mathcal{S}_{y}$ is slightly different. This gives us more flexibility to connect our results with practical applications, where GRTs arise.

Assume that the data $g=\mathcal{R}f$ are known on a regular grid $y_{j}$ with step-sizes $O(\epsilon)$ along each axis. Suppose $\mathcal{S}=\text{singsupp}(f)$ is a piecewise smooth surface. Similarly to [10, 11], we obtain explicitly the leading singular behavior of $f_{\epsilon}$ in an $O(\epsilon)$-neighborhood of a generic point $x_{0}\in\mathcal{S}$, where $f$ has a jump discontinuity. We also prove that under some generic conditions on $\mathcal{S}$ (which include, e.g. a restriction on the order of tangency of $\mathcal{S}_{y}$ and $\mathcal{S}$), the singularities of $f$ do not lead to non-local artifacts. For both computations, a connection with the uniform distribution theory [13] turns out to be important. It is possible that violation of the imposed conditions leads to artifacts. Analysis of such artifacts and analysis of more general surfaces $\mathcal{S}$ will be the subject of future research.

The reconstruction formula $g\to\check{f}$, which contains a suitably adapted adjoint $\mathcal{R}^{*}$, is one specific example of an FIO. Here $\check{f}$ is such that $\check{f}-f$ is smoother than $f$. Thus, the reconstruction algorithm can be viewed as an application of an FIO to discrete data $(\mathcal{R}f)(y_{j})$. A number of methods for computing the action of FIOs on discrete data have been proposed, see e.g. [4, 5, 1, 22] and references therein. To the best of the author’s knowledge, proposed is the first method to compute the resolution of the reconstruction obtained by applying an FIO to discrete data that comes from an image with classical singularities. Extension of the method to more general FIOs will also be the subject of future work.

The paper is organized as follows. In Section 2 we define the GRT $\mathcal{R}$ via an incidence relation $\mathcal{C}\subset\mathcal{U}\times\mathcal{V}$, list the properties of the function $\Phi(x,y)$ that defines the incidence relation, define generic points, and specify the continuous and discrete inversion formulas that are used in the analysis. The main result is formulated in Section 3, which also contains the beginning of the proof. The entire proof spans Sections 3–7. In Section 3 we obtain the behavior of $g$ near its singular support, which generalizes one of the results of [16, 17] from the CRT to the GRT. The behavior of the interpolated data $g_{\epsilon}$ near $\text{singsupp}(g)$ is obtained in Section 4. The contribution of the leading singular term to the edge response at a generic point $x_{0}\in\mathcal{S}$ is computed in Section 5. In Section 6 we show that lower order terms do not contribute to the edge reponse. In Section 7 we prove that, under some assumptions, remote singularities do not contribute to the edge response as well. Results of a numerical experiment, which show a good match between the theoretically predicted behavior and actual reconstruction, are in Section 8.

2. Preliminary construction

Let $\mathcal{U},\mathcal{V}\subset\mathbb{R}^{3}$ be two open connected sets, where $\mathcal{U}$ is the image domain, and $\mathcal{V}$ is the data domain. Each $y\in\mathcal{V}$ determines a smooth surface $\mathcal{S}_{y}\subset\mathcal{U}$. Let $\mathcal{C}$ be the corresponding incidence relation $\mathcal{C}\in\mathcal{U}\times\mathcal{V}$, which is defined in terms of a smooth function $\Phi(x,y)\in C^{\infty}(\mathcal{U}\times\mathcal{V})$:

[TABLE]

Another way to state (2.1) is that $x\in\mathcal{S}_{y}$ if and only if $\Phi(x,y)=0$. Define the submanifold:

[TABLE]

Thus, $\mathcal{T}_{x}$ is the collection of all $y\in\mathcal{V}$ such that $\mathcal{S}_{y}$ contains $x$. The main assumptions about $\Phi$ are as follows (“DF” stands for Defining Function):

• DF1.

$\Phi$ is real-valued and non-degenerate, i.e.

[TABLE]

• DF2.

For each $x\in\mathcal{U}$, the map $\mathcal{T}_{x}\to S^{2}$ defined by $y\to\pm\Phi_{x}^{\prime}(x,y)/|\Phi_{x}^{\prime}(x,y)|$, $y\in\mathcal{T}_{x}$, is surjective;

• DF3.

For each $y\in\mathcal{V}$, the vectors $\Phi^{\prime}_{y}(x,y)$ and $\Phi^{\prime}_{y}(z,y)$ are not parallel whenever $x,z\in\mathcal{S}_{y}$, $x\not=z$; and

• DF4.

The mixed Hessian of $\Phi$ is non-degenerate,

[TABLE]

where $\partial/\partial x_{i}$, $i=1,2$, and $\partial/\partial y_{j}$, $j=1,2$, are basis vectors in the tangent spaces to the submanifolds $\mathcal{S}_{y_{0}}$ and $\mathcal{T}_{x_{0}}$ at $x_{0}$ and $y_{0}$, respectively.

Condition DF2 means that the tomographic data are complete, i.e. any singularity is visible. Condition DF3 says that there are no conjugate points. Condition DF4 is a local version of the Bolker condition. Conditions DF3 and DF4 imply the (global) Bolker condition. Conditions DF1–DF4 are analogous to Conditions (I)-(IV) in [2]. Conditions DF1 and DF4 combined are equivalent to the condition (cf. eq. (4.23), [20], p. 335) that at every point $(x,y)\in\mathcal{C}$:

[TABLE]

In the paper we consider functions, which can be represented as a finite sum

[TABLE]

where $\chi_{D_{j}}$ is the characteristic function of the domain $D_{j}\subset\mathcal{U}$. For each $j$:

1. (1)

$D_{j}$ is bounded, 2. (2)

The boundary of $D_{j}$ is piecewise $C^{\infty}$, 3. (3)

$f_{j}$ is $C^{\infty}$ in a domain containing the closure of $D_{j}$.

Denote $\mathcal{S}:=\cup_{j}\partial D_{j}$. By construction, $\text{singsupp}(f)\subset\mathcal{S}$.

The GRT of $f$ is given by:

[TABLE]

where the weight $b$ is smooth (i.e., $C^{\infty}$) and non-vanishing, and $dx$ is the area element on $\mathcal{S}_{y}$. The discrete data are given by

[TABLE]

for some $r\in\mathbb{R}^{3}$.

Even though (2.8) assumes that the stepsize along each data axis equals $\epsilon$, this is a non-restrictive assumption. Indeed, consider a smooth diffeomorphism $\psi$: $\mathcal{V}\to\tilde{\mathcal{V}}$ for some open $\tilde{\mathcal{V}}\subset\mathbb{R}^{2}$, so that $\psi$ maps an irregular grid covering $\mathcal{V}$ into a regular, square grid covering $\tilde{\mathcal{V}}$. Introducing a new defining function $\tilde{\Phi}(x,\tilde{y}):=\Phi(x,\psi^{-1}(\tilde{y}))$, we can transform any smoothly sampled data set into the one with a square grid. Clearly, if $\Phi$ satisfies DF1–DF4, then $\tilde{\Phi}$ satisfies DF1–DF4 as well.

Conditions DF1–DF4 imply that (cf. [2] and [20], Sections VIII.5 and VIII.6):

1. (1)

The GRT $\mathcal{R}$ is a Fourier Integral Operator (FIO) with phase function $\lambda\Phi(x,y)$; 2. (2)

The corresponding canonical relation is

[TABLE]

which is a local canonical graph; 3. (3)

Any suitably modified adjoint of $\mathcal{R}$, denoted $\mathcal{R}^{*}$, is also an FIO, whose canonical relation $C^{*}$ is obtained from (2.9) by switching the $(x,\xi)\in T^{*}\mathcal{U}$ and $(y,\eta)\in T^{*}\mathcal{V}$ variables; and 4. (4)

The composition $\mathcal{R}^{*}(\dots)\mathcal{R}$, where the dots denote a cut-off combined with a suitable differential operator, is a pseudo-differential operator ($\Psi$DO), i.e. $C^{*}\circ C$ is a subset of the diagonal in $T^{*}\mathcal{U}$.

Given a point $x\in\mathcal{S}$, find $y=y(x)$ (which is smooth locally) such that $\mathcal{S}_{y}$ is tangent to $\mathcal{S}$ at $x$. Denote

[TABLE]

where $\text{II}_{\mathcal{S}}(x)$ is the matrix of the second fundamental form of $\mathcal{S}$ at $x\in\mathcal{S}$ written in an orthonormal basis of $T_{x}\mathcal{S}$.

For any $x\in\mathcal{U}$, introduce the set

[TABLE]

Definition 1**.**

A pair $(x_{0},y_{0})\in\mathcal{C}$ is globally generic if whenever $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at some $z\not=x_{0}$ the following conditions hold:

• GG1.

$\mathcal{S}$* is smooth at $z$, and $N(z)$ is either positive or negative definite; and*

• GG2.

Let $\dot{\Gamma}_{x_{0}}$ be a non-vanishing at any point tangent vector field along $\Gamma_{x_{0}}$. There exists an open set $\mathcal{V}_{1}$, $y_{0}\in\mathcal{V}_{1}\subset\mathcal{V}$, such that for each $m\in\mathbb{Z}^{3}$, $|m|>0$, and all $\delta>0$ sufficiently small,

1. (1)

The set $\{y\in\Gamma_{x_{0}}\cap\mathcal{V}_{1}:\,|m\cdot\dot{\Gamma}_{x_{0}}(y)|\leq\delta\}$ is contained in a finite number of segments of $\Gamma_{x_{0}}$ (this number may depend on $m$ and $\delta$), and 2. (2)

The sum of the lengths of these segments goes to zero as $\delta\to 0$.

As is shown in Section 7, Conditions GG1 and DF3 imply that $\Gamma_{x_{0}}$ is a smooth curve, so Condition GG2 makes sense.

An example when Condition GG2 is violated is when $\Gamma_{x_{0}}$ contains a straight line segment and $m\cdot\dot{\Gamma}_{x_{0}}(y)\equiv 0$ on this segment for some $m\in\mathbb{Z}^{3}$, $|m|>0$.

Definition 2**.**

A pair $(x_{0},y_{0})\in\mathcal{C}$, $x_{0}\in\mathcal{S}$, is locally generic if whenever $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}$ the following conditions hold:

• LG1.

$S$* is smooth at $x_{0}$, and $N(z)$ is either positive or negative definite; and*

• LG2.

There is no $\lambda\not=0$ such that $\lambda\Phi^{\prime}_{y}(x_{0},y_{0})\in\mathbb{Z}^{3}$.

Definition 3**.**

A pair $(x_{0},y_{0})\in\mathcal{C}$ is generic if it is both locally and globally generic.

Let $\varphi$ be an interpolating kernel (IK), i.e. $\varphi(0)=1$ and $\varphi(j)=0$ for all $j\in\mathbb{Z}$, $j\not=0$. Suppose also that $\varphi$ satisfies the following assumptions:

• IK1.

$\varphi$ is exact up to the order $2$, i.e.

[TABLE]

• IK2.

$\varphi$ is compactly supported;

• IK3.

All partial derivatives of $\varphi$ up to the order $2$ are continuous;

• IK4.

All partial derivatives of $\varphi$ of order $3$ are piecewise continuous and bounded; and

• IK5.

$\varphi$ is normalized, i.e. $\int\varphi(y)dy=1$.

The interpolated version of $g$ can be written in the form

[TABLE]

First, we derive a microlocal inversion formula for the GRT $\mathcal{R}$, which reconstructs exactly the leading singularities of $f$. Pick any $(x_{0},y_{0})\in\mathcal{C}$. Let $\alpha_{0}$ be a unit vector normal to $\mathcal{S}_{y_{0}}$ at $x_{0}$. For $(x,\alpha)\in\mathcal{U}\times S^{2}$ close to $(x_{0},\alpha_{0})$ and for $t$, $|t|\ll 1$, find the local solution $y=Y(\alpha,t;x)$ such that $x+t\alpha\in\mathcal{S}_{y}$ and $\alpha$ is normal to $\mathcal{S}_{y}$ at $x+t\alpha$. By construction, $y_{0}=Y(\alpha_{0},t=0;x_{0})$. Here we use the assumption that the data are complete, i.e. such a solution exists. It is shown below (see (4.2)) that the map $(\alpha,t)\to y=Y(\alpha,t;x)$ is a local diffeomorphism that depends smoothly on $x$.

Let $\mathcal{V}_{1}$ be a small neighborhood of $y_{0}$. Pick any $\chi\in C_{0}^{\infty}(\mathcal{V}_{1})$ such that $\chi\equiv 1$ near $y_{0}$. The inversion formula with continuous data is given by

[TABLE]

This inversion formula emulates the CRT inversion formula by backprojecting a second order derivative of the GRT. The affine variable $t$ is computed relative to $x$ as opposed to the origin, as is the case with the CRT. Hence the GRT analogue of the usual term $\alpha\cdot x$ is missing from (2.14), because it is absorbed by the function $Y$. Due to the symmetry $g(Y(\alpha,t;x))=g(Y(-\alpha,-t;x))$, in (2.14) we integrate over half of the unit sphere $S_{+}^{2}$.

Using the argument following (2.8), it is easy to show that the map $f\to f_{\chi}$ is a $\Psi$DO of degree zero with principal symbol 1 microlocally near $(x_{0},\alpha_{0})$ (see e.g. [2, 9]). Thus, the singularities of $f$ and $f_{\chi}$ are the same to leading order (e.g., in the scale of Sobolev spaces) microlocally near $(x_{0},\alpha_{0})$. An inversion formula that recovers all the singularities of $f$ can be obtained by combining (2.14) with a microlocal partition of unity. In the case of discrete data, we use the same inversion formula (2.14), but replace $g$ with $g_{\epsilon}$. The corresponding reconstruction is denoted $f_{\chi\epsilon}$.

3. Statement of main result. Beginning of proof

3.1. Statement of main result.

Pick $(x_{0},y_{0})\in\mathcal{C}$ such that $\mathcal{S}$ is smooth at $x_{0}\in\mathcal{S}$, $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}$, and $N(x_{0})$ is either positive definite or negative definite. Fix some orthonormal basis in the common plane tangent to $\mathcal{S}$ and $\mathcal{S}_{y}$ at $x_{0}$. Let $\alpha_{0}$ be the unit vector normal to $\mathcal{S}$ at $x_{0}$. For convenience of calculations, we sometime use the Cartesian coordinates $(x_{1},x^{\perp})$ determined by

[TABLE]

where $\alpha_{0}^{\perp}$ is the plane through $x_{0}$ and normal to $\alpha_{0}$. This plane is tangent to both $\mathcal{S}$ and $\mathcal{S}_{y_{0}}$ at $x_{0}$. The direction of $\alpha_{0}$ is chosen so that $N(x_{0})$ is negative definite. The side of $\mathcal{S}$ where $\alpha_{0}$ points is called interior. The other side of $\mathcal{S}$ is called exterior. If necessary, multiply $\Phi$ by $(-1)$ so that $\Phi^{\prime}_{x}(x_{0},y_{0})/|\Phi^{\prime}_{x}(x_{0},y_{0})|=-\alpha_{0}$.

Consider the point

[TABLE]

where $\tilde{\mathcal{U}}$ is a bounded set. Denote:

[TABLE]

Here and in what follows the convention is that if the arguments of $\Phi$ and its derivatives are omitted, then they are evaluated at $x_{0},y_{0}$.

We also introduce local $y$-coordinates with the origin at $y_{0}$:

[TABLE]

Thus, equation $y_{1}=0$ determines the plane tangent to the submanifold $\mathcal{T}_{x_{0}}$ at $y_{0}$. We frequently denote this plane $\beta_{0}^{\perp}$.

Finally, we use the extension of the CRT to all of $\mathbb{R}^{3}\setminus 0$ according to:

[TABLE]

for a sufficiently smooth $f$.

The main result of the paper is the following

Theorem 1**.**

Pick a generic pair $(x_{0},y_{0})\in\mathcal{C}$ such that $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}\in\mathcal{S}$. Then

[TABLE]

where $h=\tilde{x}\cdot\alpha_{0}$, and $\hat{\varphi}$ is the CRT of $\varphi$.

By IK5, $\int_{\mathbb{R}}\hat{\varphi}(\beta_{0},s)ds=1$. The inversion formula (2.14) reconstructs jumps of $f$ accurately, so $f_{\chi}(x_{0^{+}})-f_{\chi}(x_{0^{-}})=f_{0}$ and (3.8) can be written as follows

[TABLE]

By linearity, the proof of Theorem 1 can be split into two parts: local and global. The local part is formulated as follows.

Theorem 2**.**

Pick a locally generic pair $(x_{0},y_{0})\in\mathcal{C}$ such that $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}\in\mathcal{S}$. Suppose $\text{supp}(f)$ is contained in a sufficiently small neighborhood of $x_{0}$. Then

[TABLE]

where $h=\tilde{x}\cdot\alpha_{0}$, and $\hat{\varphi}$ is the CRT of $\varphi$.

Since $\mathcal{R}$ is an FIO with canonical relation (2.9), $g=\mathcal{R}f$ is singular only when $\mathcal{S}_{y}$ is tangent to $\mathcal{S}$. Therefore, we are interested in the behavior of $g$ in a neighborhood of $y_{0}$.

Remark**.**

Strictly speaking, one has to distinguish between the original coordinates that describe points $x\in\mathcal{U}$, $y\in\mathcal{V}$ and those in (3.1), (3.4), respectively. For example, one should write $x=\hat{x}_{1}\alpha_{0}+\hat{x}^{\perp}$ instead of $x=x_{1}\alpha_{0}+x^{\perp}$. Such notation would emphasize that $\hat{x}_{1}$ is not the first component of $x$ in the original coordinates, i.e. $\hat{x}_{1}\not=x_{1}$. Similarly, a derivative like $\partial\Phi(x,y)/\partial x_{1}$, if written in full, becomes $\partial\Phi(x(\hat{x}_{1},\hat{x}^{\perp}),y(\hat{y}_{1},\hat{y}^{\perp}))/\partial\hat{x}_{1}$. To avoid burdensome notations, whenever the coordinates $(\hat{x}_{1},\hat{x}^{\perp})$ and $(\hat{y}_{1},\hat{y}^{\perp})$ are used, we will stick with the simplified notation and assume that the above convention holds.

3.2. Behavior of $g$ near its singular support

Because $\mathcal{S}$ is smooth in a neighborhood of $x_{0}\in\mathcal{S}$, there is a smooth local diffeomorphism $x\to(z,p)$ so that

[TABLE]

The normal $n(z)$ is chosen so that $N(z)$ is negative definite. Thus $n(x_{0})=\alpha_{0}$. Clearly, we can extend the function $n(z)$, $z\in\mathcal{S}$, to $n(x)$ defined in a neighborhood of $\mathcal{S}$ by the formula $n(z+pn(z)):=n(z)$. With a slight abuse of notation, the extended function will also be denoted $n(\cdot)$.

Using (3.9), define $\Psi(z+pn(z)):=p$. Then $\Psi(x)=0$ is the equation of $\mathcal{S}$ near $x_{0}$, and $\Psi$ is smooth. By construction, $\Psi(x)>0$ on the interior side of $\mathcal{S}$.

Consider the system of equations

[TABLE]

If we set $p=0$ and solve (3.10) for $y$, we find submanifolds $\mathcal{S}_{y}$ tangent to $\mathcal{S}$. We also need to solve these equations for $z,p$, and $\lambda$ in terms of $y$.

Lemma 1**.**

Pick $(x_{0},y_{0})\in\mathcal{C}$ such that (1) $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}\in\mathcal{S}$, (2) $\mathcal{S}$ is smooth at $x_{0}$, and (3) $N(x_{0})$ is negative definite. There exists an open set $\mathcal{V}_{1}$, $y_{0}\in\mathcal{V}_{1}\subset\mathcal{V}$, such that

1. (1)

The set of $y\in\mathcal{V}_{1}$ such that $\mathcal{S}_{y}$ is tangent to $\mathcal{S}$ in a neighborhood of $x_{0}$ is a smooth submanifold of $\mathcal{V}$ through $y_{0}$. The vector $\Phi^{\prime}_{y}(x_{0},y_{0})$ is normal to this submanifold at $y_{0}$. 2. (2)

The solutions $z=Z(y)$, $p=P(y)$, and $\lambda=\Lambda(y)$ to (3.10) depend smoothly on $y\in\mathcal{V}_{1}$, and

[TABLE] 3. (3)

Equations (3.10) determine a smooth function $y=\overline{Y}(z,p)$, $(z,p)\in\mathcal{S}\times\mathbb{R}$, in a neighborhood of $(x_{0},0)$.

Proof.

Differentiate (3.10) with respect to $z$, $p$, $\lambda$, and $y$, and set $z=x_{0}$, $p=0$, $y=y_{0}$ to obtain the $5\times 8$ matrix

[TABLE]

Here $n_{x}^{\prime}$ is the derivative of the function $n(x)$ extended to a neighborhood of $\mathcal{S}$ as described following (3.9). Since $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}\in\mathcal{S}$, we have $\Phi^{\prime}_{x}\parallel\alpha_{0}$, so $\Phi^{\prime}_{x}\cdot\alpha_{0}=\lambda=-|\Phi^{\prime}_{x}|\not=0$. This also gives the value of $\lambda$ to be used in (3.12).

To prove the first part of the first assertion we need to show that $z,\lambda$, and $y_{1}$ are smooth functions of $y^{\perp}$. Remove the columns corresponding to the derivatives with respect to $p$ (because $p=0$ is fixed) and $y^{\perp}$ to obtain a $5\times 5$ submatrix. Calculation in coordinates shows that $\Psi_{x}^{\prime}(z+pn(z))\equiv n(z)$ (cf. (3.9)). By applying elementary row and column operations, it is clear that this submatrix is full-rank if and only if the following matrix has rank two:

[TABLE]

Here $\Phi^{\prime\prime}_{x^{\perp}x^{\perp}}:\alpha_{0}^{\perp}\to\alpha_{0}^{\perp}$ is the appropriate submatrix of $\Phi^{\prime\prime}_{xx}$ in the coordinates (3.1), and $n^{\perp}$ is the projection of $n$ onto $\alpha_{0}^{\perp}$. It is easy to see that

[TABLE]

The desired assertion follows from Condition LG1 (see Definition 2).

Next, set $p=0$ in (3.10) and assume that $z,y_{1}$ are functions of $y^{\perp}$. Differentiating the last two equations in (3.10) with respect to $y^{\perp}$ and using that $\Phi^{\prime}_{y}\cdot\beta_{0}=|\Phi_{y}^{\prime}|\not=0$ gives $\Phi^{\prime}_{x}z^{\prime}_{y^{\perp}}=0$ and $\partial y_{1}/\partial y^{\perp}=0$. This proves the second part of the first assertion.

The first part of the second assertion follows by retaining the columns corresponding to the derivatives with respect to $z,p$, and $\lambda$. As before, the resulting $5\times 5$ submatrix is full-rank because the matrix in (3.13) has rank two. The second part of the second assertion follows by considering $z,p$, and $\lambda$ as functions of $y$, differentiating the last two equations in (3.10) with respect to $y$, and using that $\Phi^{\prime}_{x},\Phi^{\prime}_{y}\not=0$.

The third assertion follows immediately from data completeness and the Bolker condition (Conditions DF2 and DF4, respectively). ∎

We need the following lemma, which generalizes one of the results of Ramm and Zaslavsky [16, 17] from the CRT to the GRT.

Lemma 2**.**

Pick $(x_{0},y_{0})\in\mathcal{C}$ such that (1) $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $x_{0}\in\mathcal{S}$, (2) $\mathcal{S}$ is smooth at $x_{0}$, and (3) $N(x_{0})$ is negative definite. Suppose $\text{supp}(f)$ is contained in a small neighborhood of $x_{0}$. For any $z\in\mathcal{S}$ and $p$ in small neighborhoods of $x_{0}$ and [math], respectively, one has:

[TABLE]

for some smooth $G(z,p)$, $G_{1}(z,p)$.

Proof.

Recall that $\overline{Y}(z,p)$ is the smooth function of $(z,p)\in\mathcal{S}\times\mathbb{R}$ defined by the conditions that $z+pn(z)\in\mathcal{S}_{y}$ and $n(z)$ be normal to $\mathcal{S}_{y}$ at the point $z+pn(z)$, see Assertion (3) of Lemma 1. Recall also that in coordinates (3.1), $N(x_{0})$ is negative definite. By linearity, we may assume that $f\equiv 0$ on the exterior side of $\mathcal{S}$. In particular, $f(x_{0}-\epsilon\alpha_{0})\equiv 0$, $\epsilon>0$, in (3.3). In this case we have to prove (3.15) with $G_{1}\equiv 0$. By construction,

[TABLE]

where $\theta$ is the unit step function (Heaviside function). Consider the system

[TABLE]

which we solve for $x$ and $\mu$ in terms of $y$. Equations (3.17) determine the stationary point $x_{*}(y)$ of $\Psi(x)$ on the surface $\mathcal{S}_{y}$ (the parameter $\mu$, which corresponds to $1/\lambda$ in (3.10), is the Lagrange multiplier). From (3.10) and the property of $\Psi^{\prime}_{x}$, the solution $x_{*}(y)$, $\mu(y)$ to (3.17) can be obtained from the solution $z=Z(y),p=P(y),\lambda=\Lambda(y)$ to (3.10): $x_{*}(y)=Z(y)+P(y)n(Z(y))$, $\mu(y)=1/\Lambda(y)$. By Assertion (2) of Lemma 1, $x_{*}(y)$ is smooth near $y_{0}$.

As is easily checked, when the matrix $\Psi^{\prime\prime}_{xx}(x_{0})-\mu(y_{0})\Phi^{\prime\prime}_{xx}(x_{0},y_{0})$ is viewed as a quadratic form on $\alpha_{0}^{\perp}$, it coincides with $N(x_{0})$. The latter is negative definite, hence $x_{*}(y)$ is the local maximum of $\Psi(x)$ on $\mathcal{S}_{y}$. By the Morse lemma, find local coordinates $\omega$ on $\mathcal{S}_{y}$, which depend smoothly on $y$, such that $\omega(x_{*}(y))=0$ and $\Psi(x)=\Psi(x_{*}(y))-|\omega|^{2}$, $x=x(\omega;y)\in\mathcal{S}_{y}$. Since $f,b$, and $x(\omega;y)$ are all smooth, we get from (3.16)

[TABLE]

for some smooth $F$. Expand $F$ in the Taylor series around $\omega=0$ and integrate in spherical coordinates $\omega=r\Theta$. Integration with respect to $\Theta$ removes all the odd powers of $r$, i.e. only the even powers of $r$ remain. By construction,

[TABLE]

and the first statement in (3.15) follows.

To prove the second statement, we use the local coordinates (3.1). In these coordinates, the local equation of $\mathcal{S}_{y}$ becomes

[TABLE]

for some smooth $\phi$, $a$, and $A$. Suppose $y$ is such that $x_{*}(y)=x_{0}+p\alpha_{0}$. Substitute such a pair $(x_{*}(y),y)$ into (3.17) and use (3.20) to conclude that $\phi(y)\equiv p$ and $a(y)\equiv 0$. Let $x_{1}=Q(x^{\perp})$ be the equation of $\mathcal{S}$ in the coordinates (3.1). By construction, $Q^{\prime}(0)=0$. Clearly, $A(y_{0})$ and $Q^{\prime\prime}(0)$ are the matrices of the second fundamental form of $\mathcal{S}_{y_{0}}$ and $\mathcal{S}$, respectively, at $x_{0}$ in the coordinates (3.1). From (3.16),

[TABLE]

Integrating in (3.21) by diagonalizing $N=\psi^{\prime\prime}_{x^{\perp}x^{\perp}}(0,y_{0})$ and changing variables, we find

[TABLE]

which finishes the proof. ∎

4. Local behavior of interpolated data.

Similarly to (3.10), consider the equations

[TABLE]

which we solve to find $y=Y(\alpha,t;x)$ assuming $(\alpha,t,x)$ is in a neighborhood of $(\alpha_{0},0,x_{0})$. Differentiating (4.1) with respect to $\lambda$ and $y$ and setting $x=x_{0}$, $t=0$, $\alpha=\alpha_{0}$, $\lambda=-|\Phi^{\prime}_{x}|$, we obtain similarly to (3.12) a matrix, which is non-degenerate. As opposed to (3.12), the key reason why it is non-degenerate is the Bolker condition. Hence $Y(\alpha,t;x)$ is a smooth function of $\alpha,t$ and $x$. Next we substitute $y=Y(\alpha,t;x)$ into (4.1) and obtain a few useful properties of $Y$. The first one is that $\partial Y_{1}/\partial\alpha^{\perp}=0$. Indeed, differentiate the second equation in (4.1) with respect to $\alpha^{\perp}$ and set $x=x_{0}$, $t=0$, $\alpha=\alpha_{0}$ to obtain $\Phi^{\prime}_{y}\partial Y/\partial\alpha^{\perp}=0$. The desired assertion follows from (3.4). In a similar fashion, we have

[TABLE]

The first result is obtained by differentiating the first equation in (4.1) with respect to $\alpha^{\perp}$ and using the Bolker condition and that $\partial Y_{1}/\partial\alpha^{\perp}=0$. The second result is obtained by differentiating the second equation in (4.1) with respect to $t$ and using the properties of the selected coordinates (3.1), (3.4). The last result is an obvious consequence of the first two and that $\partial Y_{1}/\partial\alpha^{\perp}=0$.

In view of (4.2), given any small $\omega>0$, we can find a sufficiently small open set $\mathcal{V}_{1}$, $y_{0}\in\mathcal{V}_{1}\subset\mathcal{V}$, such that $Y(\alpha,t;x_{\epsilon})\in\mathcal{V}_{1}$ implies $|\alpha^{\perp}|<\omega$ for all $t=O(\epsilon)$ and $x_{\epsilon}$ provided that $\epsilon$ is sufficiently small. This relationship between $\mathcal{V}_{1}$ and $\omega$ is assumed in what follows.

Since $\text{supp}(\chi)\subset\mathcal{V}_{1}$, the integral with respect to $\alpha$ in (2.14) can be split into two sets:

[TABLE]

for some small (but fixed) $\omega>0$. Here $A>0$ is a large parameter. Let $f_{\chi\epsilon}^{(j)}$ denote the result of integrating in (2.14) (with $g$ replaced by $g_{\epsilon}$) over $\Omega_{j}$, $j=1,2$. The behavior of $f_{\chi\epsilon}^{(1)}$ is investigated first. This is done in Section 5. In the remainder of this section we lay the groundwork for that investigation by deriving the behavior of $Y(\alpha,t;x_{\epsilon})$ and $g_{\epsilon}(y)$ in a neighborhood of $(\alpha_{0},0,x_{0})$ and $y_{0}$, respectively.

4.1. Local behavior of $Y(\alpha,t;x_{\epsilon})$.

The first step is to obtain the leading term behavior of the function $Y(\alpha,t;x_{\epsilon})$ for $t=O(\epsilon)$ and $\alpha\in\Omega_{1}$, i.e. for $|\alpha^{\perp}|=O(\epsilon^{1/2})$. In this section we continue using the coordinates (3.1).

Expanding $y=Y(\alpha,t;x_{\epsilon})$ in the Taylor series around $x=x_{0}$, $t=0$, $\alpha=\alpha_{0}$ and using that $|x_{\epsilon}-x_{0}|=O(\epsilon)$, $t=O(\epsilon)$, $|\alpha^{\perp}|=O(\epsilon^{1/2})$, and $\partial Y_{1}/\partial\alpha^{\perp}=0$ gives

[TABLE]

For $x$ in an $O(\epsilon)$ neighborhood of the origin (i.e., $x_{0}$) and for $y$ in an $O(\epsilon^{1/2})$ neighborhood of the origin (i.e., $y_{0}$) we have

[TABLE]

To find $y=Y(\alpha,t;x_{\epsilon})$, substitute $x=x_{\epsilon}$ into (4.5) and solve

[TABLE]

Recall that $h=\tilde{x}\cdot\alpha_{0}$. Switching to the coordinates (3.1), (3.4), using (4.4), and keeping only the terms of order $O(\epsilon)$ in the first equation in (4.6) gives

[TABLE]

Projecting the second equation in (4.6) onto $\alpha_{0}$ and onto $\alpha_{0}^{\perp}$ implies

[TABLE]

leading to

[TABLE]

By the Bolker condition (2.4), $\Phi^{\prime\prime}_{x^{\perp}y^{\perp}}$ is nondegenerate. Substitution into (4.7) now yields:

[TABLE]

Recall that $\nu$ is defined in (3.3).

4.2. The leading local behavior of $g(Y(\alpha,t;x_{\epsilon}))$

We plan to substitute $y=Y(\alpha,t=0;x_{\epsilon})$ into (2.13). Hence the second step is to find the leading behavior of $g(y)$ when $|y-Y(\alpha,0;x_{\epsilon})|=O(\epsilon)$ and $|\alpha^{\perp}|=O(\epsilon^{1/2})$. This is done by finding the asymptotics of $z=Z(y)$ and $p=P(y)$, which are determined by solving (3.10).

Recall that the local equation of $\mathcal{S}$ in the coordinates (3.1) and the interior unit normal are given by

[TABLE]

By (4.9), (4.10), $|y_{1}|=O(\epsilon)$, $|y^{\perp}|=O(\epsilon^{1/2})$. From (3.11), $\partial p/\partial y^{\perp}=0$ at $y=y_{0}$, so this implies $|p|=O(\epsilon)$, $x_{1}=O(\epsilon)$, and $|z^{\perp}|,|x^{\perp}|=O(\epsilon^{1/2})$, where $x=z+pn(z)$. The equation (4.11) leads to:

[TABLE]

From (4.12), $x^{\perp}=z^{\perp}+O(\epsilon^{3/2})$.

Given that now $|x^{\perp}|=O(\epsilon^{1/2})$, the expansion in the first line in (4.5) should include additional terms. The second equation in (3.10) becomes

[TABLE]

Solving for $x_{1}$ we find:

[TABLE]

Using (4.9) and that $|y-Y(\alpha,0;x_{\epsilon})|=O(\epsilon)$, we have

[TABLE]

and the unit normal vector to $\mathcal{S}_{Y}$ is thus:

[TABLE]

The big-$O$ terms in (4.16) follow by noticing that differentiation with respect to $x_{1}$ in (4.13), (4.14) converts $O(\epsilon^{3/2})$ into $O(\epsilon^{1/2})$, and differentiation with respect to $x^{\perp}$ converts $O(\epsilon^{3/2})$ into $O(\epsilon)$. This follows by looking at the terms absorbed by $O(\epsilon^{3/2})$.

From (4.11) and (4.16), the two normal vectors are parallel (the first equation in (3.10)) if

[TABLE]

which implies

[TABLE]

Matching the first components in (4.12) and using (4.14), (4.15), (4.18) gives after simple transformations:

[TABLE]

where $M$ is defined in (4.10). Summarizing (4.18) and (4.19) we have:

[TABLE]

Substituting into (3.15) we find

[TABLE]

where $y^{*}=\overline{Y}(z,p=0)$ whenever $y=\overline{Y}(z,p)$.

From (4.11), (4.18) and (4.9), (4.10) it follows that $|Z(y)-x_{0}|=O(\epsilon^{1/2})$ and $|y-y_{0}|=O(\epsilon^{1/2})$ whenever $|y-Y(\alpha,0;x_{\epsilon})|=O(\epsilon)$, hence

[TABLE]

4.3. Local behavior of the interpolated data.

In this subsection we find the behavior of the interpolated data near $y_{0}=0$. Combining (4.20)–(4.22) we get

[TABLE]

where $M_{1}$ is defined in (4.19). From (4.10) and (4.20), $|P(\epsilon j)|=O(\epsilon)$. Using (4.9), (4.10) and (4.23) gives:

[TABLE]

Since $\varphi$ is compactly supported, the number of terms in the sum in (4.24) is bounded. Using additionally that $|\tilde{\alpha}^{\perp}|\leq A<\infty$, the sum itself is bounded as well. Finally, combining with the fact that $|(a+O(\epsilon^{1/2}))_{+}-a_{+}|=O(\epsilon^{1/2})$ uniformly in $a\in\mathbb{R}$ and using (4.9) and (4.10) gives

[TABLE]

In view of (4.25), denote

[TABLE]

Here we will need a higher order approximation of $Y(\alpha,t;x_{\epsilon})$ than the one in (4.9), (4.10):

[TABLE]

where $A_{1}$ is a bilinear map $\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$, and $A_{2},A_{3}\in\mathbb{R}^{2}$. Moreover, the two $O(\epsilon^{1/2})$ terms in (4.27) depend smoothly on $\tilde{t}$ and $\tilde{\alpha}^{\perp}$. In particular, differentiation with respect to $\tilde{t}$ does not change the order of these terms as $\epsilon\to 0$. Hence we rewrite (4.25) as follows

[TABLE]

We also have

[TABLE]

5. Estimating the term $f_{\chi\epsilon}^{(1)}$.

To study $f_{\chi\epsilon}^{(1)}$ we need the following lemma, which follows immediately from (4.26) and the properties IK1–IK3 of $\varphi$.

Lemma 3**.**

Partial derivatives of $\psi(q,u)$ with respect to $u$ up to the order two are continuous. Also, one has

[TABLE]

and, for some $c>0$,

[TABLE]

Denote

[TABLE]

The following lemma is a direct consequence of Lemma 3 (see also properties IK3, IK4 of $\varphi$).

Lemma 4**.**

The function $U(q,u)$ has piecewise continuous bounded first order partial derivatives with respect to $q$ and $u$. Also,

[TABLE]

and, for some $c>0$,

[TABLE]

Note that the derivative in the inversion formula (2.14) is with respect to $t$. Using that $t=\epsilon\tilde{t}$ (cf. (4.24)) and taking into account (4.29), (5.3), in the formula below we will acquire the factor $(\epsilon/\nu)^{2}$. Thus, in terms of $U$, the expression for $f_{\chi\epsilon}^{(1)}$ becomes after changing variables $\alpha^{\perp}\to\tilde{\alpha}^{\perp}$ (this brings the factor $\epsilon$), setting $\tilde{t}=0$, and using (4.28), (4.29), (5.3):

[TABLE]

Two simplifications have been made in deriving (5.6). First, using that first and second order derivatives of $\psi$ are bounded, it follows from (4.29) that

[TABLE]

Second, since the derivatives of $U$ are bounded, the integral with respect to $\tilde{\alpha}^{\perp}$ is over a bounded set, and the coefficient in front of the integral is bounded, using (5.7) and then replacing $u(\tilde{\alpha}^{\perp},0)$ with $u_{0}(\tilde{\alpha}^{\perp})$ in the arguments of $U$ leads to the term $O(\epsilon^{1/2})$ outside the integral.

Approximate the domain $|\tilde{\alpha}^{\perp}|<A$ by a union of non-overlapping small squares of size $\delta$. Let these squares be denoted $B_{k}$, $k=1,2,\dots,O(\delta^{-2})$. By (5.6),

[TABLE]

where $\tilde{\alpha}^{\perp}_{k}$ is the center of $B_{k}$. By Lemma 4,

[TABLE]

Therefore

[TABLE]

where we have used that $\Delta u_{k}(\tilde{\alpha}^{\perp})\cdot\beta_{0}=0$. In (5.10) and below the fractional part of a vector is computed component wise: $\{u\}=(\{u_{1}\},\{u_{2}\},\{u_{3}\})$, where $\{u_{i}\}=u_{i}-\lfloor u_{i}\rfloor$ and $\lfloor u_{i}\rfloor$ is the largest integer not exceeding $u_{i}$.

Pick any $m\in\mathbb{Z}^{3}$, $m\not=0$. Let $m^{\perp}$ be the projection of $m$ onto the plane $\beta_{0}^{\perp}$. Condition LG2 in Definition 2 implies that $m^{\perp}\not=0$. By the local Bolker condition DF4, $\Phi^{\prime\prime}_{x^{\perp}y^{\perp}}$ is non-degenerate, so the vector $(\Phi^{\prime\prime}_{x^{\perp}y^{\perp}})^{-T}m^{\perp}\in\mathbb{R}^{2}$ is not zero. Using (5.8), the standard Weyl-type argument (cf. [13]) implies that

[TABLE]

Indeed, consider the function $U_{1}(q,\omega):=U(q-\omega\cdot\beta_{0},\omega)$. Clearly, $U_{1}(q,\omega)$ is periodic: $U_{1}(q,\omega)=U_{1}(q,\omega+m)$, $m\in\mathbb{Z}^{3}$. Expand $U_{1}(q,\omega)$ in a Fourier series:

[TABLE]

By the argument preceding (5.11),

[TABLE]

Therefore,

[TABLE]

In other words, only the term corresponding to $m=0$ survives, and the desired assertion follows from the standard approximation argument.

By (4.23), (4.28) and (5.6),

[TABLE]

Add the integrals over all the squares $B_{k}$ and use (5.6), (5.11), and (5.15):

[TABLE]

Since $\delta>0$ can be as small as we like, (5.16) implies

[TABLE]

The second argument of $U$ is bounded and $N$ is negative definite, so by (5.5) the integral with respect to $\tilde{\alpha}^{\perp}$ over the set $|\tilde{\alpha}^{\perp}|<A$, when $A>0$ is large enough (but fixed), can be replaced by the integral over all $\mathbb{R}^{2}$. Changing variables and integrating in spherical coordinates gives

[TABLE]

The integral with respect to $\omega$ can be evaluated explicitly. Recall that in our coordinates, $\beta_{0}=(1,0,0)$ (cf. (3.4)). From (4.26) and (5.3),

[TABLE]

Substitute (5.19) into (5.18)

[TABLE]

Since $\varphi$ is normalized and compactly supported,

[TABLE]

for some $c>0$.

Using the definition (3.5) and some simple transformations, we can rewrite the integral in (5.20) in two different forms

[TABLE]

6. Analysis of the term $f_{\chi\epsilon}^{(2)}$.

Lemma 5**.**

One can find $\omega>0$, $\epsilon_{0}>0$ small enough and $A>0$ large enough so that $g(y)$ is smooth in a neighborhood of all $y$ such that $(y-Y(\alpha,0;x_{\epsilon}))/\epsilon\in\text{supp}(\varphi)$ for any $\alpha\in\Omega_{2},\tilde{x}\in\tilde{U}$ (cf. (3.2)), and $0<\epsilon<\epsilon_{0}$.

Proof.

Fix some $c>0$ sufficiently large. Using that $N(x_{0})$ is negative definite, equation (4.20) implies that we can find $A>0$ large enough and $\omega>0$ small enough so that $P(Y(\alpha,0;x_{\epsilon}))>c\epsilon$ for all $\alpha\in\Omega_{2}$ and all $\tilde{x}\in\tilde{U}$ provided that $\epsilon$ is small enough. Since $\varphi$ is compactly supported and $P(y)$ is smooth, $P(y)>0$ for all $y$ such that $(y-Y(\alpha,0;x_{\epsilon}))/\epsilon\in\text{supp}(\varphi)$ (this is where we use that $c>0$ is sufficiently large). Therefore, $g$ is a smooth function in a neighborhood of all such $y$ because in this case we can drop the subscript ${}^{\prime}+^{\prime}$ from $P_{+}(y)$ in (3.15). ∎

Lemma 5 implies that the inversion formula (2.14) and its discrete analogue do not see singularities in the data when $\alpha\in\Omega_{2}$ provided that $\omega,A$ are selected as in the proof of Lemma 5. Since $|\Omega_{1}|=O(\epsilon)$, clearly the limit $\lim_{\epsilon\to 0^{+}}f_{\chi\epsilon}^{(2)}(x_{\epsilon})$ exists and is independent of $\tilde{x}\in\tilde{\mathcal{U}}$.

Using again that $N(x_{0})$ is negative definite, there exist sufficiently small neighborhoods $\mathcal{U}_{1}$ of $x_{0}$ and $\mathcal{V}_{1}$ of $y_{0}$ such that $\mathcal{S}_{y}\cap\mathcal{U}_{1}$ is on the exterior side of $\mathcal{S}$ whenever $\mathcal{S}_{y}$ is tangent to $\mathcal{S}$ and $y\in\mathcal{V}_{1}$. This implies that if $x\in\mathcal{U}_{1}$ and $x$ is on the interior side of $\mathcal{S}$, then there is no $y\in\mathcal{V}_{1}$ such that $\mathcal{S}_{y}$ contains $x$ and is tangent to $\mathcal{S}$. In turn, this implies that the data $g(y)$, $y=Y(\alpha,t=0;x)\in\mathcal{V}_{1}$, which is used to compute $f_{\chi}(x)$ is also smooth. Then, clearly, $\lim_{\epsilon\to 0^{+}}(f_{\chi\epsilon}^{(2)}(x_{0}+\epsilon h\alpha_{0})-f_{\chi}(x_{0}+\epsilon h\alpha_{0}))=0$ for any $h>0$, and

[TABLE]

This concludes the proof of the first part of the theorem.

7. Contribution of remote singularities

Suppose $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at some $z_{0}\in\mathcal{S}_{y_{0}}$, $z_{0}\not=x_{0}$, $\mathcal{S}$ is smooth at $z_{0}$, and $N(z_{0})$ is either positive definite or negative definite. Set $\alpha_{0}=\Phi^{\prime}_{x}(z_{0},y_{0})/|\Phi^{\prime}_{x}(z_{0},y_{0})|$, so that $y_{0}=Y(\alpha_{0},0;x_{0})$. As before, $\mathcal{V}_{1}$ is a small neighborhood of $y_{0}$, and $\text{supp}(\chi)\subset\mathcal{V}_{1}$.

As follows from assertion (1) of Lemma 1 (with $x_{0}$ replaced by $z_{0}$ as the point of tangency), the set of $y\in\mathcal{V}_{1}$ such that $\mathcal{S}_{y}$ is tangent to $\mathcal{S}$ near $z_{0}$ is a smooth submanifold of $\mathcal{V}_{1}$ through $y_{0}$, and the vector $\Phi^{\prime}_{y}(z_{0},y_{0})$ is normal to it at $y_{0}$. The local equation of the manifold is $P(y)=0$, where the function $P$ is the same as in Lemma 1. By assertion (2) of the lemma, $P^{\prime}_{y}(y_{0})\not=0$.

Additionally, $\mathcal{T}_{x_{0}}$ is another smooth submanifold through $y_{0}$, and $\Phi^{\prime}_{y}(x_{0},y_{0})$ is normal to it at $y_{0}$. By the assumption DF3 of no conjugate points, $\Phi^{\prime}_{y}(x_{0},y_{0})$ and $\Phi^{\prime}_{y}(z_{0},y_{0})$ are not parallel, so the intersection of the two submanifolds is a smooth curve $\Gamma_{x_{0}}$ through $y_{0}$. This is the same curve $\Gamma_{x_{0}}$, which was introduced in (2.11). From this argument it is easy to see that $\Gamma_{x}$ depends smoothly on $x$ near $x=x_{0}$.

Theorem 3**.**

Pick a globally generic pair $(x_{0},y_{0})\in\mathcal{C}$ such that $x_{0}\not\in\mathcal{S}$ and $\mathcal{S}_{y_{0}}$ is tangent to $\mathcal{S}$ at $z_{0}\in\mathcal{S}$. Suppose $\text{supp}(f)$ is contained in a sufficiently small neighborhood of $z_{0}$. One has

[TABLE]

Proof.

Suppose first that the reconstruction point is $x_{0}$. Since the reconstruction point is fixed, the dependence of various quantities on $x_{0}$ is omitted from notations in most places when there is no risk of confusion. Then

[TABLE]

By Lemma 2,

[TABLE]

Since $\varphi$ is compactly supported, we can expand the factor $P(y)$ in (7.3) in the Taylor series centered at $Y(\alpha,0)$. Let $L(y)$ be its linear term:

[TABLE]

We begin by looking at the expression, which is obtained by ignoring the second and higher order terms in the expansion of $P$:

[TABLE]

Clearly,

[TABLE]

Consider the most singular part of $J_{\epsilon}^{(1)}$, which is obtained by using the first term on the right in (7.6) and replacing $G(\epsilon j)$ with $G(Y(\alpha,0))$:

[TABLE]

In view of (7.4) and (7.7), similarly to (4.26) and (5.3), introduce the function

[TABLE]

Clearly,

1. (1)

$\psi$ is compactly supported in $q$ (by (2.12)), 2. (2)

$\psi$ has bounded first order partial derivatives, and 3. (3)

$\psi(q,v;\alpha)=\psi(q,v-m;\alpha)$ for any $m\in\mathbb{Z}^{3}$.

Using (7.8) in (7.7) yields:

[TABLE]

Introduce local coordinates $s=(s_{1},s_{2})$ on $S_{+}^{2}$ so that $s_{1}\equiv P(Y(\alpha,0))$ in a neighborhood of $\alpha_{0}$. As is shown at the beginning of this section, $\Gamma$ is the transverse intersection of the submanifolds $\mathcal{T}_{x_{0}}$ and $\{y\in\text{supp}(\chi):\,P(y)=0\}$. By the first equation in (4.2), $\alpha\to Y(\alpha,0)\in\mathcal{T}_{x_{0}}$ is a regular parametrization near $\alpha_{0}$. In (4.2), $y_{0}=Y(\alpha_{0},0;x_{0})$, and $Y^{\perp}$ is determined by the projection onto the plane $(\Phi_{y}^{\prime}(x_{0},y_{0}))^{\perp}$ (as opposed to $(\Phi_{y}^{\prime}(z_{0},y_{0}))^{\perp}$). Hence $P^{\prime}_{y}(y_{0})\not=0$ (cf. (3.11)) implies $\partial P(Y(\alpha,0))/\partial\alpha^{\perp}\not=0$ near $\alpha_{0}$. Therefore the preimage of $\Gamma\cap\mathcal{V}_{1}$, given by $\{\alpha\in S_{+}^{2}:\,P(Y(\alpha,0))=0,\,Y(\alpha,0)\in\mathcal{V}_{1}\}$ is also a smooth curve, and local coordinates $(s_{1},s_{2})$ with the required property do exist. Then

[TABLE]

In the first line, the integral is over the bounded set $\{s\in\mathbb{R}^{2}:\,Y(\alpha(s),0)\in\text{supp}(\chi)\}$. In the second line, the integral can be confined to a bounded set $\{(\tilde{s}_{1},s_{2})\in\mathbb{R}^{2}:\,|\tilde{s}_{1}|<\tilde{A},Y(\alpha(0,s_{2}))\in\text{supp}(\chi)\}$ for some $\tilde{A}>0$ large enough.

Lemma 6**.**

Let $D$ be a rectangle $D:=[a_{1},b_{1}]\times[a_{2},b_{2}]$. Consider a function $\psi\in C(D\times\mathbb{R}^{3})$. Suppose $\psi$ is periodic: $\psi(s,y)=\psi(s,y+m)$ for any $m\in\mathbb{Z}^{3}$ and $(s,y)\in D\times\mathbb{R}^{3}$. Let $Y:[a_{2},b_{2}]\to\mathbb{R}^{3}$ be a $C^{1}$ function with the following property. For any $m\in\mathbb{Z}^{3}$, $|m|>0$,

1. (1)

The set $\{s_{2}\in[a_{2},b_{2}]:\,|m\cdot Y^{\prime}(s_{2})|\leq\delta\}$ is contained in a finite number of intervals for all $\delta>0$ sufficiently small (this number may depend on $m$ and $\delta$), and 2. (2)

The sum of the lengths of these intervals goes to zero as $\delta\to 0$.

Then one has

[TABLE]

Proof.

Pick any $\delta_{1}>0$. Let $\Psi_{m}(s)$ denote the coefficients of the Fourier expansion of $\psi(s,y)$ with respect to $y$. We can find $M>0$ large enough and a partition of $D$ into sufficiently small rectangles such that

[TABLE]

Here $\tilde{\Psi}_{m}$ is an approximation of $\Psi_{m}$, which is constant on each rectangle of the partition. Thus, the lemma will be proven if we show that

[TABLE]

Using assumptions (1) and (2) of the lemma, partition $[a,b]$ into a finite collection of non-overlapping intervals so that (i) their union is as close to $[a,b]$ as we like, and (ii) in each of these intervals $m\cdot Y^{\prime}(s_{2})$ is bounded away from zero. The result now follows immediately. ∎

Clearly, Condition GG2 in Definition 1 is independent of the choice of the vector field $\dot{\Gamma}$ as long as it does not vanish at any point of $\Gamma$. In the $s$-coordinates, $Y(\alpha(0,s_{2}),0)$ is a regular parametrization of $\Gamma\cap\mathcal{V}_{1}$ because $|\partial Y^{\perp}/\partial\alpha^{\perp}|_{\alpha=\alpha_{0}}\not=0$ and $|\partial\alpha^{\perp}/\partial s|\not=0$ (see also the argument following (7.9)). The latter determinant is computed at $s$ such that $\alpha(s)=\alpha_{0}$. Therefore $\partial Y/\partial s_{2}$ never vanishes on $\Gamma\cap\mathcal{V}_{1}$. Condition GG2 implies that $Y(\alpha(0,s_{2}))$ satisfies Conditions (1), (2) in Lemma 6. Set

[TABLE]

Using the properties (1)–(3) of $\psi$, we see that Lemma 6 applies to $\psi_{1}$. Also, $\psi_{1}$ is compactly supported. Compact support along $\tilde{s}_{1}$ is due to the property (1) of $\psi$, and along $s_{2}$ – due to the cut-off $\chi$. Substituting $\psi_{1}$ into (7.10), using Lemma 6, and then expressing $\psi_{1}$ in term of $\psi$ yields

[TABLE]

By (7.8), similarly to (5.19),

[TABLE]

As $|e|$ not necessarily equals one, (7.16) assumes the extended definition of the CRT, cf. (3.5). By (7.6) and (7.8),

[TABLE]

With $\varphi$ normalized, using (7.16) with $q=\tilde{s}_{1}$ and (7.17) in (7.15) gives

[TABLE]

Consequently, from (7.10) we get

[TABLE]

In the second line we used that $s=(s_{1},s_{2})$ and $s_{1}\equiv P(Y(\alpha(s),0))$.

Next, consider the second part of $J_{\epsilon}^{(1)}$, which is obtained by using the second term on the right in (7.6) and replacing $G(\epsilon j)$ with $G(Y(\alpha,0))$:

[TABLE]

The function $L_{+}(y)$ has bounded first derivatives, hence the limit of $J_{\epsilon}^{(1b)}$ can be easily found:

[TABLE]

Recall that $P>0$ on the interior side of $\mathcal{S}$.

The final piece of $J_{\epsilon}^{(1)}$ is

[TABLE]

In an $O(\epsilon)$ neighborhood of $\Gamma$, we have $|L(y)|=O(\epsilon)$. For any fixed $y\in\mathcal{V}_{1}$, $P(y)>0$, we compute by dropping the subscript ‘$+$’ from $L$:

[TABLE]

Substitution into (7.22) gives

[TABLE]

We can apply the limit as $\epsilon\to 0$ inside the integral in (7.22) to obtain (7.24) because the integrand is uniformly bounded. This follows because the function $(G(y)-G(Y(\alpha,0)))L(y)$ is smooth away from an $O(\epsilon)$ neighborhood of $\Gamma$, and the integrand is $O(1)$ within that neighborhood.

The final term to be considered arises because of the difference between $P(y)$ (cf. (7.3)) and its linear approximation $L(y)$ (cf. (7.4)):

[TABLE]

In the domain where $P(\epsilon j)$ and $L(\epsilon j)$ are both positive, we have

[TABLE]

This difference is zero if $P(\epsilon j)$ and $L(\epsilon j)$ are both negative. Thus,

[TABLE]

In the region where $P(Y(\alpha,0))<0$, the limit is obviously zero. Hence

[TABLE]

As before, we can apply the limit as $\epsilon\to 0$ inside the integral in (7.25) to obtain (7.28) because the integrand is uniformly bounded. Indeed,

[TABLE]

so the integrand in (7.25) remains bounded as $\epsilon\to 0$. The domain where $P(y)$ and $L(y)$ are of different signs is a shrinking $O(\epsilon)$ neighborhood of $\Gamma$, and the desired result follows.

Combining (7.19), (7.21), (7.24), and (7.28) gives the result, which, in compact form, can be written as follows

[TABLE]

Here $G_{0}:=G(Y(\alpha,0))$, $G:=G(Y(\alpha,t))$, $P:=P(Y(\alpha,t))$, $P_{0}:=P(Y(\alpha,0))$, and the derivatives with respect to $t$ are evaluated at $t=0$. This coincides with what we get by substituting $g=P_{+}G$ into the continuous inversion formula (2.14). Indeed, representing $P_{+}G=(GP)\theta(P)$, we have

[TABLE]

The coefficients in front of the delta-function in (7.30) and (7.31) match. Subtracting the coefficients in front of the Heaviside function gives:

[TABLE]

Here we have used that $G_{0}$ is independent of $t$, and the expression under the derivative has a zero of third order at $t=0$. Thus the theorem is proven in the case $x=x_{0}$.

Next, consider the case of a general $x_{\epsilon}:=x_{0}+\epsilon\tilde{x}$ (cf. (7.1)). We begin by repeating the steps (7.3)–(7.9), where all the auxiliary functions, such as $Y$, are computed using $x_{\epsilon}$ instead of $x_{0}$. It is clear that in any place where an auxiliary function is not divided by $\epsilon$, e.g. $B(\alpha)$ in (7.5) and $e(\alpha)$, $u(\alpha)$ in (7.8), replacing $x_{\epsilon}$ with $x_{0}$ introduces an error of magnitude $O(\epsilon)$. Here we also used the property (2) of $\psi$. Consequently, the analogue of (7.9) for $x_{\epsilon}$ becomes:

[TABLE]

where only $Y$ is different from the analogous function in (7.9). Note that $P(y)$ depends only on the shape of $\mathcal{S}$ in a neighborhood of $x_{0}$ and, therefore, is independent of $x_{\epsilon}$. We have

[TABLE]

for some smooth and bounded $W$. Here $Y(\alpha,0)$ is the same as in (7.9). Substituting into (7.33) gives

[TABLE]

Similarly to (7.14), introduce

[TABLE]

The point $\tilde{x}$ is fixed, so we do not need to list it in the arguments of $\psi_{2}$. Clearly, $\psi_{2}$ satisfies the same properties (1)–(3) as $\psi$. Hence Lemma (6) applies to $\psi_{2}$ as well, and we get similarly to (7.10), (7.14), and (7.15):

[TABLE]

Here we have used that the integrals with respect to $\tilde{s}_{1}$ and $v$ are unaffected by the constant (with respect to $\tilde{s}_{1}$ and $v$) shifts in (7.36). Therefore, (7.19) holds with $J_{\epsilon}^{(1a)}(x_{\epsilon})$ on the left.

To find the limit of $J_{\epsilon}^{(1b)}(x_{\epsilon})$, consider the key step in (7.21):

[TABLE]

which is rewritten with $x_{0}$ replaced by $x_{\epsilon}$. This equality holds everywhere except in an $O(\epsilon)$ neighborhood of $\Gamma$ ($=\Gamma_{x_{0}}$). Here we use that the curve $\Gamma_{x_{\epsilon}}$, which is obtained by solving $P(Y(\alpha,0;x_{\epsilon}))=0$, depends smoothly on $x_{\epsilon}$, and $\text{dist}(\Gamma_{x_{\epsilon}},\Gamma)=O(\epsilon)$ (see the argument preceding the statement of Theorem 3). Similarly to (7.21), the integrand is uniformly bounded, and we get

[TABLE]

The fact that the limits of $J_{\epsilon}^{(1c)}(x_{\epsilon})$ and $J_{\epsilon}^{(2)}(x_{\epsilon})$ as $\epsilon\to 0$ are independent of $\tilde{x}\in\tilde{\mathcal{U}}$ can be established in a similar way, and the theorem is proven. ∎

8. Numerical experiment

We start by constructing an interpolation kernel with the required properties. To obtain $\varphi$, we first obtain an interpolation kernel $\varphi_{\text{1D}}$ that has properties IK1–IK5 in $\mathbb{R}$, and then extend it to $\mathbb{R}^{3}$ in a separable fashion. To obtain $\varphi_{\text{1D}}$ we use the result of [11], where such a kernel is obtained following the method in [3]:

[TABLE]

Here $B_{n}$ is the cardinal B-spline of degree $n$ supported on $[0,n+1]$. Then the kernel $\varphi$ becomes

[TABLE]

where $\Delta_{k}$ is the data stepsize along the $k$-th axis. For simplicity, in this paper all the $\Delta_{k}$ are equal, i.e. $\Delta_{k}=\epsilon$, $k=1,2,3$.

The GRT we consider here integrates a function supported in the half-space $x_{3}>0$ over spheres that are tangent to the plane $x_{3}=0$. The family of such spheres is three-dimensional. We parametrize the spheres (and, consequently, the GRT) by the coordinates of their center $y$. Thus, the surfaces $\mathcal{S}_{y}$ are spheres, and the defining function $\Phi$ in (2.7) becomes:

[TABLE]

Clearly,

[TABLE]

The test object is the ball with center $x_{c}=(0,0,11)$, radius $R=5$, and uniform density 1. The point on the boundary $x_{0}$, in a neighborhood of which we compute resolution, is given by

[TABLE]

In agreement with our convention, $\alpha_{0}$ points into the interior of the ball.

There can be two spheres that are tangent to the ball at $x_{0}$. As an example, we consider the sphere whose center $y_{0}$ satisfies $(x_{0}-y_{0})\cdot\alpha_{0}>0$. Thus, for reconstruction near $x_{0}$ we use the data in a neighborhood of $y_{0}$. With this choice of $y_{0}$, the condition $\Phi^{\prime}_{x}(x_{0},y_{0})/|\Phi^{\prime}_{x}(x_{0},y_{0})|=-\alpha_{0}$ (see the text following (3.1)) is satisfied with $\Phi$ given by (8.3). For the selected $x_{0}$, $\alpha_{0}$, and $y_{0}$, we compute using (8.4):

[TABLE]

To compute the GRT, we use the formula for the area of the spherical cap:

[TABLE]

where $R$ is the radius of the sphere, and $h$ is the height of the cap. The values of $R$ and $h$ can be computed once the center of the sphere $\mathcal{S}_{y}$ is chosen (e.g., $R=y_{3}$). To simulate discrete data, the GRT is computed at the points $y=r+\epsilon j$. The interpolated data $g_{\epsilon}$ is computed using (2.13), where the kernel is given by (8.2) with $\Delta_{k}=\epsilon$, $k=1,2,3$.

To apply the inversion formula (2.14), we numerically integrate $g_{\epsilon}$ over a neighborhood of $\alpha_{0}$ on the unit sphere. To compute $(\partial/\partial t)^{2}g_{\epsilon}(Y(\alpha,t;x))$ at $t=0$ we use (2.13) and the chain rule as in (7.6). Given $x$, $\alpha$, and $t$, the center of the sphere containing the point $x+t\alpha$ and normal to $\alpha$ at that point (cf. the paragraph following (2.13)) is easily found to be:

[TABLE]

Consequently,

[TABLE]

and

[TABLE]

The cut-off function $\chi$ in (2.14) is constructed as follows. Let $\alpha^{\perp}$ run through the unit sphere in the plane $\alpha_{0}^{\perp}$. Then any $\alpha\in S_{+}^{2}$ ($\alpha_{0}$ is the North pole of $S_{+}^{2}$) can be represented in the form $\alpha=(\cos\omega)\alpha_{0}+(\sin\omega)\alpha^{\perp}$, $0\leq\omega\leq\pi/2$. In the code we use

[TABLE]

Finally, the predicted response is computed using (3.8). The results corresponding to $\epsilon=0.01$ are shown in Figure 1. We see a good match between the predicted and actual transition curves.