I have to calculate a 2-dimensional radially symmetric distribution from a single projection.

I know that InverseRadon should actually do the job, but I get the wrong scaling and also the wrong spatial frequencies. I know that it is possible to define a filtering function for InverseRadon by choosing a method, e.g. by Method -> (# &). Nevertheless, I tried quite many filter functions and never get the proper scaling or proper spatial frequencies.

In the following I have tried to find out what's wrong using one somewhat complex but well-defined radially symmetric function from which I generate the projection and then apply the inverse Radon on it:

komplFkt =
ArrayResample[
Table[5 Cos[2π Sqrt[x^2+y^2]] - (x^2 + y^2) + 6,{x,-3,3,0.1},{y,-3,3,0.1}],
{128, 128}];

projectionAlongYDir = Total/@komplFkt;

GraphicsRow[
{ListPlot3D[komplFkt, PlotRange -> All],
Automatic,
ColorFunction -> ArrayPlot[komplFkt, PlotLegends -> "TemperatureMap"],
ListPlot[projectionAlongYDir, PlotLabel -> "Projection along Y-direction"]}]


which gives

Left and middle: radially symmetric example function/data. Right: projected data along Y-direction

Even though I have only a single projection, I know that for a radially symmetric function I should get always the same projection, even if I look at it from a different angle. That's why I construct now a "virtual" Radon transform of my original data, by simply repetitively placing my projected data into the columns of an image:

makeVirtaulRDNfromSingleProj[projectedDataList_] :=
Image @ Transpose @ ConstantArray[projectedDataList,Length[projectedDataList]];



this gives then (watch out for the // ImageAdjust):

Now I perform the backprojection, i.e. the InverseRadon[] on the previously generated 'virtual' transform (that I get from my measurements):

result =
InverseRadon[makeVirtaulRDNfromSingleProj[#], {Length @ #, Length @ #}] & @ projectionAlongYDir;

ArrayPlot[ImageData @ result,
PlotLegends -> Automatic, ColorFunction -> "TemperatureMap"]


Which gives:

Recovered data from single projection

As we can see, the scale is completely different than the function I started with. Moreover, when we look at the cross-sections through the middle (i.e. data's origin) we also see different spatial frequencies by calling the following functions:

GraphicsRow[
{ListPlot[{komplFkt[[64]], (ImageData @ result)[[64]]},
Joined -> True,
PlotLabel -> "unscaled cross-sections\nthrough center",
PlotLegends -> {"original function", "after InverseRadon[]"}],
ListPlot[Rescale /@ {komplFkt[[64]], (ImageData @ result)[[64]]},
Joined -> True,
PlotLabel -> "rescaled cross-sections\nthrough center",
PlotLegends -> {"original function", "after InverseRadon[]"}]}]


Tested in Mathematica v11.3 and v12.0.

Just for comparison, I tried to do the same in MATLAB. Neglecting some small deviations in the middle and some larger deviations at the boundaries I get there almost the same cross section I would expect:

Ok, the obvious answer is then probably: "do your calcs in MATLAB then". But why does Mathematica's InverseRadon behave so completely differently from iradon? Is it a bug?

I know that the default filter setting in Mathematica is (1 + Cos[# Pi])/2&. Hann filter (default) and in MATLAB is Ram-Lak (which should correspond in Mathematica to setting of Method -> (#&)).

But no matter what kind of filters I apply in Mathematica, I don't get anywhere near the original data (which means the result is several magnitudes off).

I know this is not a MATLAB forum, but for completeness, I put here also my MATLAB code:

numOfElts=128;
x=linspace(1,128,128);
startMat=zeros(numOfElts,numOfElts);
for i = 1:numOfElts
for j = 1:numOfElts
startmat(j,i) = 6+5*cos(2*pi*sqrt((3*(i-numOfElts/2-0.5)/(numOfElts/2))^2+(3*(j-numOfElts/2-0.5)/(numOfElts/2))^2))-(3*(i-numOfElts/2-0.5)/(numOfElts/2))^2-(3*(j-numOfElts/2-0.5)/(numOfElts    /2))^2;
end
end
middleVecOrigin=startmat(64,:);
orthoProj=sum(startmat);
virtualRDN=transpose(repmat(orthoProj,180,1));

middleVecRes=quasiOriginalData(64,:);

figure
subplot(1,3,1)
contourf(startmat)
colorbar

subplot(1,3,2)
plot(x,middleVecOrigin,x,middleVecRes)
subplot(1,3,3)
contourf(quasiOriginalData)
colorbar

• "Flagged it now as bug, because no one responded with an answers since August 2020" - but did you report it to Support? Feb 4, 2021 at 10:56
• The frequency is easy to match, just pad the projectionAlongYDir: makeVirtaulRDNfromSingleProj[projectedDataList_] := Image@Transpose@ ConstantArray[ ArrayPad[projectedDataList, (Sqrt[2] - 1) Length[projectedDataList]/2 // Floor, "Fixed"], Sqrt[2] Length[projectedDataList] // Floor]. Feb 5, 2021 at 5:47
• @J.M. No I didn't report it as bug, since I wanted some interaction from the community first to make sure that, I haven't made a mistake of overlooking something ;-) Feb 5, 2021 at 15:17
• @xzczd: Thanks, that seems to work at least with the frequencies, but it doesn't get the amplitdues right. Feb 5, 2021 at 15:33
• Actually the source code of InverseRadon can be checked with GeneralUtilitiesPrintDefinitions@InverseRadon. (The output is quite readable! ) But I'm not that familiar with the underlying math so it's a bit hard for me to go deeper. Mar 14, 2021 at 7:41

As it explained in tutorials functions Radon and InverseRadon are supposed to be used with images only and not with arbitrary matrix, array or list. Even direct use ImageData[InverseRadon[Radon[Image[data]]]] does not return data as it is, since Radon returns an image in which values are normalized so that the highest possible value is 1. Let see this example (as in Matlab code)

numOfElts = 128; startmat = ConstantArray[0, {numOfElts, numOfElts}];

Do[
Do[
startmat[[j, i]] =
6 + 5*Cos[
2*Pi*Sqrt[(3*(i - numOfElts/2 - 0.5)/(numOfElts/
2))^2 + (3*(j - numOfElts/2 - 0.5)/(numOfElts/
2))^2]] - (3*(i - numOfElts/2 - 0.5)/(numOfElts/
2))^2 - (3*(j - numOfElts/2 - 0.5)/(numOfElts/2))^2;, {j,
1, numOfElts}];, {i, 1, numOfElts}]


First we apply Radon, then InverseRadon and get this images

im0 = Image[startmat]


Therefore we can't return to startmat, and if we try compare cross section we need to use some scaling

ListPlot[{ImageData[im0][[64]],
ImageData[im1][[64]]}]

ListLinePlot[
Rescale /@ {(ImageData@im0)[[64]], (ImageData@im1)[[64]]}]


The best scaling what I got with @xzczd padding is not so differ from Matlab iradon, the corresponding code is given by

n = numOfElts = 2 128; startmat =
ConstantArray[0, {numOfElts, numOfElts}];

Do[
Do[
startmat[[j, i]] =
6 + 5*Cos[
2*Pi*Sqrt[(3*(i - numOfElts/2 - 0.5)/(numOfElts/
2))^2 + (3*(j - numOfElts/2 - 0.5)/(numOfElts/
2))^2]] - (3*(i - numOfElts/2 - 0.5)/(numOfElts/
2))^2 - (3*(j - numOfElts/2 - 0.5)/(numOfElts/2))^2;, {j,
1, numOfElts}];, {i, 1, numOfElts}];
orthoProj = Sum[startmat[[j]], {j, numOfElts}]; max =
Max[orthoProj]; min = Min[orthoProj]; uv =
ConstantArray[1, {Length[orthoProj]}];
makeVirtaulRDNfromSingleProj[projectedDataList_] :=
Image@Transpose@
ConstantArray[
projectedDataList, ((Sqrt[2] - 1) Length[projectedDataList]/2 //
Floor), "Fixed"], Sqrt[2] Length[projectedDataList] // Floor]

Method -> {(# &), "CutoffFrequency" -> .5}] &@orthoProj;

{ListPlot[Rescale /@ {startmat[[n/2]], ImageData[result][[n/2]]},
Joined -> True,
PlotLabel -> "Rescaled cross-sections\nthrough center",
PlotLegends -> {"original function", "after InverseRadon[]"}],

ListPlot[{startmat[[n/2]],
Pi ((max - min) Rescale@(ImageData[result // ImageAdjust][[n/2]]) +
min uv)/n}, Joined -> True,
PlotLabel -> "unscaled cross-sections\nthrough center",
PlotLegends -> {"original function", "after InverseRadon[]"}]}


The final step is to redefine uv so that it minimize difference of original function and result.

• So, you think it's not possible to reproduce the behavior of iradon of MATLAB with InverseRadon? Mar 30, 2021 at 2:29
• @xzczd Actually we can and my example demonstrates how we should prepare image for InverseRadon with Radon method. Also we can't compare final result without scaling since InverseRadon returns image, not array or matrix as iradon in Matlab. Mar 30, 2021 at 10:49
• @AlexTrounev This scaling is in my opinion an unnecessary operation, as even the MMa help for Radon transform uses AdjustImage[]. I need, however, the absolute amplitude of the not known initial distribution. A further problem in my opinion is that the the scale factor depends on the number of sampling points in combination with the fact that for my problem it is necessary to do zero padding to recover the most important frequencies. But if the scale depends on the number of zeros used which changes the frequencies, I have a real problem in efficiently estimating an optimal number of zeros. Mar 30, 2021 at 20:01
• @Quit007 Probably your problem is not linked to InverseRadon directly, since this function has to be applied to images only. It looks like your problem is mostly mathematical problem. In this regard your example with komplFkt also not directly linked to your problem. Mar 30, 2021 at 22:35
• @AlexTrounev Of course komplFkt is only an example function, which, nevertheless, has the same symmetry than my actual problem. Maybe I just fail to see, why InverseRadon in Mma shouldn't be able to recover the absolute values - which is basically the task of inverse radon transforms used in life science e.g. phase contrast computer tomography for which it is essential to calculate the absolute value of the objects. This might be also caused by the fact that all scientific literature I have looked at dealt specifically with radon transforms restoring the absolute values. Mar 30, 2021 at 23:04

As this has come up in the comments, somebody might be interested at least in the discrete version as described here: discrete radon transform:

getRadonLinesComped[] is the (compiled) function that extracts the discrete lines from the image. The uncompiled version is getRadonLines[].

(*rekursive Funktion, die aus einer Steigung und Vektorlänge Daten \
für eine entprechende "digital line" generiert*)
calcSlopeRLAndMiddleElem[structure_(*structure \[Equal] {slope_,
oldVectorSize_(odd Number),...other stuff*)] := {(*total slope in \
the new halves*)(structure[[1]] - Mod[structure[[1]], 2])/2,
(*length of new halves*)(structure[[2]] - 1)/
2,(*Value of the middle element*)Mod[structure[[1]], 2]}

(*erzeugt die Steigungspunkte für eine "digital line" für eine \
gesamtsteigung und eine Bildgröße
ACHTUNG Bildgröße muss sein: 2^Nx2^N *)
constructDLine[totalSlope_, imgWidth_] := Module[{recursiveStruct},
(*Die Rekursion erfolgt solange,
bis die Länge der Subliste 1 ist  \[Rule] #[[2]]>1&.*)
recursiveStruct =
Drop[NestWhileList[
calcSlopeRLAndMiddleElem, {totalSlope,
imgWidth - 1}, #[[2]] > 1 &], 1];
(*structure of recursiveStruct is as follows:
{{total slope in the new halves,length of new halves,
value of the middle element which is not in any one of the \
halves}},...}
*)
(*Das hier würde den parametrischen y-
Verlauf der Linien im Bild ergeben:
Accumulate@Prepend[makeSlopeList[recursiveStruct],0]
*)
(*das hier gibt den relativen Anstieg der Linien im Bild:*)
Prepend[makeSlopeList[recursiveStruct], 0]
]

(*erzeugt aus den Daten in recursiveResultStruct die \
"Steigungspunkte" einer "digital line"*)
makeSlopeList[recursiveResultStruct_] :=
Module[{startingList, elementaryList, middleElements},
startingList = ConstantArray[recursiveResultStruct[[-1, 1]], 3];
elementaryList =
ReplacePart[startingList, 2 -> recursiveResultStruct[[-1, -1]]];
(*extrahiert die Mittenelemente und ordnet sie aufsteigend*)
middleElements =
Reverse[#[[3]] & /@ Drop[recursiveResultStruct, -1]];
(*das hier erzeugt aus den mittenelementen die Gesamtliste*)
Fold[Join[#1, {#2}, #1] &, elementaryList, middleElements]
]

parametricDLines[totalSlope_, imgWidth_, yOffset_] := Module[{},
Transpose[{Range[imgWidth],
1 + Accumulate@constructDLine[totalSlope, imgWidth] + yOffset}]
]

Module[{grayImg, bild0, tempImg, imgWidth, tempImgData, offsetVar,
(*Operationen für den Fall, dass der Input ein Bild ist:
(*das eigentliche Bild wird unten und oben mit einem 0-
grayImg=ColorConvert[#,"Grayscale"]&@img;
bild0=0 grayImg;
tempImg=ImageAssemble[{{bild0},{grayImg}, {bild0}}];
imgWidth=ImageDimensions[img][[1]];
tempImgData=Transpose@ImageData@tempImg;*)
imgWidth = Dimensions[imgData][[1]];

(*Möglichkeit zur Anpassung des Prozesses:
tempImgData=Transpose@Join[0 imgData,imgData,0 imgData];
ist näher am Originalpaper,
verursacht aber bei der Rücktrafo einen zum Rand anteigenden "Kegel".
tempImgData=Transpose@Join[imgData,imgData,imgData];
Verursacht keinen ansteigenden Kegel
*)
(*tempImgData=Transpose@Join[0 imgData,imgData,0 imgData];*)
tempImgData = Transpose@Join[imgData, imgData, imgData];
testImg1,64]] ;
dLineDef=parametricDLines[30,imgWidth,offsetVar];
Image@Reverse@Table[Extract[tempImgData,
dLineDef/.(offsetVar\[Rule]offset)],{offset,2 imgWidth,1,-1}]
*)

dLineDef = parametricDLines[slope, imgWidth, offsetVar];
(*das Extract extrahiert die Datenwerte entlang einer digital \
line. Das Table enthält deshalb diese Linieneinträge zeilenweise.
Die Summe über jede Zeile erfolgt daher mit Total/@*)
Total /@
Table[Extract[tempImgData,
dLineDef /. (offsetVar -> offset)], {offset, 2 imgWidth, 1, -1}]
, {slope, 0, imgWidth - 1}];
(*das Reverse@Transpose@ bewirkt eine Anordnung so,
dass die Radontrafo mit erster Spalte den Winkel 0 hat,
letzte Spalte den Winkel 45\[Degree] und das Bild aufrecht unten \
steht*)
]

Compile[{{imgData, _Real, 2}},
Module[{grayImg, bild0, tempImg, imgWidth, tempImgData, offsetVar,
arrayLength, i, structure, temp, recursiveStruct,
startingList, elementaryList, middleElements, slopeList,
tempSlope, tempSize, tempMiddle,
xVals, rawDlineShape},

imgWidth = Dimensions[imgData][[1]];

tempImgData =(*Transpose@*)Join[imgData, imgData, imgData];
testImg1,64]] ;
dLineDef=parametricDLines[30,imgWidth,offsetVar];
Image@Reverse@Table[Extract[tempImgData,
dLineDef/.(offsetVar\[Rule]offset)],{offset,2 imgWidth,1,-1}]
*)

arrayLength = Log2[imgWidth] - 1;
xVals = Range[imgWidth];
tempSlope = slope;
tempSize = imgWidth - 1;
tempMiddle = tempSlope;
recursiveStruct = ConstantArray[0, {arrayLength, 3}];
i = 1;
While[tempSize > 1 ,
tempMiddle = Mod[tempSlope, 2];
tempSlope = Round[(tempSlope - Mod[tempSlope, 2])/2];
tempSize = Quotient[(tempSize - 1), 2];
recursiveStruct[[i]] = {tempSlope, tempSize, tempMiddle};
i++;
];
(*recursiveStruct wird korrekt aufgebaut!*)

startingList = ConstantArray[recursiveStruct[[-1, 1]], 3];
elementaryList = {startingList[[1]], recursiveStruct[[-1, -1]],
startingList[[1]]};
(*extrahiert die Mittenelemente und ordnet sie aufsteigend*)
middleElements = Reverse[#[[3]] & /@ Drop[recursiveStruct, -1]];
(*das hier erzeugt aus den mittenelementen die Gesamtliste*)
slopeList =
Fold[Join[#1, {#2}, #1] &, elementaryList, middleElements];
rawDlineShape = 1 + Accumulate@Prepend[slopeList, 0];

Total /@ Table[

tempImgData[[rawDlineShape[[z]] + offset,
xVals[[z]]]], {offset, 2 imgWidth, 1, -1}, {z, imgWidth}]
, {slope, 0, imgWidth - 1}];
(*
(*das Extract extrahiert die Datenwerte entlang einer digital \
line. Das Table enthält deshalb diese Linieneinträge zeilenweise.
Die Summe über jede Zeile erfolgt daher mit Total/@*)
Total/@Table[Extract[tempImgData,
Transpose[{Range[imgWidth],1+ Accumulate@constructDLineComped[
slope,imgWidth]+offsetVar}]
],{offset,2 imgWidth,1,-1}]
,{slope,0,imgWidth-1}];
(*das Reverse@Transpose@ bewirkt eine Anordnung so,
dass die Radontrafo mit erster Spalte den Winkel 0 hat,
letzte Spalte den Winkel 45\[Degree] und das Bild aufrecht unten \
steht*)
], RuntimeAttributes -> {Listable}, Parallelization -> True,
CompilationTarget -> "C"];


the Radon data from an image can be determined then by this function, where imgData is the matrix containing the image data:

backProjection[radonData_] :=
Module[{imgWidth, height, transRDNdata, dLineDef, orthogonalProj,
dLineJumps, backProjectedData, slope},

backProjectedData = Sum[
(*das hier bestimmt die Orte,
an denen die Matrix getiltet werden muss*)
dLineDef = parametricDLines[slope, imgWidth, 0];
(*das hier erzeugt eine Seed Matrix aus den zeilenweisen \
Der Zeileneintrag wird hier verschmiert und deshalb müsste der \
Wert im "Schmierer" 1/imgWidth betragen,
aber diese Operation wird später durchgeführt.*)
orthogonalProj = ConstantArray[transRDNdata[[slope]], imgWidth];
dLineJumps = Transpose[dLineDef][[2]];
Do[orthogonalProj =
MapAt[RotateRight[#, dLineJumps[[i]] - 1] &, orthogonalProj,
i], {i, 2, imgWidth}];
orthogonalProj
, {slope, imgWidth}];
1/imgWidth^2 Transpose@(Take[#, -imgWidth] & /@ backProjectedData)
]

(*erzuegt die verschiedenen Geometrien,
die für die diskrete Radon Trafo mittels digital lines benötigt \
werden*)
(*falls der Input ein Bild ist,
kann das durch folgenden Code bewerkstelligt werden:
ImageReflect@#}&@img;
*)
(*falls der Input ein DatenArray ist,
wird das durch folgenden Code bewerkstelligt:*)
Reverse@Transpose@Reverse@#} &@(imgData);
]


The backtransform function is this one:

discreteRadonBackProject[radonDataSet_] :=
rawBackProjectionData, correctedBPdata, finalBP},
(*erzuegt die verschiedenen Geometrien,
die für die diskrete Radon Trafo mittels digital lines benötigt \
werden*)
(*falls der Input ein Bild ist,
kann das durch folgenden Code bewerkstelligt werden:
ImageReflect@#}&@img;
*)
(*averageImgBrightness=Mean[Flatten[imgData]];
(*falls der Input ein DatenArray ist,
wird das durch folgenden Code bewerkstelligt:*)
Reverse@#}&@(imgData-0 averageImgBrightness);
*)

(*ArrayPlot[#,PlotLegends\[Rule]Automatic]&/@
rawBackProjectionData*)

correctedBPdata = {#[[1]], Transpose@Reverse@#[[2]], Reverse@#[[3]],
Reverse@Transpose@Reverse@#[[4]]} &@rawBackProjectionData;
(*ArrayPlot[#,PlotLegends\[Rule]Automatic]&/@correctedBPdata*)

finalBP = 1/4  Total[correctedBPdata];

finalBP
]


The exact inverse is then found by a multi-grid approach with these functions:

shrinkArrays[arrayList_] :=
ArrayResample[#, 1/2 Dimensions[#]] & /@ (1./4 arrayList);
enlargeArray[array_] := ArrayResample[array, 2 Dimensions[array]];

sharpening[array_] :=
ListConvolve[{{-1./16, -1./8, -1./16}, {-1./8,
3./4, -1./8}, {-1./16, -1./8, -1./16}}, array, 1, 0]

Module[{i = 1, numOfDownRecursions, initialBackTransform,
shrunkResidual, forwardList, approxSolutions, finalApproxSol},
forwardList = {initialBackTransform};
While[i <= numOfDownRecursions - 2,
forwardList[[1]]]];*)
shrunkResidual =
PrependTo[forwardList,
i++;}
];
h*)},numOfDownRecursions-2];
approxSolutions=Reverse@Transpose[forwardList][[2]];
finalApproxSol=Fold[#2+enlargeArray[#1]&,approxSolutions]*)
Fold[#2 + enlargeArray[#1] &, forwardList]
]

(*approximation by iterative comparison of \
Module[{initialBackTransform, resList},
resList =
NestList[# +
numOfIterations];
resList
]

convergenceCheck[initialImgData_, improvementDataList_] :=
StandardDeviation[Flatten[(-initialImgData + #)]] & /@
improvementDataList


For testing I have used these functions with this image:

It is necessary that the actual image is just single channel (e.g. grayscale) and the size is NxN with N~2^k with k an integer, e.g. defined as: ImageData@ImageResize[ColorConvert[(*image here*), "Grayscale"],{128,128}]

DRTC256 = discreteRadonC[startImgData2];
ArrayPlot[#, PlotLegends -> Automatic] & /@ DRTC256


400 iterations are calculated for the approximate backtransform and every 80th is plotted

t1 = AbsoluteTime[];

• I'm now with a mobile phone so can't check further, but what if you use Table instead of ConstantArray in Compile? Apr 3, 2021 at 2:30
• Function backProjection is not defined. Also function getRadonLinesComped is not really compiled to "C" due to several errors. What version you are run on? Apr 3, 2021 at 7:28
• This code not solves the problem of 2-dimensional radially symmetric distribution reconstruction from a single projection calculation and also it has no any advantage against Radon and InverseRadon`. Nevertheless it is good that Quit007 made this efforts (+1). May be we need next bounty with 1000 points to solve this problem:) Apr 11, 2021 at 22:20