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.
projectionAlongYDir
:makeVirtaulRDNfromSingleProj[projectedDataList_] := Image@Transpose@ ConstantArray[ ArrayPad[projectedDataList, (Sqrt[2] - 1) Length[projectedDataList]/2 // Floor, "Fixed"], Sqrt[2] Length[projectedDataList] // Floor]
. $\endgroup$InverseRadon
can be checked withGeneralUtilities`PrintDefinitions@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. $\endgroup$