# Computing the dual Radon transform and its inverse

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?

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]]

Module[{r, w, h, lhalf},
{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["http://i.stack.imgur.com/68Zv0.jpg"]


idr = inverseDualRadon[image]


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

recovery = dualRadon[idr, ImageDimensions[image]]


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