InverseRadon behaves differently from iradon of MATLAB

Question

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

makeVirtaulRDNfromSingleProj[projectionAlongYDir] // ImageAdjust

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

"virtual" Radon transform

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));
quasiOriginalData = iradon(virtualRDN,0:179,128);

middleVecRes=quasiOriginalData(64,:);

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

subplot(1,3,2)
plot(x,middleVecOrigin,x,middleVecRes)
legend('original cross-section','inverse radon c.-s.')
subplot(1,3,3)
contourf(quasiOriginalData)
colorbar

Answer 1

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]
im0R = ImageResize[Radon[im0], {128, 128}]
 im1 = ImageResize[InverseRadon[im0R], {128, 128}]

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[
    ArrayPad[
     projectedDataList, ((Sqrt[2] - 1) Length[projectedDataList]/2 // 
       Floor), "Fixed"], Sqrt[2] Length[projectedDataList] // Floor]

result = InverseRadon[makeVirtaulRDNfromSingleProj[#], {n, n}, 
     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.

Answer 2

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

getRadonLines[imgData_] := 
 Module[{grayImg, bild0, tempImg, imgWidth, tempImgData, offsetVar, 
   dLineDef, offset, radonData},
  (*Operationen für den Fall, dass der Input ein Bild ist:
  (*das eigentliche Bild wird unten und oben mit einem 0-
  Bild gepaddet*)
  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];
  (*Funktionstest mit: \[Rule] Input: getRadonLines[ImageResize[
  testImg1,64]] ;
  dLineDef=parametricDLines[30,imgWidth,offsetVar];
  Image@Reverse@Table[Extract[tempImgData,
  dLineDef/.(offsetVar\[Rule]offset)],{offset,2 imgWidth,1,-1}]
  *)

  radonData = Table[
    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*)
  (*1/imgWidth*) Reverse@Transpose@radonData
  ]

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

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

   tempImgData =(*Transpose@*)Join[imgData, imgData, imgData];
   (*Funktionstest mit: \[Rule] Input: getRadonLines[ImageResize[
   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];
   radonData = Table[
     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}];
   Reverse@Transpose@radonData
   (*
   radonData=Table[
   (*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*)
   (*1/imgWidth*) Reverse@Transpose@radonData*)
   ], 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},
  {height, imgWidth} = Dimensions[radonData];
  transRDNdata = Transpose[radonData];

  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 \
Projektionen der Radon-Transformierten des Objekts.
    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)    
  ]


discreteRadonC[imgData_] := 
 Module[{averageImgBrightness, radonInitImgs, radonInitData, 
   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:
  radonInitImgs={#,ImageRotate@#,ImageReflect@#,ImageRotate@
  ImageReflect@#}&@img;
  radonDataSet=getRadonLines/@ImageData/@radonInitImgs;
  *)
  (*falls der Input ein DatenArray ist, 
  wird das durch folgenden Code bewerkstelligt:*)
  radonInitData = {#, Reverse@Transpose@#, Reverse@#, 
      Reverse@Transpose@Reverse@#} &@(imgData);
  (*ArrayPlot[#,PlotLegends\[Rule]Automatic]&/@radonDataSet*)
  radonDataSet = ParallelMap[getRadonLinesComped, radonInitData];
  radonDataSet
  ]

The backtransform function is this one:

discreteRadonBackProject[radonDataSet_] := 
 Module[{averageImgBrightness, radonInitImgs, radonInitData, 
   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:
  radonInitImgs={#,ImageRotate@#,ImageReflect@#,ImageRotate@
  ImageReflect@#}&@img;
  radonDataSet=getRadonLines/@ImageData/@radonInitImgs;
  *)
  (*averageImgBrightness=Mean[Flatten[imgData]];
  (*falls der Input ein DatenArray ist, 
  wird das durch folgenden Code bewerkstelligt:*)
  radonInitData={#,Reverse@Transpose@#,Reverse@#,Reverse@Transpose@
  Reverse@#}&@(imgData-0 averageImgBrightness);
  radonDataSet=getRadonLines/@radonInitData;
  (*ArrayPlot[#,PlotLegends\[Rule]Automatic]&/@radonDataSet*)
  *)

  rawBackProjectionData = backProjection /@ radonDataSet;
  (*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];
  (*{ArrayPlot/@radonDataSet,ArrayPlot/@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]

approxRadonInverseByFullMultiGrid[radonData_] := 
 Module[{i = 1, numOfDownRecursions, initialBackTransform, 
   shrunkResidual, forwardList, approxSolutions, finalApproxSol},
  numOfDownRecursions = Log2[Dimensions[radonData[[1]]][[2]]];
  initialBackTransform = discreteRadonBackProject[radonData];
  forwardList = {initialBackTransform};
  shrunkResidual = radonData;
  While[i <= numOfDownRecursions - 2,
   {(*shrunkResidual=shrinkArrays[shrunkResidual-discreteRadon[
    forwardList[[1]]]];*)
    shrunkResidual = 
     shrinkArrays[shrunkResidual - discreteRadonC[forwardList[[1]]]];
    PrependTo[forwardList, 
     sharpening@discreteRadonBackProject[shrunkResidual]];
    i++;}
   ];
  (*forwardList=NestList[{shrinkArrays[#[[1]]-discreteRadon[#[[2]]]],
  sharpening@discreteRadonBackProject[shrinkArrays[#[[1]]-
  discreteRadon[#[[2]]]]]}&,{radonData(*f^h*),initialBackTransform(*u^
  h*)},numOfDownRecursions-2];
  approxSolutions=Reverse@Transpose[forwardList][[2]];
  finalApproxSol=Fold[#2+enlargeArray[#1]&,approxSolutions]*)
  Fold[#2 + enlargeArray[#1] &, forwardList]
  ]

(*approximation by iterative comparison of \
radon[radonInversed[approximation]] to original radon data*)
approxRadonInverseFMG[radonData_, numOfIterations_] := 
 Module[{initialBackTransform, resList},
  initialBackTransform = approxRadonInverseByFullMultiGrid[radonData];
  resList = 
   NestList[# + 
      approxRadonInverseByFullMultiGrid[
       radonData - discreteRadonC[#]] &, initialBackTransform, 
    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[];
iterTest64 = approxRadonInverseFMG[DRTC64, 400];
t2 = AbsoluteTime[];
Print[(t2 - t1)/3600. "h"];
ArrayPlot[#, PlotLegends -> Automatic] & /@ 
 iterTest64[[1 ;; 400 ;; 80]]

Depending on the size of the initial image size there is a iteration number dependent minimum of the residual error. Unfortunately, with this discrete version (besides the fact that this is quite slow) I wasn't able to calculate the virtual (discrete) radon transform of the single backprojection of my radially symmetric problem I have described in my above question. And that's why the discrete transform doesn't really work for my problem.