Analysis of reconstruction from discrete Radon transform data in $\mathbb R^3$ when the function has jump discontinuities

TL;DR

This paper analyzes how well discrete Radon transform data can reconstruct functions with jump discontinuities in three dimensions, providing explicit asymptotic behavior and proving accurate reconstruction at generic points.

Contribution

It provides the first explicit asymptotic analysis of reconstruction near discontinuities in 3D and establishes conditions under which the discrete inversion accurately recovers functions.

Findings

Explicit leading term behavior of reconstructed function near discontinuities

Discrete inversion formula converges at generic points with jumps as data density increases

Connection between reconstruction accuracy and uniform distribution theory

Abstract

In this paper we study reconstruction of a function $f$ from its discrete Radon transform data in $\mathbb R^3$ when $f$ has jump discontinuities. Consider a conventional parametrization of the Radon data in terms of the affine and angular variables. The step-size along the affine variable is $\epsilon$, and the density of measured directions on the unit sphere is $O(\epsilon^2)$. Let $f_\epsilon$ denote the result of reconstruction from the discrete data. Pick any generic point $x_0$ (i.e., satisfying some mild conditions), where $f$ has a jump. Our first result is an explicit leading term behavior of $f_{\epsilon}$ in an $O(\epsilon)$-neighborhood of $x_0$ as $\epsilon\to0$. A closely related question is why can we accurately reconstruct functions with discontinuities at all? This is a fundamental question, which has not been studied in the literature in dimensions three and higher. We prove that the discrete inversion formula `works', i.e. if $x_0\not\in S:=\text{singsupp}(f)$ is generic, then $f_{\epsilon}(x_0)\to f(x_0)$ as $\epsilon\to0$. The proof of this result reveals a surprising connection with the theory of uniform distribution (u.d.). This is a new phenomenon that has not been known previously. We also present some numerical experiments, which confirm the validity of the developed theory.

Equations

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

$$g(\alpha ,p)={\int }_{\alpha \cdot x=p}f(x)dx,$$

$$g({\alpha }_i,p_j),p_j=ϵ(\check{p}({\alpha }_i)+j),j\in \mathbb{Z},$$

$${\alpha }_{\vec{i}}=H(ϵ(i_1+r_1),ϵ(i_2+r_2))\text{ for any }{\alpha }_{\vec{i}}\in \Omega .$$

$$\text{det}(H^{\prime }(\vec{t}))=\text{det}\left(\frac{\partial \alpha }{\partial (t_1,t_2)}\right)>0$$

$$\sum _{j\in \mathbb{Z}}{j}^m\phi (t-j)=t^m,0\le m\le 2,t\in \mathbb{R},$$

$$g_{ϵ}({\alpha }_i,p):=\sum _{j\in \mathbb{Z}}g({\alpha }_i,ϵ(\check{p}({\alpha }_i)+j))\phi \left(\frac{p-ϵ(\check{p}({\alpha }_i)+j)}{ϵ}\right).$$

$$f(x)=-\frac{1}{4{\pi }^2}{\int }_{S_{+}^2}{(\partial /\partial p{)}^{2}g(\alpha ,p)|}_{p=\alpha \cdot x}d\alpha ,$$

$$f_{ϵ}(x)=-\frac{1}{4{\pi }^2}\sum _{i:{\alpha }_i\in S_{+}^2}{c}_i{(\partial /\partial p{)}^2{g}_{ϵ}({\alpha }_i,p)|}_{p={\alpha }_i\cdot x},$$

$$a{ϵ}^2\le c_i\le b{ϵ}^2$$

$$x_h:=x_0+ϵh{\Theta }_0,$$

$$f_0:=\underset{ϵ\to 0^{+}}{\lim }f(x_0+ϵ{\Theta }_0),{\rho }_0:=\underset{ϵ\to 0^{+}}{\lim }(f(x_0+ϵ{\Theta }_0)-f(x_0-ϵ{\Theta }_0)).$$

$$\alpha ={\alpha }^{\perp }+{\alpha }_3{\Theta }_0\in S^2,{\alpha }^{\perp }=({\alpha }_1,{\alpha }_2)\in {\Theta }_0^{\perp },{\alpha }_3=\sqrt{1-|{\alpha }^{\perp }{|}^2}.$$

$$\begin{array}{r}q(\alpha ):=\alpha \cdot x_0,v(\alpha ):=\alpha \cdot x_0-p_0(\alpha ).\end{array}$$

$$\underset{ϵ\to 0}{\lim }{f}_{ϵ}(x_h)=f_0-{\rho }_0{\int }_h^{\infty }\phi (s)ds.$$

$$\begin{array}{rl}{\Omega }_1& :=\{\alpha \in S_{+}^2:|{\alpha }^{\perp }|<A\sqrt{ϵ}\},{\Omega }_2:=\{\alpha \in S_{+}^2:A\sqrt{ϵ}<|{\alpha }^{\perp }|<\omega \},\\ {\Omega }_3& :=S_{+}^2\setminus ({\Omega }_1\cup {\Omega }_2),\end{array}$$

$$g(\alpha ,p)=(p-p_0(\alpha ){)}_{+}G(\alpha ,p-p_0(\alpha )),$$

$$G(\alpha ,0)=2\pi {\rho }_0(\alpha )/\sqrt{\text{det}Q(\alpha )},$$

$$p_0(\alpha )=\alpha \cdot x_0-\frac{Q^{-1}({\Theta }_0){\alpha }^{\perp }\cdot {\alpha }^{\perp }}{2}+O({ϵ}^{3/2})$$

$$\begin{array}{rl}{g}_{ϵ}({\alpha }_i,p)=& \sum _{j\in \mathbb{Z}}\left[\frac{2\pi {\rho }_0({\alpha }_i)}{\sqrt{\text{det}Q({\alpha }_i)}}{\left(ϵ(\check{p}+j)-p_0({\alpha }_i)\right)}_{+}+O\left((ϵ(\check{p}+j)-p_0({\alpha }_i){)}_{+}^2\right)\right]\\ & \times \phi \left(\frac{p-ϵ(\check{p}+j)}{ϵ}\right).\end{array}$$

$$\psi (q,s):=\sum _{j\in \mathbb{Z}}(j-q+s{)}_{+}{\phi }^{\prime \prime }(q-j),q,s\in \mathbb{R}.$$

$$\text{the derivatives ψq′ and ψs′ are piecewise-continuous and bounded},$$

$$\psi (q,s)=\psi (q+n,s)\text{ for any }n\in \mathbb{Z},$$

$$\psi (q,s)=0\text{ if ∣s∣ is large enough.}$$

$$\begin{split}&f_{\epsilon}^{(1)}(x_{h})\\ &=-\frac{1}{2\pi}\sum_{|\alpha_{i}^{\perp}|<A\sqrt{\epsilon}}c_{i}\sum_{j\in\mathbb{Z}}\\ &\qquad\biggl{[}\frac{\rho_{0}(\alpha_{i})}{\sqrt{\text{det}Q(\alpha_{i})}}\left(\epsilon(\check{p}+j)-\alpha_{i}\cdot x_{0}+\frac{Q^{-1}\alpha_{i}^{\perp}\cdot\alpha_{i}^{\perp}}{2}+O(\epsilon^{3/2})\right)_{+}+O(\epsilon^{2})\biggr{]}\\ &\qquad\times\epsilon^{-2}\varphi^{\prime\prime}\left(\frac{\alpha_{i}\cdot x_{0}-\check{p}}{\epsilon}+h+O(\epsilon)-j\right)\\ &=-\frac{\rho_{0}}{2\pi\sqrt{\text{det}Q}}\sum_{|\alpha_{i}^{\perp}|<A\sqrt{\epsilon}}\frac{c_{i}}{\epsilon}\psi\left(\frac{\alpha_{i}\cdot x_{0}}{\epsilon}-\check{p}_{0}+h,h+\frac{Q^{-1}\alpha_{i}^{\perp}\cdot\alpha_{i}^{\perp}}{2\epsilon}\right)+O(\epsilon^{1/2}),\end{split}$$

$${\rho }_0:={\rho }_0({\Theta }_0),{\check{p}}_0:=\check{p}({\Theta }_0),Q:=Q({\Theta }_0).$$

$${\tilde{c}}_i=\frac{c_i}{ϵ}\left(1+O({ϵ}^{1/2})\right).$$

$$\tilde{q}(\vec{t}):=H(\vec{t})\cdot x_0,{\mathcal{H}}_0^{\prime }:={\tilde{q}}^{\prime }(\vec{0}),{\mathcal{H}}_0^{\prime \prime }:={\tilde{q}}^{\prime \prime }(\vec{0}).$$

$$\frac{{\alpha }_{\vec{i}}\cdot x_0}{ϵ}=\frac{{\Theta }_0\cdot x_0}{ϵ}+(i_1+r_1,i_2+r_2)\cdot {\mathcal{H}}_0^{\prime }+\frac{{\mathcal{H}}_0^{\prime \prime }{\tilde{\alpha }}_{\vec{i}}^{\perp }\cdot {\tilde{\alpha }}_{\vec{i}}^{\perp }}{2}+O({ϵ}^{1/2}).$$

$$\begin{array}{rl}{f}_{ϵ}^{(1)}(x_h)=& -\frac{{\rho }_0}{2\pi \sqrt{\text{det}Q}}{J}_{ϵ}+O({ϵ}^{1/2}),\\ J_{ϵ}:=& \sum _{|{\tilde{\alpha }}_{\vec{i}}^{\perp }|<A}{\tilde{c}}_{\vec{i}}\psi \left(A_{ϵ}+\vec{i}\cdot {\mathcal{H}}_0^{\prime }+\frac{{\mathcal{H}}_0^{\prime \prime }{\tilde{\alpha }}_{\vec{i}}^{\perp }\cdot {\tilde{\alpha }}_{\vec{i}}^{\perp }}{2},h+\frac{Q^{-1}{\tilde{\alpha }}_{\vec{i}}^{\perp }\cdot {\tilde{\alpha }}_{\vec{i}}^{\perp }}{2}\right),\end{array}$$

Full text

A closely related question is why can we accurately reconstruct functions with discontinuities at all? This is a fundamental question, which has not been studied in the literature in dimensions three and higher. We prove that the discrete inversion formula “works”, i.e. if $x_{0}\not\in S:=\text{singsupp}(f)$ is generic, then $f_{\epsilon}(x_{0})\to f(x_{0})$ as $\epsilon\to 0$. The proof of this result reveals a surprising connection with the theory of uniform distribution (u.d.). This is a new phenomenon that has not been known previously. We also present some numerical experiments, which confirm the validity of the developed theory.

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

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

1. Introduction

In this paper we consider the question: “Why and how well does tomographic reconstruction from discrete data work?” The answer is known if the function $f$ to be reconstructed is sufficiently smooth. Related problems are usually studied by sampling theory, where the goal is to estimate how closely the reconstructed function $f_{\epsilon}$ approximates $f$ in some global norm [9, 3, 11, 12, 15, 2, 5, 4, 20, 7]. Usually, $f$ is required to be essentially bandlimited, which imposes a smoothness requirement on $f$. In the case when $f$ is not smooth (e.g. has jump discontinuities) the question has not been studied much. It consists of two parts:

• Q1.

What does reconstruction look like near the singularities of $f$?

• Q2.

What is the effect of “remote” singularities, i.e. the part of $\text{singsupp}(f)$ located at a distance from the reconstruction point?

Frequently, one is less interested in how $f_{\epsilon}$ approximates $f$ in some global norm. Instead, one would like to know how accurately and with what resolution the singularities of $f$ are reconstructed. Convergence of reconstruction algorithms in the case of objects with discontinuities has been studied as well [6, 16, 14, 17]. However, in these works the discontinuities of the object are a complicating factor rather than the object of the study. The first paper focusing specifically on the behavior of $f_{\epsilon}$ near a jump discontinuity of $f$ is [8], where the author considers inversion of the Radon transform in $\mathbb{R}^{2}$ in the parallel beam setting. The object can be static or change with time. The parametrization of the data is conventional, 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 descrete data. Pick a point $x_{0}$, where $f$ has a jump. We suppose that $x_{0}$ is generic, i.e. it satisfies some mild conditions. The result of [8] is an explicit leading term behavior of $f_{\epsilon}$ in an $O(\epsilon)$-neighborhood of $x_{0}$ as $\epsilon\to 0$. The obtained behavior, which we call transition behavior or, equivalently, edge response provides the desired resolution of the reconstruction algorithm. No such results were known in dimensions greater than two. In this paper we obtain the edge response in $\mathbb{R}^{3}$.

A closely related question is why can we accurately reconstruct functions with discontinuities at all? This is a fundamental question, which has not been studied in the literature in dimensions three and higher. The answer to the question is complicated. In $\mathbb{R}^{2}$, convergence of reconstruction algorithms in the case of objects with discontinuities has been studied in [6, 13, 16, 17]. In this paper we answer the question in $\mathbb{R}^{3}$. The key reason why reconstruction works comes from a surprising connection with the theory of uniform distribution (u.d.). This is a completely new phenomenon, which has not been noticed previously. The first hint of such a connection appeared in [8], where the transition behavior is derived using the u.d. property of a relevant sequence.

Let us discuss the main results of this paper in more detail. Similarly to the 2D case, we assume here that the step-size along the affine variable is $\epsilon$, and the density of measured directions on the unit sphere is $O(\epsilon^{2})$. First, we extend the notion of a generic point. Second, we obtain the transition behavior of $f_{\epsilon}$ in an $O(\epsilon)$-neighborhood of a generic point $x_{0}$ where $f$ has a jump discontinuity. This answers the question Q1 above. A connection with the u.d. theory here is similar to the one in [8].

Third, we answer the question Q2. We show that if $x_{0}\not\in S:=\text{singsupp}(f)$ is generic, then $f_{\epsilon}(x_{0})\to f(x_{0})$ as $\epsilon\to 0$. This means that if $x_{0}\in S$, then the singularities of $f$ at a distance from $x_{0}$ (i.e., “remote” singularities) do not contribute to the transition behavior at $x_{0}$. Only the behavior of $f$ near $x_{0}$ contributes to the edge response. There is no contradiction between the assumptions $x_{0}\in S$ and $x_{0}\not\in S$, since the singularities of $f$ near $x_{0}$ and away from $x_{0}$ can be separated by a partition of unity using the linearity of the Radon transform. Additionally, this result proves that the discrete inversion formula “works” (i.e., provides increasingly accurate results as $\epsilon\to 0$) and illuminates the reasons why. We show that reconstruction is guaranteed to work if certain collection of points is u.d. inside a set, whose shape depends on $S$. A major complication in the proof is that the set shrinks as $\epsilon\to 0$. This connection with the u.d. theory never appeared in the literature before (including [8]).

As was noted in [13, 14], the behavior of the reconstruction algorithm depends on the geometry of $S$. Clearly, it also depends on the strength of singularities of $f$ (e.g., in the Sobolev scale). Comprehensive analysis of all these cases is beyond the scope of this paper. Here we consider only the case when $f$ has jump discontinuities across $S$, and $S$ is the union of smooth convex surfaces.

The paper is organized as follows. In Section 2 we describe the problem set-up, specify the class of functions we consider, describe discrete Radon transform data, introduce the discrete inversion formula, define the notion of a generic point, and formulate the main result (Theorem 1). In Section 3 we obtain the edge response in a neighborhood of a generic point $x_{0}\in S$. In Section 4 we show that “remote” singularities do not contribute to the edge response. In Section 5 we illustrate on numerical experiments that when $x_{0}$ is generic, then the theoretically predicted and numerically computed edge responses match very well. We also show that if $x_{0}$ is not generic, then the match may no longer be accurate. Technical results related to the u.d. property of certain sets of points, which are needed in Section 4, are proven in Appendix A.

2. Preliminary construction

Consider a function $f(x)$, which can be represented as a finite sum

[TABLE]

where $\chi_{D_{j}}$ is the characteristic function of the domain $D_{j}$. We assume that for each $j$:

1. (1)

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

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

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

The Radon transform of $f$ is defined as follows:

[TABLE]

where $dx$ is the area element on the plane $\Pi(\alpha,p):=\{x\in\mathbb{R}^{3}:\alpha\cdot x=p\}$. The discrete data are given by

[TABLE]

where $\check{p}\in C^{1}(S^{2})$, and $i\in\mathbb{N}$, $i\leq O(\epsilon^{-2})$, is the index that enumerates the measured directions $\alpha_{i}\in S^{2}$. When considering $\alpha_{i}$ in some small open set $\Omega\subset S^{2}$, we will use a 2D multi-index $\vec{i}=(i_{1},i_{2})$ instead of $i$: $\alpha_{\vec{i}}$. The assumption is that for almost all $\alpha_{0}\in S^{2}$ there exists an open set $\Omega$, $\alpha_{0}\in\Omega\subset S^{2}$, such that

[TABLE]

Here $r_{1,2}$ are some constants, $0\leq r_{1,2}<1$, which may depend on $\epsilon$. The function $H(\vec{t}):U\to\Omega$, where $U\subset\mathbb{R}^{2}$ is some bounded domain, is a smooth diffeomorphism. Both $H$ and $U$ may depend on $\alpha_{0}$ and $\Omega$. We also assume that the determinant

[TABLE]

is bounded away from zero in each $U$.

Let $\varphi$ be an interpolating kernel, which satisfies the following assumptions:

• A1.

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

[TABLE]

• A2.

$\varphi$ is compactly supported,

• A3.

The derivatives $\varphi^{(m)}$, $1\leq m\leq 3$, exist,

• A4.

$\varphi^{(3)}$ is piecewise continuous and bounded,

• A5.

$\int\varphi(t)dt=1$.

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

[TABLE]

To simplify notations, the dependence of $\check{p}$ on $\alpha$ is omitted in most places.

Since the Radon transform is even, the classical inversion formula with continuous data is given by

[TABLE]

where $S_{+}^{2}$ is any half of the unit sphere in $\mathbb{R}^{3}$. The discrete inversion formula is given by

[TABLE]

where $c_{i}$ are integration weights. We assume that there exist $0<a\leq b$ so that

[TABLE]

for all $i$ as $\epsilon\to 0$. A convenient and useful way to think of the coefficients $c_{i}$ is that they represent the area of pieces that tessellate $S^{2}$, each piece is sufficiently regular and contains only one $\alpha_{i}$. In particular, $c_{i}=|H^{\prime}(\vec{t})|\epsilon^{2}+O(\epsilon^{3})$ for any $\vec{t}=(t_{1},t_{2})$ such that $\alpha=H(\vec{t})$ belongs to the same tessellation piece as $\alpha_{i}$.

Pick a point $x_{0}\in S:=\text{singsupp}f$, and suppose that $S$ is a smooth surface with positive principal curvatures in a neighborhood of $x_{0}$. Let $\Theta_{0}$ be the unit vector normal to $S$ at $x_{0}$ and pointing into the interior of the corresponding domain $D_{j}$. Consider the point

[TABLE]

where $h$ varies over a bounded set. Denote also:

[TABLE]

Introduce the coordinates $\alpha^{\perp}$ on $S^{2}$ in a neighborhood of $\Theta_{0}$:

[TABLE]

Let $p=p_{0}(\alpha):\,S^{2}\to\mathbb{R}$ be a function determined by the condition that the plane $\Pi(\alpha,p_{0}(\alpha))$ be tangent to $S$. Generally, $p_{0}$ can be multi-valued. In what follows we always consider one of its single-valued local branches. Pick some $x_{0}$ and introduce the functions

[TABLE]

Pick any $\alpha_{0}\in S^{2}$ such that $x_{0}\in\Pi(\alpha_{0},p_{0}(\alpha_{0}))$. If $S$ has positive principal curvatures at a point $z_{0}$, where $\Pi(\alpha_{0},p_{0}(\alpha_{0}))$ is tangent to $S$, and $x_{0}\not=z_{0}$, then the equation $v(\alpha)=0$ determines a locally smooth curve, which we denote by $\Gamma$. By construction, $\alpha_{0}\in\Gamma$. The curve is smooth, because $v^{\prime}(\alpha_{0})=x_{0}^{\perp}-z_{0}^{\perp}\not=0$. The derivative with respect to $\alpha$ is computed on the unit sphere, and $x_{0}^{\perp},z_{0}^{\perp}$ are the orthogonal projections of $x_{0},z_{0}$, respectively, onto the plane $\alpha_{0}^{\perp}$. Clearly, $\Pi(\alpha_{0},p_{0}(\alpha_{0}))$ can be tangent to $S$ at other points, but here $p_{0}(\alpha)$ is the local branch that is valid in a neighborhood of $z_{0}$. Our statements apply to any of the branches and the corresponding points of tangency.

Definition 1**.**

We say that a point $x_{0}\in S=\text{singsupp}f$ is generic if

1. (1)

$S$* is a smooth surface with positive principal curvatures in a neighborhood of $x_{0}$. In particular, $S$ does not self-intersect at $x_{0}$;* 2. (2)

Let $\Theta_{0}$ be the unit vector normal to $S$ at $x_{0}$, and let $\vec{t}^{*}$ be such that $\Theta_{0}=H(t^{*})$. The vector $\nabla(H(\vec{t})\cdot x_{0})|_{\vec{t}=\vec{t}^{*}}$ has at least one irrational component; and 3. (3)

If a plane $\Pi$ contains $x_{0}$ and is tangent to $S$ at some $z_{0}\not=x_{0}$, then the principal curvatures of $S$ at $z_{0}$ are positive.

An additional condition for a point to be generic is formulated in Section 4.

Definition 2**.**

We say that a point $x_{0}\not\in S=\text{singsupp}f$ is generic if it satisfies condition (3) above and the additional condition in Section 4.

Now we formulate the main result of the paper.

Theorem 1**.**

For a generic $x_{0}\in S$, one has

[TABLE]

The proof of the theorem is split into two sections.

3. Contribution of the local singularity

In this section we choose the hemisphere $S_{+}^{2}:=\{\alpha\in S^{2}:\ \alpha_{3}>0\}$ in (2.8) and (2.9). The range of $\alpha$ in (2.9) can be split into three sets:

[TABLE]

for some small (but fixed) $\omega>0$. Here $A>0$ is a large parameter. Let $f_{\epsilon}^{(j)}$ denote the value of the sum in (2.9) ranging over indices $i$ such that $\alpha_{i}\in\Omega_{j}$, $j=1,2,3$. Consider $f_{\epsilon}^{(1)}$ first.

It will be shown below (see Section 3.1) that the Radon transform of locally smooth components of $f$ does not contribute to the transition behavior of $f_{\epsilon}$. Hence, without loss of generality, we can assume that $f\equiv 0$ in a neighborhood of $x_{0}$ on the exterior side of $S$ (this is the side corresponding to $h<0$ in (2.11)). By linearity and using a partition of unity if necessary, we can assume that $\text{supp}(f)$ is contained in a small neighborhood of $x_{0}$. The local branch of $p_{0}$ that we use in this section is determined by the condition that the point of tangency be close to $x_{0}$. As is known, [18, 19], in this case we have

[TABLE]

where $G$ is a smooth function near $(\Theta_{0},0)$,

[TABLE]

and $Q(\alpha)$ is the matrix of the second fundamental form of $S$ in the basis $(\alpha_{1},\alpha_{2})$ at the point where $\Pi(\alpha,p_{0}(\alpha))$ is tangent to $S$. A simple calculation shows that

[TABLE]

when $|\alpha^{\perp}|=O(\sqrt{\epsilon})$. Hence the interpolated data are given by

[TABLE]

In view of (2.9) and (3.5), denote

[TABLE]

The assumptions A1–A4 imply that

[TABLE]

[TABLE]

[TABLE]

Combining (3.4)–(3.6) and using (2.9), (3.1), we obtain

[TABLE]

where we have introduced the notation

[TABLE]

The $O(\epsilon^{2})$ term in brackets in (3.10) corresponds to the second term in brackets in (3.5). This follows because $|x_{h}-x_{0}|=O(\epsilon)$, $\varphi$ is compactly supported, and $|\alpha_{i}^{\perp}|=O(\epsilon^{1/2})$ (cf. (3.1) and (3.4)). In (3.10) we used that (a) there are $O(1/\epsilon)$ terms in the sum with respect to $i$, and (b) using (2.10), the sum in (3.10) remains bounded as $\epsilon\to 0$.

Set $\tilde{\alpha}^{\perp}:=\alpha^{\perp}/\sqrt{\epsilon}$. The tessellation of $S^{2}$ based on $\alpha_{i}$ induces a tessellation of the domain $|\tilde{\alpha}^{\perp}|<A$. Let $\tilde{c}_{i}$ denote the area of the tessellation piece corresponding to $\tilde{\alpha}_{i}^{\perp}$. From (2.13),

[TABLE]

The term $O(\epsilon^{1/2})$ here also includes the Jacobian of the transformation $\mathbb{R}^{2}\ni\alpha^{\perp}\to\alpha\in S^{2}$.

Apply (2.4) to an open set containing $\Theta_{0}$. By shifting indices $i_{1,2}$ and adjusting $r_{1,2}\in[0,1)$, if necessary, we can always assume that $H(0,0)=\Theta_{0}$. Denote

[TABLE]

Then

[TABLE]

Combining the above, rewrite (3.10) as follows

[TABLE]

where $A_{\epsilon}$ depends on $\epsilon$, but is independent of $\vec{i}$. The $O(\epsilon^{1/2})$ term in (3.14) is absorbed by the $O(\epsilon^{1/2})$ term in (3.15).

To find the limit of $J_{\epsilon}$ as $\epsilon\to 0$, pick some small $\epsilon^{\prime}>0$ and break up the domain $|\tilde{\alpha}_{\vec{i}}^{\perp}|<A$ into squares $B_{l}$ of size $\epsilon^{\prime}$. Consider one $B_{l}$ and all $\tilde{\alpha}_{\vec{i}}^{\perp}$ within it. Replace all $\tilde{\alpha}_{\vec{i}}^{\perp}\in B_{l}$ with $\tilde{\alpha}_{*}^{\perp}$, where $\tilde{\alpha}_{*}^{\perp}$ is the center of $B_{l}$, everywhere in the arguments of $\psi$. This leads to an error of magnitude $O(\epsilon^{\prime})$. By $O(\epsilon^{\prime})$ we denote any quantity whose magnitude does not exceed $c\epsilon^{\prime}$, where $c$ is independent of the selected square $B_{l}$, the chosen $\tilde{\alpha}_{\vec{i}}^{\perp}\in B_{l}$, and $\epsilon$. Clearly,

[TABLE]

By (3.8) and (3.16), we can consider the set

[TABLE]

where ${\vec{i}}$ are such that $\tilde{\alpha}_{\vec{i}}^{\perp}\in B_{l}$. There are $O(1/\epsilon)$ such $\vec{i}$: $O(\epsilon^{-1/2})$ per each coordinate $i_{1,2}$. Here and in what follows, $\{x\}$ denotes the fractional part of a number $x\in\mathbb{R}$. By Theorem 2.9 and Example 2.9 of [10], the set $\left\{\vec{i}\cdot{\mathcal{H}}_{0}^{\prime}\right\}$ is u.d. as $\epsilon\to 0$ if at least one of the components of the vector ${\mathcal{H}}_{0}^{\prime}$ is irrational (see condition (2) in the definition of a generic point). Clearly, the shift by $A_{\epsilon}^{\prime}$ does not affect the u.d. property. Therefore, for the sum of the terms in (3.15) corresponding to $\tilde{\alpha}_{i}^{\perp}\in B_{l}$ we obtain

[TABLE]

Here $N_{\epsilon}(B_{l})$ is the number of terms in the sum in (3.18), and $\tilde{c}_{*}$ is the average of $\tilde{c}_{\vec{i}}$. By construction, $|\tilde{c}_{\vec{i}}-\tilde{c}_{*}|/\epsilon=O(\epsilon^{\prime})$ whenever $\tilde{\alpha}_{\vec{i}}^{\perp}\in B_{l}$.

Sum over all squares and observe that the resulting Riemann sum is within $O(\epsilon^{\prime})$ from the corresponding integral. This yields:

[TABLE]

Since $\epsilon^{\prime}>0$ can be as small as we like, combining (3.19) with (3.15) gives

[TABLE]

In view of (3.9), if $h$ is fixed, the integrand in (3.20) is zero when $|\tilde{\alpha}^{\perp}|$ is large enough (recall that $Q(\Theta_{0})$ is positive definite). Hence we can select $A>0$ large enough and obtain by changing variables

[TABLE]

Lemma 1**.**

One has

[TABLE]

for some $c>0$ large enough.

Proof.

By (3.6) and (3.9),

[TABLE]

Using A2 and A5 we finish the proof. ∎

3.1. Analysis of the terms $f_{\epsilon}^{(2)}$ and $f_{\epsilon}^{(3)}$.

The Radon transform of $f$ is smooth when $(\alpha,p)$ are such that $\alpha\in\Omega_{3}$ and $|p-\alpha\cdot x_{0}|=O(\epsilon)$, i.e. $\Pi(\alpha,p)$ is not tangent to $S$. Hence the interpolated data has derivatives with respect to $p$ up to the order $3$ uniformly bounded as $\epsilon\to 0$. Therefore,

[TABLE]

For a function $g\in C^{3}(\mathbb{R})$ we have

[TABLE]

where we have used (2.6). Combining with (3.24) this implies

[TABLE]

In a similar fashion, we will show that

[TABLE]

Choose $w>0$ in (3.1) so that $\alpha\cdot x_{0}-p_{0}(\alpha)$ is strictly convex on $\Omega_{1}\cup\Omega_{2}$. In view of (3.4), this is possible. Hence the lower bound on the values of $\alpha\cdot x_{0}-p_{0}(\alpha)$ on $\Omega_{2}$ is determined by looking at its values at the boundary of $\Omega_{1}$, where (3.4) holds. Pick $c>0$ large enough. Our argument implies that there exist sufficiently large $A>0$ and sufficiently small $\epsilon_{0}>0$ such that $\alpha\cdot x_{0}-p_{0}(\alpha)>c\epsilon$ on $\Omega_{2}$ for all $\epsilon$, $0<\epsilon\leq\epsilon_{0}$. Since $\varphi$ is compactly supported, $g(\alpha,p)$ is a smooth function in a neighborhood of all $(\alpha_{i},p_{j})$ such that $\alpha_{i}\in\Omega_{2}$ and $(\alpha_{i}\cdot x_{h}-p_{j})/\epsilon\in\text{supp}(\varphi)$. The reason is that for such $(\alpha,p)$ we can drop the subscript ${}^{\prime}+^{\prime}$ from $(p-p_{0})_{+}$ in (3.2):

[TABLE]

Hence (3.27) follows similarly to (3.24)–(3.26).

Examining the above argument more carefully, we see that on $\Omega_{1}\cup\Omega_{2}$, the part of the Radon transform of $f$ that is used in the inversion formula coincides with a globally smooth function, so we have

[TABLE]

where $f_{0}$ is defined in (2.12). The notation “$+0$” means that the derivative with respect to $p$ is evaluated on the interior side of the domain where the Radon transform of $f$ is smooth. The last equality in (3.29) is an immediate corollary of the inversion formula and the smoothness of the Radon transform.

4. Contribution of remote singularities

Suppose the plane $\Pi(\Theta_{0},p_{0}(\Theta_{0}))$ is tangent to $S$ at some $z_{0}$, $z_{0}\not=x_{0}$. In this section $x_{0}$ is generic, but otherwise arbitrary, and $\Theta_{0}\in S^{2}$ is arbitrary as well. In particular, $x_{0}$ is not necessarily on $S$. Again, using linearity, we suppose that $\text{supp}(f)$ is a subset of a small neighborhood of $z_{0}$. The surface $S$ can have a wide variety of shapes in a neighborhood of $z_{0}$. Some of these shapes may produce artifacts (e.g., if $S$ is locally flat near $z_{0}$), while others - will not. Comprehensive analysis of all such cases is beyond the scope of this paper. Here we will consider one generic situation: the principal curvatures of $S$ at $z_{0}$ are positive. In this section we prove that the singularity at $z_{0}$ does not contribute to the transition behavior of the reconstruction at $x_{0}$. More precisely, we will prove that in the limit $\epsilon\to 0$ the discretized reconstruction formula (2.9) gives the same result as the exact reconstruction formula (2.8).

To illustrate the main ideas, we start with the case when the reconstruction point is $x_{0}$ rather than $x_{h}$. As noted in Section 2, assumption (3) in the definition of a generic point implies that the equation $\alpha\cdot x_{0}=p_{0}(\alpha)$ determines a locally smooth curve denoted $\Gamma$. Let $\chi(\alpha)$ be a smooth cut-off function with small support such that $\chi\not\equiv 0$ on $\Gamma$. Using (2.14), introduce the functions

[TABLE]

where $U$ is the domain of $H$. The function $\tilde{q}$ here is the same as in (3.13).

Let $U_{t}$ and $\Gamma_{t}$ be the images of $\text{supp}(\chi)$ and $\Gamma\cap\text{supp}(\chi)$, respectively, under the map $H^{-1}:\,\alpha\to\vec{t}$. Using (2.5) and that $v^{\prime}(\alpha)\not=0$ for any $\alpha\in\Gamma_{t}$, we can assume that $\text{supp}(\chi)$ is so small that either $\tilde{v}_{1}^{\prime}(\vec{t})\not=0$ for any $\vec{t}\in\Gamma_{t}$ or $\tilde{v}_{2}^{\prime}(\vec{t})\not=0$ for any $\vec{t}\in\Gamma_{t}$. Without loss of generality, suppose $\tilde{v}_{2}^{\prime}(\vec{t})\not=0$. This implies that an equation of $\Gamma_{t}$ can be written in the form

[TABLE]

for some smooth $A$. Indeed, setting $\tilde{v}(t_{1},A(t_{1}))\equiv 0$, the Implicit Function theorem implies that $A(t_{1})$ is well-defined, and $A^{\prime}=-\tilde{v}_{1}^{\prime}/\tilde{v}_{2}^{\prime}$.

Now we can formulate the last condition that a generic point satisfies.

1. (4)

Pick any local piece $\Gamma_{t}$ constructed as above and any $M_{1},M_{2}\in\mathbb{Z}$ such that $M_{1}^{2}+M_{2}^{2}>0$. Consider the function $f(t_{1})=M_{1}\tilde{q}(t_{1},A(t_{1}))+M_{2}A(t_{1})$ defined on the domain $\{t_{1}\in\mathbb{R}:\,(t_{1},A(t_{1}))\in\Gamma_{t}\}$. We have

1. (a)

The set of $t_{1}$ such that $f^{\prime\prime}(t_{1})=0$ is the union of a finite number of distinct points and non-intersecting intervals, and 2. (b)

If $f^{\prime\prime}(t_{1})\equiv 0$ on an interval, then $f^{\prime}(t_{1})$ is irrational on that interval.

An obvious modification needs to be made for pieces where $\tilde{v}_{1}^{\prime}(\vec{t})\not=0$, $t\in\Gamma_{t}$.

Using (3.2), (3.3), define

[TABLE]

Let $f_{\chi,\epsilon}^{(j)}$ denote the result of inserting $\chi$ and $g^{(j)}$ into the discrete inversion formula (2.9), $j=1,2$. By linearity, $f_{\chi,\epsilon}=f_{\chi,\epsilon}^{(1)}+f_{\chi,\epsilon}^{(2)}$. Similarly to (3.5), (3.10):

[TABLE]

Let $[t_{1}^{(\text{min})},t_{1}^{(\text{max})}]$ be the projection of $\Gamma_{t}$ onto the $t_{1}$-axis. Partition the projection onto the intervals of length $\Delta t_{1}$: $[t_{1}^{(m)},t_{1}^{(m+1)}]$, $t_{1}^{(m+1)}-t_{1}^{(m)}=\Delta t_{1}$. Pick one of the subintervals, e.g. the $m$-th. Let $t_{1}^{*}$ be its center, and set $\vec{t}^{*}:=(t_{1}^{*},A(t_{1}^{*}))$, $\alpha^{*}:=H(\vec{t}^{*})$. For convenience, denote $\bar{i}_{j}:=i_{j}+r_{j}$, $j=1,2$, and let $\mathcal{I}_{m}$ be the set of indices $\vec{i}$ such that $\epsilon\bar{i}_{1}\in[t_{1}^{(m)},t_{1}^{(m+1)}]$. Rewrite the part of the sum in (4.4) that corresponds to the indices $\vec{i}\in\mathcal{I}_{m}$:

[TABLE]

Here we have used that the sum is $O(\Delta t_{1})$ uniformly as $\epsilon\to 0$, and that $c_{\vec{i}}=c^{*}(1+O(\Delta t_{1}))$ for all $\vec{i}$ such that $\vec{i}\in\mathcal{I}_{m}$ and $\psi(\cdot,\cdot)\not=0$, where $c^{*}$ is the area of the tessellation piece containing $\alpha^{*}$. By (3.9), there are $O(\epsilon\Delta t_{1})/\epsilon^{2}$ nonzero terms in the sum. Combined with (2.10), we get that the sum is indeed of order $O(\Delta t_{1})$.

By construction,

[TABLE]

The assumption $\tilde{v}_{2}^{\prime}(\vec{t})\not=0$, $\vec{t}\in\Gamma_{t}$, implies that the lines $t_{1}=\epsilon\bar{i}_{1},t_{2}\in\mathbb{R}$, $\epsilon\bar{i}_{1}\in[t_{1}^{(m)},t_{1}^{(m+1)}]$ intersect $\Gamma_{t}$ transversely. Then

[TABLE]

where we have used that $|\bar{i}_{2}-(A(\epsilon\bar{i}_{1})/\epsilon)|$ is bounded whenever $\psi\not=0$, and $|\epsilon\bar{i}_{1}-t_{1}^{*}|=O(\Delta t_{1})$. There are two big-$O$ terms in (4.7), because $\epsilon$ and $\Delta t_{1}$ are independent variables. Similarly,

[TABLE]

By (3.8), the non-zero terms in (4.6) can be written in the form

[TABLE]

where we denoted

[TABLE]

Because of (3.9), there are finitely many values of $\bar{j}_{2}$ for which $\psi$ in (4.9) is non-zero.

For a fixed value of $\bar{j}_{2}$, the effect of the term $-\check{p}(\alpha^{*})+c_{q}^{*}\bar{j}_{2}$ in the first argument of $\psi_{1}$ is just a shift that does not affect the u.d. property. So it can be ignored if we show that the points

[TABLE]

are u.d. in the square $[0,1]\times[-1,0]$ as $\epsilon\to 0$. Corollary 2 and Assumption (4) imply that the set

[TABLE]

is u.d. in $[0,1]\times[-1,0]$ as $\epsilon\to 0$. Then, Lemma 3 implies the desired result.

Pick one value of $\bar{j}_{2}$. The arguments of $\psi_{1}$ on the last line in (4.9) are u.d. over the rectangle $[0,1]\times[\bar{j}_{2}-1,\bar{j}_{2}]$ with area $1$. There are $N_{1}(\epsilon):=\Delta t_{1}/\epsilon$ distinct values of $i_{1}$ such that $\epsilon\bar{i}_{1}\in[t_{1}^{(m)},t_{1}^{(m+1)}]$. Thus, the area per each value is $1/N_{1}(\epsilon)$. Summation over $\bar{j}_{2}$ just shifts this rectangle along the second variable. These rectangles are stacked precisely one on top of the other without overlap. Hence

[TABLE]

Using (3.23) and then A5 gives:

[TABLE]

From (4.5), (4.6), (4.9), (4.13), and (4.14) it follows that we have to evaluate the following expression:

[TABLE]

Here we have used that the elementary area in the $t_{1},t_{2}$-coordinates is $\epsilon^{2}$, and the area of the corresponding piece on $S^{2}$ is $c^{*}$. Combining the above results gives

[TABLE]

Next we sum $J_{m}(\epsilon)$ for all $m$. Since $\Delta t_{1}>0$ can be arbitrarily small, we conclude from (4.16) by using (4.7) that

[TABLE]

The result in (4.17) is, of course, expected, since this is what one gets in the continuous case by substituting $\chi(\alpha)$ and $g^{(1)}(\alpha,p)$ (cf. (4.3)) into (2.8). The same result for $g^{(2)}(\alpha,p)$, which leads to $f_{\chi,\epsilon}^{(2)}$, follows immediately, since the second order derivative or $g^{(2)}(\alpha,p)$ with respect to $p$ is piecewise-smooth and bounded on $S^{2}$. In this case the contribution of the terms in an $\epsilon$-neighborhood of $\Gamma$ goes to zero as $\epsilon\to 0$.

Next we discuss how the above derivation changes when $x_{h}$ replaces $x_{0}$ in (4.4). Clearly, it is sufficient to consider only the leading singular term $g^{(1)}(\alpha,p)$. The analogue of (4.1) becomes

[TABLE]

From (4.8) and (4.18), we find

[TABLE]

Hence the analogue of (4.9) becomes:

[TABLE]

The new arguments of $\psi_{1}$ are obtained from the old ones by a constant shift. Clearly, this does not affect the averaging argument in (4.10)–(4.16). The interval $[t_{1}^{(\text{min})},t_{1}^{(\text{max})}]$ needs to be enlarged slightly to account for the dependence of $\Gamma_{t}$ on $h$. The derivation works the same way, and we get the same result as in (4.17) even if $x_{h}$ replaces $x_{0}$ on the left side of (4.17). Thus, we proved the following result.

Theorem 2**.**

Pick a generic point $x_{0}\not\in S$. For any $\Theta_{0}\in S^{2}$, $h$ confined to a bounded set, and any $\chi\in C^{\infty}(S^{2})$ one has

[TABLE]

Combining (3.21), (3.23), (3.26), (3.27), (3.29), and (4.21), we finish the proof of Theorem 1.

5. Numerical experiments

We start by constructing an interpolation kernel $\varphi$ with the required properties. We use the method of [1]. In [1] the authors construct piecewise-polynomial kernels $\varphi$, which are characterized by a quadruple of numbers $\{N,W,R,L\}$. Here

1. (1)

$N$ is the maximal degree of polynomial pieces that make up $\varphi$, 2. (2)

$W$ is the support: $\text{supp}(\varphi)=[0,W]$, 3. (3)

$R$ is regularity, i.e. $\varphi$ is $R$ times continuously differentiable, and 4. (4)

$L$ is order, i.e. the approximation error behaves like $O(\epsilon^{L})$, where $\epsilon$ is step-size.

Additional requirements, e.g. that $\varphi$ be symmetric and interpolating can be imposed too.

Assumption A3 implies that $\varphi\in C^{2}(\mathbb{R})$. Since $\varphi$ is piecewise-polynomial, $\varphi^{\prime\prime\prime}$ is automatically piecewise continuous and bounded, so Assumptions A3 and A4 are satisfied if $R=2$. As is easy to see, Assumption A1 is satisfied if $L=3$. Assumption A2 means that $W<\infty$. Hence we should have $R=2$, $L=3$, and $W<\infty$. The construction in [1] involves a number of free parameters (degrees of freedom) $P$. In our case there are

[TABLE]

degrees of freedom. The requirement that $\varphi$ be interpolating: $\varphi((W/2)+n)=\delta_{n}$ imposes $W$ (respectively, $W-1$) additional conditions if $W$ is odd (respectively, even). Consequently, $(N-2)(W-3)\geq W$ (resp., $W-1$) if $W$ is odd (resp., even). Hence, $N\geq 4$ and $W\geq 4$. Choosing $N=4$, gives the minimal acceptable value $W=6$. Thus, $\varphi$ can be selected in the class $\{N=4,W=6,R=2,L=3\}$.

Imposing the condition that $\varphi$ be symmetric and using equations (20), (22), and (23) in [1] implies that $\varphi$ can be represented in the form

[TABLE]

where $a_{1,2,3}$ are some constants, and $B_{n}$ denotes the cardinal B-spline of degree $n$ supported on $[0,n+1]$. Finding $a_{1,2,3}$ by solving the system $\varphi((W/2)+n)=\delta_{n}$ and simplifying gives the final expression

[TABLE]

Using the properties of B-splines, we see that $\int_{\mathbb{R}}\varphi(t)dt=1$.

The test object consists of two balls with centers $c_{1}=(0,0,-5)$ and $c_{2}=(-5.52,0,-7.36)$. Each ball has radius $R=4$ and uniform density 1. The point on the boundary $x_{0}$, in a neighborhood of which we study resolution, is given by

[TABLE]

In agreement with our convention, $\Theta_{0}$ points into the interior of the first ball. The second ball is chosen in such a way that the plane through $x_{0}$ and normal to $\Theta_{0}$ is tangent to its boundary. This way we demonstrate that remote singularities do not contribute to edge response.

To simulate discrete data, Radon transform is computed at the points

[TABLE]

Thus, $\epsilon=\Delta p=0.04$. The predicted response is computed using (2.15), where $f_{0}=\rho_{0}=1$. The results are shown in Figure 1. The $x$-axis in the figure shows the $h$-values. The $y$-axis shows the reconstructed values at the points $x_{h}=x_{0}+\epsilon h\Theta_{0}$ (cf. (2.11)) as well as the predicted transition curve (edge response). The results show a good match.

Next we illustrate the situation where condition (2) does not hold and Theorem 1 does not hold either. The test object consists of a single ball of radius 1 with center $c=(0,0,0)$ in the first experiment, and center $c=(0,0,1)$ - in the second experiment. The point on the boundary $x_{0}$ is now given by

[TABLE]

As is easily seen, in both cases the components of the vector $\mathcal{H}_{0}^{\prime}$ (cf. (3.13)) are integers, so $x_{0}$ is not generic. The results are shown in Figure 2. The plot on the left is for the case $c=(0,0,0)$, and on the right - for the case $c=(0,0,1)$. The match between the predicted and actual transition behaviors is no longer accurate.

Appendix A Auxiliary results

Lemma 2**.**

Given an interval $[a,b]$, consider a real-valued function $f\in C^{2}([a,b])$ such that either (1) $f^{\prime\prime}\equiv 0$ on $[a,b]$ and $f^{\prime}$ equals to an irrational number, or (2) $f^{\prime\prime}\not=0$ on $[a,b]$, and

[TABLE]

for some $n\in\mathbb{Z}$. Then the points $\{f(\epsilon i)/\epsilon\}$, $\epsilon i\in[a,b]$, are u.d. as $\epsilon\to 0$.

Proof.

Case (1) is well-known (see [10]), so we only consider case (2). Using Weyl’s criterion, we wish to show that

[TABLE]

for all integers $M\not=0$, where $e(t):=\exp(2\pi it)$. Let $M\geq 1$ be fixed in what follows. The case $M\leq-1$ can be reduced to the case $M\geq 1$ by complex conjugation.

Without loss of generality, we may assume that $f^{\prime\prime}>0$ and $f^{\prime}(a)<f^{\prime}(b)$. Let $a_{m}$ be the points such that

[TABLE]

Suffice it to show that

[TABLE]

for each $m$, $0\leq m<M$.

Partition $[a_{m},a_{m+1}]$ into the intervals of length $\epsilon^{2/3}$:

[TABLE]

By assumption, $f^{\prime\prime}(t)$ is bounded away from zero on $[a,b]$, so each subinterval $\Delta_{k}$ contains $O\left(\epsilon^{-1/3}\right)$ points from the original set $\epsilon i$.

First, we show that

[TABLE]

Since

[TABLE]

it follows from (A.6)

[TABLE]

The big-$O$ terms in (A.7) and (A.8) are uniform in $k$, because $\max_{t\in[a,b]}f^{\prime\prime}(t)<\infty$ uniformly as $\epsilon\to 0$.

To prove (A.4), partition the sum and then use (A.6), (A.8) to obtain an estimate:

[TABLE]

In $J_{m}^{(l)}$, the sum is over $k$ such that

[TABLE]

In $J_{m}^{(r)}$, the sum is over $k$ such that

[TABLE]

Each of the sums is over $O(\epsilon^{-2/3})$ consequtive values of $k$. In particular, the number of terms in each sum is bounded from below by $c\epsilon^{-2/3}$ for all $\epsilon>0$ sufficiently small, where $c>0$ is some constant.

Consider $J_{m}^{(l)}$. We use the inequality

[TABLE]

where $\|r\|$ denotes the distance from $r$ to the nearest integer. By (A.1), (A.3), (A.5), and (A.10),

[TABLE]

for all $k$ used in $J_{m}^{(l)}$. Using (A.12) in (A.9) with

[TABLE]

gives

[TABLE]

This implies,

[TABLE]

The quantity $J_{m}^{(r)}$ can be estimated in a similar fashion (using (A.11) instead (A.10)), and the lemma is proven. ∎

Corollary 1**.**

Given an interval $[a,b]$, pick a real-valued function $f\in C^{2}([a,b])$. Suppose that

1. (1)

The set $\{x\in[a,b]:\,f^{\prime\prime}(x)=0\}$ is the union of a finite number of distinct points and non-intersecting intervals, and 2. (2)

If $f^{\prime\prime}(x)\equiv 0$ on an interval, then $f^{\prime}(x)$ is irrational on that interval.

Then the points $\{f(\epsilon i)/\epsilon\}$, $\epsilon i\in[a,b]$, are u.d. as $\epsilon\to 0$.

Proof.

Let $X$ be the set consisting of the isolated points where $f^{\prime\prime}(x)=0$ and of the endpoints of the intervals where $f^{\prime\prime}(x)\equiv 0$. By assumption, $X$ contains finitely many points. Let $X_{\delta}$ denote the $\delta$-neighborhood of $X$. Clearly, $[a,b]\setminus X_{\delta}$ is a finite union of non-overlapping intervals $[a_{i},b_{i}]$ such that (1) $|b-a|-\sum_{i}|b_{i}-a_{i}|$ can be made as small as we like by choosing $\delta>0$ small, and (2) each of the intervals $[a_{i},b_{i}]$ satisfies the assumptions of Lemma 2. Thus the points $\{f(\epsilon i)/\epsilon\}$, $\epsilon i\in\cup_{i}[a_{i},b_{i}]$, are u.d. as $\epsilon\to 0$. The fraction of the remaining points can be made arbitrarily small, which proves the desired assertion. ∎

Corollary 2**.**

Given an interval $[a,b]$, pick two functions $g,h$. Suppose the linear combination $f:=M_{1}g+M_{2}h$ satisfies the assumptions of Corollary 1 for any pair $M_{1},M_{2}\in\mathbb{Z}$, $M_{1}^{2}+M_{2}^{2}>0$. Then the set of points $(\{g(\epsilon i)/\epsilon\},\{h(\epsilon i)/\epsilon\})$, $\epsilon i\in[a,b]$, is u.d. in $[0,1]\times[0,1]$ as $\epsilon\to 0$.

Proof.

The assertion follows immediately from [10], Theorem 6.3, and Corollary 1. ∎

Lemma 3**.**

Let the set $(x_{\epsilon}(i),y_{\epsilon}(i))$, $i\in I_{\epsilon}$, be u.d. in the square $[0,1]\times[0,1]$ as $\epsilon\to 0$. Here $I_{\epsilon}$ is some set of indices that depends on $\epsilon$. Then the set $(\{x_{\epsilon}(i)+ay_{\epsilon}(i)\},y_{\epsilon}(i))$, $i\in I_{\epsilon}$, is also u.d. in $[0,1]\times[0,1]$ as $\epsilon\to 0$ for any $a\in\mathbb{R}$.

Proof.

Pick some $r_{0},y_{0},\Delta r,\Delta y>0$ so that $r_{0}+\Delta r,y_{0}+\Delta y\in[0,1)$. Consider the domain

[TABLE]

As is seen from elementary geometry (proof by picture), the area of $D$ is independent of $r_{0},y_{0}$ and equals $\Delta r\Delta y$. Since $(x_{\epsilon}(i),y_{\epsilon}(i))$ is u.d. in $[0,1]\times[0,1]$, the desired assertion follows. ∎