Author: Theo Fuchs
All simulations done within the FORBILD projects have to include a realistic and comparable noise level of attenuation values. Thereby comparability is more important than the absolute level of the noise. The algorithm listed in the appendix should be implemented without changes into the image reconstruction software, e.g. directly into the routines for reading noise-free raw data from the disk.
The resulting image noise depends on various factors: for the image reconstruction algorithm and the object, we will assume the following simple case:
The dependency of noise on the x-ray tube power, time of rotation, slice thickness, tube voltage, prefiltration etc. is implicitly taken into account by determining the number of incident quanta n0 of a single x-ray in each of the projections. This number is valid only for the assumptions made above.
The number of incident quanta is found by evaluation of a simulated image of a water phantom like the one defined above. The standard deviation (in HU) of pixel noise within a circular ROI (1cm radius) relates to n0 as:
with (cf. Fig).
It follows, that
Thereby noise for the attenuation values
is calculated by adding gaussian distributed noise
to the mean attenuated number of quanta
We assumed the variance of to equal the mean . is the realization of a normal distribution with zero mean and unit variance.
Figure: Evaluation of a central ROI with 1cm radius yields standard deviation s = 12.38HU (pixel noise). The simulation of raw data was performed with 2.244× 105 incident quanta. A Shepp-Logan convolution kernel was used for reconstruction.
Appendix: The following algorithm adds noise to given attenuation values
by using a gaussian random generator
The case of resulting intensity <= 0 causes numerical problems and makes no physical sense, so we use a simple trick: We throw the dice again until the resulting intensity is positiv.
double sigmaHU; // user defined image noise in HU
double k=3.44e+7; // thi is by definition our constant
// R() is a gaussian distributed randomizer
// with zero mean and standard deviation 1/sqrt(n0)
// from n0=k/sigmaHU^2 we get 1/sqrt(n0)=sigmaHU/sqrt(k)
RandomGauss R(0, sigmaHU/sqrt(k));
double add_noise(double p)
p = exp(-p);
noisy_p = p + sqrt(p)*R();
while(noisy_p <= 0);