MLEM and normalization factors

General Discussions
1

Dear CASToR users and developers,

I am having trouble understanding something related to MLEM and how normalization factors are applied in CASToR.

Basically, the MLEM update steps looks like this:

Update = 1 / Sensitivity × Backprojection(Data / (Projection(Estimate) + Random + Scatter))

If I understand correctly, normalization factors that are present in the input datafile intervene twice in this equation:

  1. In the Projection step (oProjectionLine.cc:759)
  2. In the Backprojection step (oProjectionLine.cc:864)

We are currently working on simulation data which do not contain any random or scatter events. When we do the math from the equation, in which we set both Random and Scatter to 0, we reach the conclusion that normalization factors cancel out in the projection-backprojection sequence.

In order to confirm this, we added some print statements into CASToR’s source code and indeed, the normalization values seem to cancel out. The only way to influence results with normalization factors without random or scatter seem to be through the sensitivity image.

Have we reached correct conclusions? If yes, do you have any reference that would hint at why the equation is designed in such a way? I find it somewhat counterintuitive.

Thank you very much for any piece of information that could help us understand what is going on!

Best,

Aurélien.

PS: is there a way to enable math rendering on the forum?

2

Hi Aurélien,
The ML-EM updates derives from the Poisson likelihood function. You could find its derivation in the following book (section MLEM and OSEM), but there are multiple other references : Positron Emission Tomography: Basic Sciences | SpringerLink

The normalization (and attenuation) factors cancel out only because there is no background events in your simulation data.

Regarding math formatting, actually there is ! I just didn’t activate the math plugin.
Now $ or $$ can be used for math rendering in the forum (not sure it works well with email though):

Example:

$$
x_j^{k+1} = \frac{x_j^k}{\sum_i^n P_{ij}} \sum_i^n \frac{P_{ij} \space m_i}{ \sum_j^n P_{ij} \space x_j^k} 
$$

becomes

x_j^{k+1} = \frac{x_j^k}{\sum_i^n P_{ij}} \sum_i^n \frac{P_{ij} \space m_i}{ \sum_j^n P_{ij} x_j^k}

Best,
Thibaut

1 Like