Simulation of
NoiseAuthor: 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:

- 20cm Water phantom located at the center of the field of measurement.
- Monochromatic simulation of the attenuation integrals.
- Fan beam geometry of a Somatom Plus4 A scanner (Siemens): RZ=1, 768 detectors within each of 1056 projections.
- Shepp-Logan convolution kernel.
- Image reconstruction on a PCTomo Station (MIR) with 20cm field of view.

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 *n*_{0} 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 *n*_{0} 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× 10 ^{5} 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

*RandomGauss(*Mean, Sigma):

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)

{

double noisy_p;

p = exp(-p);

do

{

noisy_p = p + sqrt(p)*R();

}

while(noisy_p <= 0);

return -log(noisy_p);

}