# Selecting Rays for Simulated Radon Transform

I have the task of determining approximations of a 2D function $f: (x,y)\in \mathbb{R}^2\mapsto\mathbb{R}$ from integrals along lines, i.e. from its Radon transform $R(\phi,\tau)[f(x,y)]$ and, because of the greater flexibility, I want to solve it by interpreting it as a least squares problem, which I define as follows:

• the variables correspond to the cells of a rectangular grid, that covers the relevant range of $f$'s arguments

• each equation corresponds to a different line(-segment), where each variable is weighted with the length of the intersection of that line(-segment) with the rectangular cell that corresponds to the variable and the right hand side corresponds to the integral of $f(x,y)$ along that line(-segment).
If the line misses a cell, then the length of the intersection is defined to be $0$.

Now, the problem and freedom I have, is to determine a set of line(-segments), from which the cell-weights that approximate the function values can be determined.

It is clear, that the number of equations must at least equal the number of cells but, as there are no other restrictions, there are various options for choosing the line(-segments) that yield the equations, with different implications on the numerical properties of the resulting least squares problem.

Question:
what can be recommended for generating potentially infinite sequences of lines for setting up the least squares problems, whose solutions will yield increasingly better approximations of $f(x,y)$,
i.e. how to select the angles and anchor points, that define the line(segment)s.

For the directions of the lines, the following options came to my mind:

• distribute them evenly over the range $[-\pi,\pi)$
• generate them via a low-discrepancy sequence over $[-\pi,\pi)$ e.g., van der Corput sequence.

• generate them via a sequence of Pythagorean triples $(a,b,c)\in\mathbb{N}^3: a^2+b^2=c^2$, setting $\phi := arctan(\frac{a}{b})$

For the anchor points the following strategies seem reasonable:

• use the center of cell

• use a corner of a cell

if the line integrals happen to be integers, then combining the "Pythagorean directions" with cell-corners as anchor points guarantees, that the tile weights are also integers, which may increase numeric precision.

• If you want to get increasingly good approximations of $f$, do you want to increase the resolution where $f$ seems to have interesting features or everywhere? I would add lines that go through the points where you want more information, trying to cover the normal bundle of the wavefront set of $f$. This could mean adding lines tangent to the level sets of $f$ where its gradient seems large. Do you want your approach to be somehow dynamic or independent of the actual data? Oct 13, 2014 at 9:14
• @JoonasIlmavirta for the first iteration, I would like to increase resolution everywhere. The dynamic improvement would be the further steps; that would consist either sticking to the fixed grid, but shooting more ray through interesting places or, switching to a quad-tree decomposition or, a combination of both. I want the approach to independent of the actual data, but I would also appreciate recommendations for special situations. Oct 13, 2014 at 9:45
• Then I think it would be best to distribute angles evenly and use the same spacing between parallel lines independent of the direction. If the lines are not evenly distributed, I think the reconstruction could be distorted. Uneven distribution can also be compensated by weighting different lines differently in the least squares method. Further details seem to call for experience with numerics rather than the continuous (theoretical) Radon transform, so I can only suggest experimenting and trying to cover points and directions evenly. Oct 13, 2014 at 10:47
• @JoonasIlmavirta my concern is primarily a numeric one and, evenly spaced sampling is the simplest strategy, but not the best as can be learned from approximating functions by polynomials. Oct 14, 2014 at 4:53

If your question is motivated by practical applications of the Radon transform such as computerized tomography in medical imaging or non-destructive testing in industry, then I guess that your question will not be regarded as "well posed" in these communities. For one, the statement "It is clear, that the number of equations must at least equal the number of cells." is not totally true in the applied area. There you have a lot of prior knowledge on the function $f$ which should be exploited, such as non-negativity, a good estimate of the integral (or 1-norm) of the function, and probably further assumptions such as smoothness, piecewise smoothness, regularity of the jump set and/or sparsity assumptions. All these prior knowledge makes it possible to obtain good or reasonable reconstruction of the objects under consideration even from underdetermined measurements.
Another problem is that of "mathematical ill-posedness". The inversion of the Radon transform is ill-posed in the sense that the forward operator $R$ (i.e. the Radon transform) is smoothing (roughly increases the Sobolev smoothness by 1/2). Hence, the operator $R$ does not have a continuous inverse. In consequence, there will be problems with "approximating the continuous Radon transform finer and finer" since you are approximating an infinite matrix with infinite condition number. The condition number of you approximating problems will grow unbounded. It is natural in this setting to use "regularization by discretization", that is, to balance accuracy of the approximation of the operator with computational stability/noise amplification.