Exact and fast inversion of the approximate discrete Radon transform from partial data

TL;DR

This paper presents an exact and computationally efficient inversion method for a specific approximate discrete Radon transform, capable of reconstructing images from partial data with minimal computational cost.

Contribution

It introduces a novel exact inversion formula for the approximate discrete Radon transform that operates efficiently on partial data with O(N log N) complexity.

Findings

Exact inversion formula for the approximate discrete Radon transform.

Computational complexity of O(N log N) for 2D images.

Effective reconstruction from partial data.

Abstract

We give an exact inversion formula for the approximate discrete Radon transform introduced in [Brady, SIAM J. Comput., 27(1), 107--119] that is of cost $O(N \log N)$ for a square 2D image with $N$ pixels and requires only partial data.

Equations

$$A_{n,m-1}^{\ell }(i,j):=\{\begin{array}{ll}A(i,j+(\ell -1)2^{m-1})& \text{ if }i=1,\cdots ,2^n,j=1,\cdots ,2^{m-1},\\ 0& \text{ otherwise,}\end{array}$$

$$\{\begin{array}{rl}{D}_m(h,2t+1)& :=D_{m-1}^1(h,t)\cup D_{m-1}^2(h+t+1,s),\\ D_m(h,2t)& :=D_{m-1}^1(h,t)\cup D_{m-1}^2(h+t,s),\end{array}$$

$$D_m^1(h,s):=D_m(h,s),D_m^2(h,s):=\{(i,j+2^m):(i,j)\in D_m(h,s)\},$$

$$R_{n,m}^{\ell }(h,s):=\sum _{(i,j)\in D_m(h,s)}{A}_{n,m}^{\ell }(i,j),\text{ where }\ell =1,\cdots ,2^{n-m}.$$

$${\Delta }_{n,m-1}^{\ell }(h,s):=R_{n,m-1}^{\ell }(h+1,s)-R_{n,m-1}^{\ell }(h,s),$$

$${\Delta }_{n,m-1}^{2\ell -1}(h,s)$$

$${\Delta }_{n,m-1}^{2\ell }(h,s)$$

$$R_{n,m}^{\ell }(h+1,2s)$$

$$=\sum _{(i,j)\in D_{m-1}^1(h+1,s)}{A}_{n,m-1}^{2\ell -1}(i,j)+\sum _{(i,j)\in D_{m-1}^2(h+s+1,s)}{A}_{n,m-1}^{2\ell }(i,j)$$

$$=R_{n,m-1}^{2\ell -1}(h+1,s)+\sum _{(i,j)\in D_{m-1}(h,s)}{A}_{n,m-1}^{2\ell }(i+s+1,j),$$

$$R_{n,m}^{\ell }(h,2s+1)$$

$$=\sum _{(i,j)\in D_{m-1}^1(h,s)}{A}_{n,m-1}^{2\ell -1}(i,j)+\sum _{(i,j)\in D_{m-1}^2(h+s+1,s)}{A}_{n,m-1}^{2\ell }(i,j)$$

$$=R_{n,m-1}^{2\ell -1}(h,s)+\sum _{(i,j)\in D_{m-1}(h,s)}{A}_{n,m-1}^{2\ell }(i+s+1,j),$$

$$R_{n,m-1}^{\ell }(h,s)=\sum _{k=-\infty }^{h-1}{\Delta }_{n,m-1}^{\ell }(k,s)=\sum _{k=-s}^{h-1}{\Delta }_{n,m-1}^{\ell }(k,s),$$

$$\sum _{m=1}^n{2}^{n-m}\cdot [2\cdot 2^{m-1}(2^{n+1}+2^m+1)]=(2N+\sqrt{N})\log N+4N-2\sqrt{N},$$

Full text

1 Introduction

The Radon transform is a linear transform that integrates a real-valued function on a $d$-dimensional Euclidean space over all possible $(d-1)$-dimensional hyperplanes [3]. The transform is a widely used model in tomography, and many inverse problems involving the transform or its variant have been carefully studied [4].

We will refer to the Radon transform as the continuous Radon transform, to distinguish it from the approximate discrete Radon transform (ADRT), which is the focus of this paper; a discrete version of the transform that replaces the smooth hyperplanes by broken pixelated lines, called digital lines [1, 2]. ADRT is computed by summing the pixel values that lie in these lines. The digital lines are computed recursively, yielding a fast transform of computational cost $\mathcal{O}(N\log N)$ for a square 2D image with $N$ number of pixels.

We will show that the exact inverse can be computed with the same complexity, in contrast to other discrete versions of the continuous transform. This is somewhat surprising: despite the known inversion formula for the continuous case, computing the inverse requires more effort than the forward transform, and the digital lines in ADRT converges to straight lines upon repeated refinement. Moreover, the inverse can be computed using only partial data, that is, one quadrant of the ADRT that correspond to angles in $[0,\pi/4]$ for the continuous transform.

It was observed that the inverse of ADRT can be computed to numerical precision by employing a multigrid method [5], and this was further validated using the conjugate gradient method [6]. However, the simple formula in this work is new, to the best of our knowledge.

2 Main Result

Our main assertion is that the ADRT is exactly invertible from partial data, with the same computational complexity as the forward transform.

Theorem 2.1**.**

A 2D square image with $N$ pixels can be computed exactly from its single-quadrant ADRT in $\mathcal{O}(N\log N)$ operations.

Let us be given a 2D square image $A$ containing $N$ pixels taking on real values, with the dimensions $2^{n}\times 2^{n}$. We will first define the ADRT for a rectangular sub-images of size $2^{n}\times 2^{m}$ for $m<n$, called sections.

Definition 2.2** (Section of an image).**

Let the $\ell$-th section of a $2^{n}\times 2^{n}$ image $A$ be defined on $\mathbb{Z}^{2}$ given by

[TABLE]

where $\ell=1,2,\cdots,2^{n-m}$.

Next, we define the broken lines used in ADRT.

Definition 2.3** (Digital line).**

A digital line $D_{m}(h,s)$, for $h\in\mathbb{Z}$, $s=1,\cdots,2^{m}$ is a subset of $\mathbb{Z}^{2}$ that is defined recursively. Letting $s=2t$ or $s=2t+1$,

[TABLE]

for $t=1,\cdots,2^{m-1}$ and where

[TABLE]

and the relation is initialized by $D_{0}(h,s):=\{(h,1):h\in\mathbb{Z}\}$.

Then ADRT is a sum of pixels that lie on the digital lines.

Definition 2.4** (Single-quadrant ADRT).**

We will denote the single-quadrant approximate discrete Radon transform (ADRT) applied to the $\ell$-th section of the image $A$ by $R^{\ell}_{n,m}$, defined to be the sum

[TABLE]

In particular, when $m=n$ we call the sums simply the single-quadrant ADRT of $A$, and we denote it by $R_{m}:=R^{1}_{n,n}$.

A full ADRT is obtained when one computes a single-quadrant ADRT four times upon flipping and rotating the image [1].

We state some basic properties of digital lines and the ADRT, omitting the proof.

Lemma 2.5** (Properties).**

1. (i)

$\{D_{m}(h,s):h\in\mathbb{Z}\}$ form a partition of the set $\{(i,j):i\in\mathbb{Z},j=1,\cdots,2^{m}\}$. 2. (ii)

$D_{m}^{1}(h,s)\cap D_{m}^{2}(k,t)=\emptyset$, for all $h,k\in\mathbb{Z}$ and $s,t=1,\cdots,2^{m}$. 3. (iii)

$D_{m}(h,s)\subset\{(i,j)\in\mathbb{Z}^{2}:h\leq i\leq h+s,1\leq j\leq 2^{m}\}$. 4. (iv)

$R^{\ell}_{n,0}(h,1)=A(h,\ell)$ for $h\in\mathbb{Z},\ell=1,\cdots,2^{m}$. 5. (v)

$R^{\ell}_{n,m}(h,s)=0$ when $h<-s+1$ or $h>2^{n}$, $s=1,\cdots,2^{m}$.

The key lemma below states that the recursive definition in ADRT is reversible: the ADRT of a section of an image can be computed from that of the whole.

Note that by Lemma 2.5, $R^{\ell}_{n,m}$ can have at most $2^{n}\cdot 2^{m}+2^{m-1}(2^{m}+1)$ non-zero entries.

Lemma 2.6**.**

Given $R^{\ell}_{n,m}$, $\{R^{2\ell-1}_{n,m-1},R^{2\ell}_{n,m-1}\}$ can be computed in $\mathcal{O}(M)$ operations, where $M:=2^{m-1}(2^{n+1}+2^{m}+1)$ is the size of $R^{\ell}_{n,m}$.

Proof 2.7**.**

Let us define the differences

[TABLE]

then these differences can be computed by the relations

[TABLE]

They follow from direct computation. For example,

[TABLE]

and also

[TABLE]

so we have Eq. 2.6 upon subtraction. Since we now have the differences,

[TABLE]

since part of the sum vanishes due to Lemma 2.5.

For each fixed $s$, this sum allows one to compute all values of $h$ in one sweep with one addition per one $h$: so the number of additions or subtractions for computing $R^{\ell}_{n,m-1}$ for all $(h,s)$ is $2M=2\cdot 2^{m-1}(2^{n+1}+2^{m}+1)=\mathcal{O}(N)$.

The main theorem follows.

Proof 2.8** (Theorem 2.1).**

By Lemma 2.6, we compute $\{R^{\ell}_{n,m-1}:\ell=1,\cdots,2^{n-m+1}\}$ from $\{R^{\ell}_{n,m}:\ell=1,\cdots,2^{n-m}\}$. Repeating this procedure for $m=1,\cdots,n$, one computes $R^{\ell}_{n,0}$ and by doing so one obtains the original image $A$.

The cost is then

[TABLE]

that is, $\mathcal{O}(N\log N)$.

3 Conclusion

We showed that the inverse of ADRT can be computed from $\mathcal{O}(N\log N)$ operations, using data from only one quadrant of the ADRT. This could be useful for applications in computational applications of the ADRT, e.g. for use in dimensional splitting or in approximating the continuous transform. It also has implications for implementation, especially when decomposing the domain for parallelization. This result indicates that the relation between the ADRT and the continuous transform is more involved than it may first appear.