The function implements the Acceptance Rejection (AR) procedure to simulate a marked Poisson process with spatial varying intensity. Therewith it is possible to create noisy data from a phantom, as generated in Positron Emission Tomography (PET).
markPoisson(DataInt, nSample = 200000, ThetaSamples = 181,
RhoSamples = 2*round(sqrt(sum((dim(DataInt))^2))/2)+1,
RhoMin = -0.5*sqrt(2),
vect.length = 100000, image = TRUE, DebugLevel = "Normal")(matrix) A two dimensional image in which the matrix elements comply with intensities. That means the intensity of the decay of positrons in a certain tissue.
(integer) nSample determines the number of accepted events that will be generated with AR method. Defaults to nSample = 200000.
(integer) Specifies the number of samples in the angular parameter \(\theta\) in the sinogram. The sinogram is sampled linearly from 0 to (approximately) \(\pi\) radians. Defaults to ThetaSamples=181.
(integer) Specifies the number of samples in the distance parameter \(\rho\) in the sinogram. Defaults to RhoSamples=2*round(sqrt(sum((dim(DataInt))^2))/2)+1.
(double) Specifies the minimum sample position in the sinogram on the second axis. It is assumed that DataInt is scaled to \([-0.5,0.5]^2\). Defaults to RhoMin=-0.5*sqrt(2).
(integer) Determines a bound of number of generated accepted events in each iteration. That means, if nSample > vect.length then in each iteration vect.length/(0.1217) events will be generated and roughly vect.length accepted events. If you scale down vect.length the computation will be more time-consuming and if you increase the value of vect.length, more memory will be necessary. Defaults to vect.length = 100000.
(logical) If image=FALSE, only a generated sinogram will be returned, i.e. a noisy Radon transformation of the phantom will be returned. Otherwise, if image=TRUE an additional image will be returned, where each pixel of the image corresponds to a number of accepted events generated with the AR method. Defaults to image=TRUE.
(character) This parameter controls the level of output. Available are: The default DebugLevel="Normal" for standard level output or alternative "Detail" if it desirable to logged almost all output to screen or "HardCore" for no information at all.
A matrix, that contains the Radon transformed data.
Will be returned only in case of image = TRUE. A matrix, where each pixel of the image corresponds to a number of accepted events generated with the AR method.
A list of following values:
The dimension of the rData.
The minimum x- and y-position in rData.
Sampling distance on the x- and y-axis in rData.
Arguments of the call to markPoisson.
The function implements the Acceptance Rejection (AR) method to simulate a marked Poisson process with spatial varying intensity. The function was developed to simulate data of a PET scanner. Therefore, new random events are generated with the AR method from DataInt. An angle \(\theta\) will be generated to each event (uniform and iid. in \([0,\pi]\)), because a PET scanner detects two photons, who will travel in (nearly) opposite directions. The line between the two detectors has the line parameters \((\rho,\theta)\), whereas \(\rho\) is the shortest distance from the origin of the coordinate system to the line, and \(\theta\) an angle corresponding to the angular orientation of the line. The only obtainable information from the two photons is the fact, that the photon emission took place somewhere along that line.
Schulz, Joern, Diploma Thesis: Analyse von PET Daten unter Einsatz adaptiver Glaettungsverfahren, Humboldt-Universitaet zu Berlin, Institut fuer Mathematik, 2006.
Gentle, J.E., Elements of Computational Statistics, Springer-Verlag New York Berlin Heidelberg, 2002.
# NOT RUN {
P <- phantom()
rP <- radon(P)
mP1 <- markPoisson(P)
mP2 <- markPoisson(P, nSample = 1000000)
viewData(list(P, mP1$Data, mP2$Data, rP$rData, mP1$rData, mP2$rData),
list("Phantom", "nSample = 200000", "nSample = 1000000",
"Radon Transfom of Phantom", "nSample = 200000",
"nSample = 1000000"))
cat("Number of generated accepted events for mP2:",sum(mP2$Data),"\n")
rm(mP1,mP2,P,rP)
# }
Run the code above in your browser using DataLab