15

The Radon transform maps a function $f$ on points in $\mathbb R^2$ to a function $Rf$ on lines in $\mathbb R^2$, given by $$Rf(\xi) = \int_{\mathbf x\in\xi} f(\mathbf x)\,\|\mathrm d\mathbf x\|.$$ Mathematica can compute the Radon transform via the function Radon, and its inverse via InverseRadon.

The dual Radon transform maps a function $g$ on lines in $\mathbb R^2$ to a function $R^*g$ on points in $\mathbb R^2$, given by $$R^*g(\mathbf x) = \int_{\xi\ni\mathbf x} g(\xi)\,\mathrm d\mu(\xi),$$ where $\mathrm d\mu$ is the unique rotationally invariant probability measure on the set of lines through $\mathbf x$.

How can the dual Radon transform and its inverse be computed in Mathematica?

1 Answers1

15

The dual Radon transform can be computed using InverseRadon[#, Method -> None].

The inverse of the dual Radon transform is just the Radon transform filtered by the pseudodifferential operator $(-L)^{1/2}$, where $L$ is the "radial" second derivative. In the discrete Fourier domain, this amounts to multiplying by Min[#, 1 - #] & where # is the vertical coordinate.

However, since this is usually too dark to see, I will choose a different scaling and compensate for it in dualRadon.

dualRadon[g_, dim_] := Module[{w, h},
  {w, h} = ImageDimensions[g];
  InverseRadon[ImageMultiply[g, 1/h], dim, Method -> None]]

inverseDualRadon[f_] := 
 Module[{r, w, h, lhalf}, 
  r = Radon[f]; 
  {w, h} = ImageDimensions[r]; 
  lhalf = Table[Min[i, h - i], {i, 0, h - 1}, {j, 0, w - 1}]; 
  Image@Chop@InverseFourier[lhalf Fourier@ImageData@r]]

For example:

image = Import["https://i.stack.imgur.com/68Zv0.jpg"]

enter image description here

idr = inverseDualRadon[image]

enter image description here

(There are also some negative values in that image which you can't see.)

recovery = dualRadon[idr, ImageDimensions[image]]

enter image description here

It's off by some scalar multiplier, but then so is InverseRadon @* Radon, so I'm not going to worry about it.