Inverse Problems Course Notes — The X-Ray Transform

We motivated the study of the Radon transform with a tomographic problem: given the change in intensity of X-rays along all lines through a region, can we reconstruct the attenuation (think of this as density) in the region?

For our purposes in this problem, light travels along straight lines and in two dimensions, those lines are hyperplanes.  So inverting the Radon transform — which sums functions over hyperplanes — solves our problem in 2-D.  But in higher dimensions this is not what we need.  We need to integrate along lines, not hyperplanes.

So we introduce the X-ray transform.

Definition Given a function f: \mathbb{R}^n \rightarrow \mathbb{C}, the X-ray transform of f is defined as

Xf(x, \theta) = \int_{\mathbb{R}} f(x + t\theta) dt

where x\in \mathbb{R}^n, \theta \in S^{n-1}.

The first thing to notice is that there is some redundancy here. For all x, \theta, s we will have

Xf(x + s\theta, \theta) = \int_{\mathbb{R}} f(x + s\theta + t\theta) dt = \int_{\mathbb{R}} f(x + \tilde{t}\theta) d\tilde{t} = Xf(x, \theta)

So the natural domain of Xf is actually smaller, we can pick just one representative from the line \{x + s\theta | s \in \mathbb{R} \} without losing any information.   There is a distinguished point on this line: the point E_{\theta}(x) \doteq x - \langle x, \theta\rangle\theta is the unique point on the line that is perpendicular to \theta, \langle E_{\theta}(x), \theta \rangle = 0.

Since Xf(x, \theta) = Xf(E_{\theta}(x), \theta) for all x \in \mathbb{R}^n, \theta \in S^{n-1}, it is natural to think of Xf as a function \theta \in S^{n-1} and x \in \theta^{\perp} = \{ x | \langle x, \theta \rangle = 0\}, the hyperplane normal to \theta.  Formally,  the domain of Xf is

T = \{(x, \theta) \in \mathbb{R}^n\times S^{n-1} | x \in \theta^{\perp} \}

Where \theta^{\perp} is the hyperplane normal to \theta.  This is a manifold of dimension 2n - 2.

Now \theta^{\perp} \equiv \mathbb{R}^{n-1}, and we will use standard Lebesgue measure on this space, denoting it by dH_{\theta^{\perp}}.  As with the hyperplanes in the Radon transform, we will keep using the fact that

dH_{\theta^{\perp}}dt = dx

is standard Lebesgue measure on $mathbb{R}^{n}$.

Function Spaces on T

To proceed we need to define the standard function spaces on this domain.  There is nothing surprising here. Let’s start with C^{\infty}

C^{\infty}(T) = \{ f: T \rightarrow \mathbb{C} | \partial_{\theta^{\perp}}^\alpha \partial_{\theta}^\beta f \in C^0(T) \forall \alpha, \beta\}

Here \partial_{\theta^{\perp}} is a tangential derivative to \theta^{\perp}. For example, if

\theta = (1, 0, \dots, 0)

then take

\partial_{\theta^{\perp}} = \partial_{x_j}, \quad\quad j = 2, \dots, n

For any other hyperplane, we can rotate it into this position and use these operators to define \partial_{\theta^{\perp}}.

Now we will defineC_0^{\infty}.  The sphere is compact, so we only need to worry about compact support in the \theta^{\perp} directions.

C_0^{\infty}(T) = \{ f \in C^\infty(T) | \exists R s.t.  f(x,\theta) = 0  \forall |x| \geq R \}

(Sorry — having wordpress-latex trouble. Will try to get this fixed.) The Schwartz space is similar:

\mathcal{S}(T) = \{ f \in C^{\infty}(T) | \sup_{x \in \theta^{\perp}}|x^{\beta}\partial^{\alpha}_{\theta^{\perp}}\partial^{\gamma}_{\theta} f | \leq C_{\alpha, \beta, \gamma} < \infty \}

As with the spaces defined for the radon transform, we could choose to ignore the derivatives in the angular variables without impacting the subsequent theory, but we will not do that.

In the same way, we can define \mathcal{E}^\prime, \mathcal{D}^\prime, \mathcal{S}^\prime, etc.  Topologies on these spaces are also defined in the natural way.

The Transpose of X

It is easy to see that X is injective: the Radon transform R is injective and R can be computed from X by integrating the value of Xf over all the lines that make up the hyperplanes in the definition of Rf.  But we will do better and produce an explicit inversion formula.

First we need to find the formal transpose of X, i.e. the operator X^t such that \langle Xf, g\rangle = \langle f, X^tg\rangle.

Let f \in C_0^{\infty}(\mathbb{R}^n), g \in C^{\infty}(T)

\langle Xf, g\rangle = \int_{S^{n-1}}\int_{\theta^{\perp}} Xf(x,\theta)g(x,\theta) dH_{\theta^{\perp}} d\theta

= \int_{S^{n-1}}\int_{\theta^\perp}\int_{\mathbb{R}} f(x + t\theta) g(x, \theta) dt dH d\theta

As x varies over \theta^\perp and t varies over \mathbb{R} in this integral, the variable y = x + t\theta varies over all of \mathbb{R}^n and

dy = dH_{\theta^\perp}dt

is standard Lebesgue measure. If we perform this change of variables, we must remember that g(\cdot, \theta) is only defined on \theta^\perp, but when y = x + t\theta, we have x = y - t\theta = y - \langle y, \theta\rangle\theta = E_\theta(y). So

\langle Xf, g\rangle = \int_{\mathbb{R}^n}f(y)\int_{S^{n-1}}g(E_\theta(y), \theta) d\theta dy

So we will define

X^tg(x) = \int_{S^{n-1}}g(E_\theta(x), \theta) d\theta

Proposition The maps

X: C^\infty_0(\mathbb{R}^n) \rightarrow C_0^\infty(T)

X^t: C^\infty(T) \rightarrow C^\infty(\mathbb{R}^n)

are linear and continuous.

Proof Exercise.

The Normal Operator

We will use a different approach to develop the inversion formula for the X-ray transform.  We could parallel our development for the Radon transform, and we could rederive the Radon inversion formula using this approach.  It is good to see more than one way to look at a problem.

In linear algebra, it is often convenient to work with an operator X^tX.  It is a normal operator and is usually “much better behaved” than X.

X^tXf(x) = \int_{S^{n-1}}Xf(E_\theta(x), \theta) d\theta

= \int_{S^{n-1}}\int_{\mathbb{R}}f(x + t\theta)dt d\theta

= 2\int_{S^{n-1}}\int_0^\infty f(x + t\theta) dt d\theta

Because in the integral defining X^tX, we  catch the integral over t \in \mathbb{R}^{-} twice: once for \theta, and once for -\theta.  This can be reinterpreted as an integral in polar coordinates. Let y = t\theta and t = |y|.  Then

X^tXf(x) = c \int_{\mathbb{R}^n} \frac{f(x-y)}{|y|^{n-1}} dy

= f * \frac{1}{|x|^{n-1}}

up to constants.  This is a singular integral, but not badly behaved.  We do not need to regularize it.

To invert this convolution, we use the Fourier transform.

\mathfrak{F}(X^tXf)(\xi) = c\hat{f}(\xi)\widehat{\frac{1}{|x|^{n-1}}}(\xi)

for f \in C_0^\infty(\mathbb{R}^n).

Claim \widehat{\frac{1}{|x|^{n-1}}}(\xi) = c_n|\xi|^{-1}

We will prove this claim (in fact, we’ll prove a more general result) in a supporting post.  This result immediately lets us invert the normal operator. If

AX^tXf(x) = f(x)

then we must have

\widehat{Af}(\xi) = |\xi|\hat{f}(\xi)

In other words,

Af = (-\bigtriangleup)^{\frac{1}{2}}

is a power of the negative Laplacian. Now we can state the main result

Theorem [The X-Ray Inversion Formula] For all f \in C_0^{\infty}(\mathbb{R}^n)

c(-\bigtriangleup)^{\frac{1}{2}}X^tXf = f

So we gain half a derivative by applying X. Notice that this formula does not depend on the dimension, because the dimension of the manifolds we are integrating over to define X do not change.  Also note that the inversion operator is never local.

We will see that for n > 2 this has too many lines, and the formula is not very useful.

Parallel Development for the Radon Transform

We could develop the inversion formula for the Radon transform using its normal operator in a completely analogous way.

R^tRf(x) = c_n\int\frac{f(y)}{|x-y|} dy = c_n f * \frac{1}{|x|}

So the same manipulations show that

Theorem [Radon Inversion Formula II] For f \in C_0^{\infty}(\mathbb{R}^n)

(-\bigtriangleup)^{\frac{n-1}{2}}R^tRf = f

If n is even, we get a square root and the operator is not local.  If n is odd, the inversion operator is local.

At first glance, this looks different, but it is really the same as what we had before.  The proof is left as an exercise.  It relies on the observation that the Radon transform intertwines (functions of the) Laplacian on \mathbb{R}^n with those on \mathbb{R}:

R(-\bigtriangleup)^\alpha = (\frac{\partial^2}{\partial s^2})^\alpha R

A Quick Look Ahead

In the coming lectures, we will

  • Generalize these results to distributions
  • Rewrite the inversion formula to make it look more like the formula for R.
  • Make stability estimates
  • Characterize the range of X

The range characterization is perhaps the most interesting part.  Xf depends on 2n - 2 variables and is overdetermined.  We can extend it to \mathbb{R}^n \times \mathbb{R}^n as follows

Xf(x, \xi) = \int f(x + t\xi) dt

and immediately see that we have a compatibility condition:

\frac{\partial^2 Xf}{\partial x_i\partial \xi_j} = \frac{\partial^2 Xf}{\partial x_j\partial \xi_i}

We will see the Fritz John’s theorem that states that this system of PDEs gives necessary and sufficient conditions for a function to be in the range of X.


2 thoughts on “Inverse Problems Course Notes — The X-Ray Transform

  1. In the section on “Parallel Development for the Radon Transform”,
    in order to get the representation that you have stated for R^{t}R, I believe one really needs to use the Fourier slice theorem (unless you have a quicker proof). So in a sense the derivation does not quite follow the one derived for the X-ray transform. What do you think?

  2. I think there is a direct proof, but I cannot produce it right now. I should try! Also, I have heard that this result is proven in Natterer’s book but need to make a trip to the library to see how it is done.

    I will grant this: I don’t think computing R^tR is as easy as computing X^tX and if I were asked to do it on the spot, my first reaction would be to use a Fourier Transform. So the parallel seems less than perfect at the lower level steps.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s