The function implements three of the major iterative reconstruction techniques: ART (Algebraic Reconstruction Technique), EM (Likelihood Reconstruction using Expectation Maximization) and LSCG (Least Squares Conjugate Method).
iradonIT(rData, XSamples, YSamples, StartImage = "None",
mode = "EM", UseFast = 1, RadonKernel = "NN",
Iterations = 20, IterationsType = "random",
SaveIterations = 0, SaveIterationsName = "",
LowestALevel = 0, ConstrainMin = -1,
ConstrainMax = -1, Alpha = 1, Beta = 1,
Regularization = 0, KernelFileSave = 0,
KernelFileName = "", RefFileName = "None",
ThetaSamples = nrow(rData), RhoSamples = ncol(rData),
ThetaMin = 0, RhoMin = -0.5*((2*round(sqrt(XSamples^2+
YSamples^2)/2)+1)-1), DeltaTheta = pi/ThetaSamples,
DeltaRho = (2*abs(RhoMin)+1)/RhoSamples,
Xmin = -0.5*(XSamples-1), Ymin = -0.5*(YSamples-1),
DeltaX = 1, DeltaY = 1, OverSamp = 0,
DebugLevel = "Normal", iniFile = NULL)(matrix) A matrix that contains the sinogram. The rows of the image correspond to the sampled angles \(\theta\) and the columns correspond to the sampled distance \(\rho\). These both parameters determine the lines.
(integer) XSamples is the number of samples on the x-axis (rows) in the reconstructed image.
(integer) YSamples is the number of samples on the y-axis (columns) in the reconstructed image.
(matrix) If provided, then the initial guess on the reconstructed image is taken from this image, otherwise a constant solution will be assumed from the start. Default is set to StartImage="None". Note, that the dimension of the StartImage has to be dim(StartImage)==c(XSamples,YSamples).
(character) The iterative reconstruction method. Default is mode="EM" for Likelihood Reconstruction using Expectation Maximization. Additional implementations are "ART" (Algebraic Reconstruction Technique) and "CG" (Least Squares Conjugate Method).
(logical) If UseFast=1 fast reconstruction is used, where the system matrix is stored (in case of KernelFileSave=1) using sparse techniques, otherwise a slower but more efficient memory reconstruction is used. Default is UseFast=1.
(character) The type of kernel is used to model the system matrix. Default is RadonKernel="NN" for two-level Nearest Neighbour approximation (Memory consuming). Additional implementations are "RNN" for ray driven Nearest Neighbour approximation discrete Radon transformation based (very fast with small system matrix), "RL" for ray driven Linear Interpolation discrete Radon transformation based (fast with small system matrix), "P1" for method based on Radon transformation of square with pre-guidance (slow but good) and "P2" for method based on Radon transformation of square with no pre-guidance (slower but better).
(integer) Using mode="EM" or "CG", it is the number of iterations before the iteration ends. Using mode="ART", it is the number of full iterations, i.e., divided by the number of rows in the system matrix. Defaults to Iterations=20.
(character) Using mode="ART", two iterations types are possible:
"cyclic" selection of the row index or "random" selection. Defaults to IterationsType = "random".
(integer) If it is set to SaveIterations=1, the current solution will be saved after each iteration under the name SaveIterationsName + CurrentIteration (also possible 2,3 ...). Defaults to SaveIterations=0.
NOTE: When Iterations/itINI.SaveIterations>300 then too many new files are requested. In this case default is used.
(character) In case of SaveIteration is non-zero, the parameter determines the name of the currently saved iteration. It contains the path to the file and the filename. The path contains the path relatively to the working-directory of your R-session or it contains the full path of the file. Defaults to SaveIterationsName="".
(double) If fast reconstruction is used, the matrix elements will be truncated to this level relatively to the sampling distance of \(x\). Defaults to LowestALevel=0.0.
(double) After each iteration the solution in each sample will be forced above this limit. Defaults to ConstrainMin=0.
(double) After each iteration the solution in each sample will be forced below this limit. Defaults to ConstrainMax=-1.
NOTE: For both limit it is assumed that negative limits imply that the feature is not used.
(double) If mode="ART", it will be the initial update weight. Defaults to Alpha=1.0.
(double) If mode="ART", it will be the multiplicative change to the weight factor, which should be less than one. Default is set to Beta=1.
(double) If the Regularization is not zero and positive and using fast reconstruction, rows will be appended to the system matrix with a simple Laplace operator. A weight factor should also be incorporated. If negative, identity matrix will be used. Defaults to Regularization=0.
(logical) If KernelFileSave=1 and fast reconstruction is used, the system matrix will be saved under the name KernelFileName. Defaults to KernelFileSave=0.
(character) If fast reconstruction is used, the system matrix will be saved and restored with this sif-name. If the system matrix is incompatible to the sampling parameters, a new one will be generated. KernelFileName contains the path to the filename and the filename. The path contains the path relatively to the working-directory of your R-session or it contains the full path to the file. Defaults to KernelFileName="".
(character) Error measures can be made between the image which is contained in RefFileName and the reconstructed image from rData. The file has to be either a fif-file or a pet-file (type ?writeData or ?readData in R to get more information about these formats). If RefFileName exists, the measures will be logged in a file with the name of RefFileName and extension 'dif'. Defaults to RefFileName="None".
(integer) Number of angular samples \(\theta\) in the sinogram.
Defaults to ThetaSamples=nrow(rData).
(integer) Number of samples in the sinogram \(\rho\) in the \(\rho\)-direction. Defaults to RhoSamples=ncol(rData).
(double) Start of the angular sampling (should be set to zero). Defaults to ThetaMin=0.
(double) Start of sampling positions in \(\rho\).
Defaults to RhoMin=-0.5* ((2*round(sqrt(XSamples^2+YSamples^2)/2)+1)-1).
(double) Angular sampling distance. Should be set to the default DeltaTheta = pi/ThetaSamples.
(double) Sampling distance in \(\rho\). Defaults to DeltaRho=(2*abs(RhoMin)+1)/RhoSamples.
(double) The minimum \(x\)-position of the reconstructed image. Defaults to Xmin=-0.5*(XSamples-1).
(double) The minimum \(y\)-position of the reconstructed image. Defaults to Ymin=-0.5*(YSamples-1).
(double) Sampling distance on the \(x\)-axis. Defaults to DeltaX=1.
(double) Sampling distance on the \(y\)-axis. Defaults to DeltaY=1.
(integer) Uses oversampling - squared number of samples are used. Defaults to OverSamp=0.
(character) This parameter controls the level of output. Defaults to DebugLevel="Normal" for a standard level output. Alternative implementations are "Detail" if it is desirable to show almost all output on screen or "HardCore" for no information at all.
(character) If iniFile!=NULL, iniFile has to be the name of an ini-file including a pathname to the file. In the case of a specified iniFile all parameters are read from the file. Note that in this case contingently setting parameters (except for DebugLevel) in R are ignored when calling iradonIT. The parameters which are not specified in iniFile are set to defaults. Default of iniFile is iniFile=NULL.
A matrix, that contains the reconstructed image (matrix) of rData.
A list of following values:
The dimension of the irData.
The minimum x- and y-position in irData.
Sampling distance on the x- and y-axis in irData.
Arguments of the call to iradonIT.
The function implements three of the major iterative reconstruction techniques: ART (Algebraic Reconstruction Technique), EM (Likelihood Reconstruction using Expectation Maximization) and LSCG (Least Squares Conjugate Method). The reconstruction methods were developed by P. Toft (1996) and implemented in R by J. Schulz (2006). For detailed theoretical information about the Radon-transformation see references.
There are a lot of different parameters which influence the reconstruction method. In common cases, it should be enough to determine rData, XSamples and YSamples. For the default parameters, it is assumed that the sinogram is standardised to a coordinate system with \(\theta\) from \(0\) to \(\pi\) and \(\rho\) from \(-0.5*((2*round(sqrt(M^2+N^2)/2)+1)-1)\) to \(0.5*((2*round(sqrt(M^2+N^2)/2)+1)-1)\). It is also assumed that the reconstructed image is standardised to a coordinate system with \(x\) from \(-0.5*(M-1)\) to \(0.5*(M-1)\) and \(y\) from \(-0.5*(N-1)\) to \(0.5*(N-1)\). \(R\), \(M\) and \(N\) are the number of sampling parameters in \(\rho\), \(x\) and \(y\). If that is not true, you will have to adapt the parameters.
Several things must be fulfilled to ensure a reasonable performance. Firstly sampling must be adequate in all parameters (see to references to get detail information). This will imply bounds on the sampling intervals. Secondly it is assumed that the fundamental function \(f(x,y)\) to be reconstructed have compact support, or more precisely is zero if \(\sqrt{x^2+y^2} > |\rho_{max}|\). This demand will ensure that \(Rf(\rho,\theta)=0\), if \(|\rho|>|\rho_{max}|\). If this cannot be fulfilled, numerical problems must be expected.
Assuming that \(f(x,y)\) has compact support, then \((x,y)=(0,0)\) should be placed to minimize \(|\rho_{max}|\). This will reduce the size of the data array used for the discrete Radon transform.
The numerical complexity of the procedure is mainly determined by Iterations, the size of rData, XSamples, YSamples and by the use of fast reconstruction method.
Fast reconstruction is used in case of UseFast=1. It is possible to store the computed system matrix (KernelFileSave=1) under the name KernelFileName. The speed of the reconstruction can obviously be increased by using of an existing system matrix. NOTE: This is only possible, if the sampling parameters are compatible with the system matrix. That means, that the fifteen parameters XSamples, YSamples, DeltaX, DeltaY, Xmin, Ymin, ThetaSamples, RhoSamples, DeltaRho, DeltaTheta, ThetaMin, RhoMin, LowestALevel, RadonKernel and OverSamp are equal to the parameter from which the system matrix is generated. If the system matrix is incompatible to the sampling parameters, a new one will be generated.
Another variation to determine the parameter is the offer with iniFile. iniFile has to be the name of a file (e.g. iniFile="/home/work/sino.ini") with the following structure. Each line contains at first the name of a parameter, then an equal sign follows and the value of the parameter. The first line must contain the parameter mode. Characters are not written in "", e.g. RadonKernel=NN and not RadonKernel="NN". Furthermore note that in an ini-file rData and StartImage are to be of type character, videlicet the name of the corresponding file. Supported file formats are ".txt", ".dat", ".fif", ".pet", ".tif", ".tiff", ".pgm", ".ppm", ".png", ".pnm", ".gif", ".jpg", ".jpeg". See ?readData to get detailed information about supported formats.
If a file of the type ".fif", ".pet" or ".dat" is specified for rData (see ?writeData or ?readData in R to get more information about these formats), the parameters ThetaSamples, RhoSamples, ThetaMin, RhoMin, DeltaTheta and DeltaRho will be read from the file-header of rData, but only in case of unspecified corresponding parameters. Just as, if a file of type ".fif", ".pet" or ".dat" is specified for StartImage, then the parameters XSamples, YSamples, XMin, YMin, DeltaX and DeltaY will be read from the file-header.
Note that in case of an ini-file, contingently setting parameters (except for DebugLevel) in R are ignored when calling iradonIT. Parameters that are not specified in the iniFile are set to defaults.
Toft, Peter, Ph.D. Thesis, The Radon Transform - Theory and Implementation, Department of Mathematical Modelling Section for Digital Signal Processing, Technical University of Denmark, 1996. http://eivind.imm.dtu.dk/staff/ptoft/ptoft_papers.html
Schulz, Joern, Diploma Thesis: Analyse von PET Daten unter Einsatz adaptiver Glaettungsverfahren, Humboldt-Universitaet zu Berlin, Institut fuer Mathematik, 2006.
# NOT RUN {
#
# Compare the results of iterative reconstruction method "EM" and
# direct reconstruction method "FB"
#
# }
# NOT RUN {
P <- phantom(design="B")
rP <- markPoisson(P, nSample=1600000 )
irP1 <- iradon(rP$rData , nrow(P), ncol(P))
irP2 <- iradonIT(rP$rData , nrow(P), ncol(P))
viewData(list(P, rP$rData, irP1$irData, irP2$irData),
list("Generated unnoisy Phantom", "Generated PET Data",
"Direct rec.: mode='FB'", "Iterative rec.: mode='EM'"))
rm(irP1,irP2,P,rP)
# }
# NOT RUN {
#
# Compare the results from the iterative reconstruction methods "EM"
# "CG" and "ART"
#
# }
# NOT RUN {
P <- phantom(n=151, addIm="blurred1")
rP <- markPoisson(P, nSample=2000000, RhoSamples=401)
irP1 <- iradonIT(rP$rData , nrow(P), ncol(P), Iterations=20,
ConstrainMin=0, ConstrainMax=10)
irP2 <- iradonIT(rP$rData , nrow(P), ncol(P), mode="CG",
Iterations=4, ConstrainMin=0, ConstrainMax=10)
irP3 <- iradonIT(rP$rData , nrow(P), ncol(P), mode="ART",
Iterations=5, ConstrainMin=0, ConstrainMax=10,
Alpha=0.1, Beta=0.5)
viewData(list(P,irP1$irData,irP2$irData,irP3$irData),
list("Generated unnoisy Phantom",
"mode='EM', Iterations=20","mode='CG', Iterations=4",
"mode='ART', Iter.=20, Alpha=0.1, Beta=0.5"))
rm(irP1,irP2,irP3,P,rP)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab