## EXACT AND FAST INVERSION OF THE APPROXIMATE DISCRETE RADON TRANSFORM

DONSUB RIM ∗

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.

- 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 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 , π/ 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 1. A 2D square image with N pixels can be computed exactly from its single-quadrant ADRT in O ( N log N ) operations.

Definition 2 (Section of an image). Let the /lscript -th section of a 2 n × 2 n image A be defined on Z 2 given by

Let us be given a 2D square image A containing N pixels taking on real values, with the dimensions 2 n × 2 n . We will first define the ADRT for a rectangular subimages of size 2 n × 2 m for m&lt;n , called sections.

$$( 2. 1 ) \quad A _ { n, m - 1 } ^ { \ell } ( i, j ) \coloneqq \begin{cases} A ( i, j + ( \ell - 1 ) 2 ^ { m - 1 } ) & i f \, i = 1, \cdots, 2 ^ { n }, j = 1, \cdots, 2 ^ { m - 1 }, \\ 0 & o t h e r w i s e, \end{cases}$$

where /lscript = 1 2 , , · · · , 2 n -m .

∗ Courant Institute of Mathematical Sciences, New York University, 251 Mercer St., New York, NY 10012 ( dr1653@nyu.edu ).

Fig. 1: Examples of digital lines (d-lines) D m ( h, s ) where m = 3 (left) and a diagram illustrating the recursion definition (2.2) (middle, right). In both figures the case when s is even is depicted in blue, odd in black, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00887v3-figure-55734abd-2bde-4636-86c6-e7d8a509fe1e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00887v3-figure-268f8eda-6275-4f76-a824-ed4f7546b591.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00887v3-figure-a3ea8cea-348a-422f-82ff-b4f7af6a8a39.jpg)

Next, we define the broken lines used in ADRT.

Definition 3 (Digital line). A digital line D m ( h, s ) , for h ∈ Z , s = 1 , · · · , 2 m is a subset of Z 2 that is defined recursively. Letting s = 2 t or s = 2 t +1 ,

$$( 2. 2 ) \quad \begin{cases} D _ { m } ( h, 2 t + 1 ) \coloneqq D _ { m - 1 } ^ { 1 } ( h, t ) \cup D _ { m - 1 } ^ { 2 } ( h + t + 1, s ), \\ D _ { m } ( h, 2 t ) \coloneqq D _ { m - 1 } ^ { 1 } ( h, t ) \cup D _ { m - 1 } ^ { 2 } ( h + t, s ), \\. \quad. \quad \dots \cdot \quad. \end{cases}$$

for t = 1 , · · · , 2 m -1 and where

$$( 2. 3 ) \quad \ D _ { m } ^ { 1 } ( h, s ) \coloneqq D _ { m } ( h, s ), \ \ D _ { m } ^ { 2 } ( h, s ) \coloneqq \{ ( i, j + 2 ^ { m } ) \colon ( i, j ) \in D _ { m } ( h, s ) \},$$

and the relation is initialized by D 0 ( h, s ) := { ( h, 1) : h ∈ Z } .

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

Definition 4 (Single-quadrant ADRT). We will denote the single-quadrant approximate discrete Radon transform (ADRT) applied to the /lscript -th section of the image A by R /lscript n,m , defined to be the sum

$$( 2. 4 ) \quad R ^ { \ell } _ { n, m } ( h, s ) \coloneqq \sum _ { ( i, j ) \in D _ { m } ( h, s ) } A ^ { \ell } _ { n, m } ( i, j ), \quad w h e r e \ \ell = 1, \cdots, 2 ^ { n - m }. \\ \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \,.$$

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 5 (Properties).

- (i) { D m ( h, s ) : h ∈ Z } form a partition of the set { ( i, j ) : i ∈ Z , j = 1 , · · · , 2 m } .

(iii) D m ( h, s ) ⊂ { ( i, j ) ∈ Z 2 : h ≤ i ≤ h + s, 1 ≤ j ≤ 2 m } .

- (ii) D 1 m ( h, s ) ∩ D 2 m ( k, t ) = ∅ , for all h, k ∈ Z and s, t = 1 , · · · , 2 m .
- (iv) R /lscript n, 0 ( h, 1) = A h, /lscript ( ) for h ∈ Z , /lscript = 1 , · · · , 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.

- (v) R /lscript n,m ( h, s ) = 0 when h &lt; -s +1 or h &gt; 2 n , s = 1 , · · · , 2 m .

Note that by Lemma 5, R /lscript n,m can have at most 2 n · 2 m +2 m -1 (2 m +1) non-zero entries.

Lemma 6. Given R /lscript n,m , { R 2 /lscript -1 n,m -1 , R 2 /lscript n,m -1 } can be computed in O ( M ) operations, where M := 2 m -1 (2 n +1 +2 m +1) is the size of R /lscript n,m .

Proof. Let us define the differences

$$\Delta _ { n, m - 1 } ^ { \ell } ( h, s ) \coloneqq R _ { n, m - 1 } ^ { \ell } ( h + 1, s ) - R _ { n, m - 1 } ^ { \ell } ( h, s ), \\ \therefore \dots \dots \dots \dots \dots \dots$$

then these differences can be computed by the relations

$$& ( 2. 6 ) \quad \quad \Delta _ { n, m - 1 } ^ { 2 \ell - 1 } ( h, s ) = R _ { n, m } ^ { \ell } ( h + 1, 2 s ) - R _ { n, m } ^ { \ell } ( h, 2 s + 1 ), \\ & \omega \ \pi \quad \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$\Delta _ { n, m - 1 } ^ { 2 \ell } ( h, s ) = R _ { n, m } ^ { \ell } ( h - s, 2 s + 1 ) - R _ { n, m } ^ { \ell } ( h - s, 2 s ). \\ \dots \quad \hat { \ } \cdots \quad \hat { \ } \cdots \quad \dots \quad \dots \quad \cdots \quad \cdots$$

They follow from direct computation. For example,

$$T \text{hey follow from direct computation. For example,} \\ R _ { n, m } ^ { \ell } ( h + 1, 2 s ) & = \sum _ { ( i, j ) \in D _ { m } ( h + 1, s ) } A _ { n, m } ^ { \ell } ( i, j ) \\ & = \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 ), \\ \text{and also}$$

and also

$$\text{and also} \\ R _ { n, m } ^ { \ell } ( h, 2 s + 1 ) & = \sum _ { ( i, j ) \in D _ { m } ( h, s ) } A _ { n, m } ^ { \ell } ( i, j ) \\ & = \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 ), \\ \text{so we have (2.6)} \, \text{upon subtraction. Since we now have the differences,}$$

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

$$& ( 2. 8 ) & R ^ { \ell } _ { n, m - 1 } ( h, s ) = \sum _ { k = - \infty } ^ { h - 1 } \Delta ^ { \ell } _ { n, m - 1 } ( k, s ) = \sum _ { k = - s } ^ { h - 1 } \Delta ^ { \ell } _ { n, m - 1 } ( k, s ), \\ & \text{since part of the sum vanishes due to Lommas 5}$$

since part of the sum vanishes due to Lemma 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 /lscript n,m -1 for all ( h, s ) is 2 M = 2 · 2 m -1 (2 n +1 +2 m +1) = O ( N ).

The main theorem follows.

Proof (Theorem 1 ). By Lemma 6, we compute { R /lscript n,m -1 : /lscript = 1 , · · · , 2 n -m +1 } from { R /lscript n,m : /lscript = 1 , · · · , 2 n -m } . Repeating this procedure for m = 1 , · · · , n , one computes R /lscript n, 0 and by doing so one obtains the original image A .

The cost is then

$$& ( 2. 9 ) \quad \sum _ { m = 1 } ^ { n } 2 ^ { n - m } \cdot [ 2 \cdot 2 ^ { m - 1 } ( 2 ^ { n + 1 } + 2 ^ { m } + 1 ) ] = ( 2 N + \sqrt { N } ) \log N + 4 N - 2 \sqrt { N }, \\ & \text{that is } \mathcal { O } ( N \log N )$$

that is, O ( N log N ).

- 3. Conclusion. We showed that the inverse of ADRT can be computed from 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.

## REFERENCES

- [1] M. L. Brady , A fast discrete approximation algorithm for the Radon transform , SIAM Journal on Computing, 27 (1998), pp. 107-119, doi:10.1137/S0097539793256673.
- [2] W. G¨ otz and H. Druckm¨ uller , A fast digital radon transform: an efficient means for evaluating the hough transform , Pattern Recognition, 29 (1996), pp. 711 -718, doi:10.1016/0031-3203(96)00015-5.
- [3] S. Helgason , Integral Geometry and Radon Transforms , Springer New York, New York, NY, 2011, doi:10.1007/978-1-4419-6055-9.
- [4] F. Natterer , The Mathematics of Computerized Tomography , Society for Industrial and Applied Mathematics, 2001, doi:10.1137/1.9780898719284.
- [5] W. H. Press , Discrete Radon transform has an exact, fast inverse and zes to operations other than sums along lines , Proceedings of the National Academy of Sciences, 103 (2006), pp. 19249-19254, doi:10.1073/pnas.0609228103.
- [6] D. Rim , Dimensional splitting of hyperbolic partial differential equations using the Radon transform , SIAM Journal on Scientific Computing, 40 (2018), pp. A4184-A4207, doi:10.1137/17M1135633.
Donsub Rim
Exact and fast inversion of the approximate discrete Radon transform from partial data
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.