Scatter rate estimation for LM-TOF-MLEM

[jbaran1]

Dear CASToR Users and Developers,

me and my colleagues found an issue with scatter correction implemented in CASToR (based on simulated data). I will start from the beginning.

We performed a simulation with GATE (NEMA IEC phantom) for the long axial FOV PET scanner and tried to reconstruct the results (true+scatter events) with scatter correction information incorporated. We estimated the scatter rate for each TOF bin independently (it is about 1-2 per TOF bin). We also specify the size of the TOF bin in the header. For the comparison, as the reference, we also reconstructed 2 other cases:
- only true coincidences
- true+scatter coincidences without any correction
For all three cases we use the same TOF resolution. All the images are done with LM-MLEM and raytracing projector. All the reconstructed images are attenuation corrected. You can find in the attached plot the reconstructed images from 20th iteration (NEMA_IEC_factor1.png):
TOP LEFT: true + scatter coincidences scatter corrected
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true coincidences
BOTTOM RIGHT: profiles along the orange line from previous plots

Unfortunately, the images reconstructed with true+scatter coincidences with scatter correction do not look like those reconstructed only with true. They are more like the one reconstructed with true+scatter events. Even the values are at the same level as for the images reconstructed with true+scatter events. Moreover, we also dig a little bit in the CASToR source code and dump forward projection values and it seems that the factor of the image estimation is three levels of magnitude greater than the additive correction factor. It suggests that probably our scatter correction does not work correctly.

So we started to play around with scatter correction rates. We did some other exercises and multiplied our scatter rate by several factors ranging from 1 to 10E6. I attached the results for the multiplicative factor equal to 20000 (NEMA_IEC_factor20000.png). The legend is as for the previous image. The results with the scatter correction applied start to look like images reconstructed with true coincidences. However the scatter effect (increased signal in the middle of NEMA IEC) inside the phantom is visible. Image quality (factor_20000 images) at the boundaries of the scanner is better. Additionally, the overall voxel units are still at the level of the voxels of images reconstructed with true+scatter coincidences. Study with other factors also does not affect image in the expected way. So, we think that at some point our scatter correction method might not work properly but we are also wondering if the scatter correction in CASToR is also implemented correctly.

To check it we performed one more study - simulation of the cuboid phantom with very high statistics. We performed the reconstruction without TOF resolution for such cases:
- true coincidences
- true + scatter coincidences
- true + scatter coincidences with scatter correction performed by CASToR internally.
The resulted reconstructed images are shown in cuboid.png:
TOP LEFT: true coincidences
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true + scatter coincidences with scatter correction performed by CASToR internally
BOTTOM RIGHT: profiles along the orange line from previous plots
The number of scatter coincidences registered per LOR is around 10-30 coincidences.
What we found is that the images reconstructed with true+scatter and true+scatter with incorporated scatter correction look almost the same in terms of image quality and voxel activity units, which is surprising.

In conclusion, we are wondering if you have ever found similar things and could help us understand why the scatter correction does not work (especially for the last case). The second thing is that maybe there is a feature of the CASToR which we do not consider.

We will appreciate your help.

All the best,
Jakub

[sstute]

Hello Jakub,

Thanks for the detailed email !

Do you take the scan duration into account in the datafile header ?

You say that you estimated the scatter rate for each TOF bin independently to be 1-2. Is it "1 or 2" (which would be insanely high) or "e-2" ?

How did you compute this estimated scatter rate ?

And how did you translate this information into your list-mode file ?

Best
Simon

[jbaran1]

Hi Simon,

many thanks for the quick response. The answers are as follows:

Hello Jakub,

Thanks for the detailed email !

Do you take the scan duration into account in the datafile header ?

Yes, I added the information about the scan duration in the header. On the other hand in the reconstruction I used -ignore-corr fdur option.

You say that you estimated the scatter rate for each TOF bin
independently to be 1-2. Is it "1 or 2" (which would be insanely high)
or "e-2" ?

I was not explicit here. The scatter rate passed to the castor lm is ranging from 0 to 2. However most values are in the range from 1e-3 to 1e-1.

How did you compute this estimated scatter rate ?

We made basic estimation based on Monte Carlo simulations - NEMA IEC study. For the cubic study we use -sc flag during the conversion of the ROOT data to CASToR. For this the scatter rate is at the level of about 1e-4.

And how did you translate this information into your list-mode file ?

We adapt the CASToR code to read the scatter correction rate factor from the matrix provided by ourselves.

Best
Simon

All the best,
Jakub

[sstute]

Hi,

OK thanks for the explanations.

Why ignoring the frame duration in the reconstruction ?

As the scatter correction is supplied as a rate, you really need the frame duration to rescale the scatters to the prompt coincidences.

Could you try to reconstruct with the frame duration taken into account ?

How did you translate the estimation of the scatter rates based on a sinogram to the list-mode data ? Did you take the binning size into account ? You need to divide your scatter rate for a TOF bin by the size of the TOF bin to get an "almost" continuous measurement.

Simon

[jbaran1]

Hi Simon,

thanks for your prompt reply. and sorry for the late response.

I did as you suggest. I remove the -ignore-corr fdur flag from the reconstruction. I attached the results for both, NEMA IEC phantom (scatter fraction calculated by ourselves) - NEMA_IEC_fdur.png - and the cuboid phantom (scatter correction calculated by CASToR) - cuboid_fdur.png.

For the NEMA IEC plot, all presented images are reconstructed with scatter correction implemented and tested by ourselves. The difference is in an applied -ignore-corr fdur flag and scaling factor of the scatter rates and the legend is as follows:
TOP LEFT: scalling factor = 20000 with -ignore-corr fdur
TOP RIGHT: scalling factor = 1 (without -ignore-corr fdur frame duration is incorporated)
BOTTOM LEFT: scalling factor = 20000 superimposed on scalling factor = 1 with -ignore-corr fdur (the images looks the same)
BOTTOM RIGHT: profiles along the orange line from previous plots

You can see that three of the reconstructed images looks very similar. The most significant difference is in the image on TOP LEFT so the one with scalling factor = 20000 with -ignore-corr fdur. But still the effect of the scatters is observable.

Regarding the way of scatter rate estimation we did in principle as follows. We calculate the scatter corrections for each TOF bin. Then we merge some LORs to calculate the let say downsampled scatter rate. Then we divided this number by the number of LORs and time. In principle you can find some similarity to the SSRB procedure.

For the cuboid study where the CASToR calculate the scatter rate factors internally, the results doesn't change. For this case I reconstruct images once again but without -ignore-corr fdur flag. Incorporation of the scatter rate does not give any improvement of the image quality. It seems that even with nonTOF data these is something wrong with scatter correction.

The legend for the plot (cuboid_fdu.png) is as follows:
TOP LEFT: true + scatter coincidences scatter corrected
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true coincidences
BOTTOM RIGHT: profiles along the orange line from previous plots.

All the best,
Jakub

[sstute]

Hi Jakub,

Happy new year to you and your project !

About real data, what is the scan duration ? Could you share the header of the datafile ?
I will try to modify as soon as possible the DatafileExplorer toolkit to display the total corrected counts in addition to the total number of counts, that that you will be able to see if everything is globally correct in the datafile.

About simulated data, CASToR does not perform any scatter correction, so there is something I do not understand.
I am not used to the Gate2CASToR toolkit, sorry.
Anyway, in the figure you attached, I cannot believe that bottom left is true coincidences only, there should be something wrong.

How many iterations/subsets are you using in the reconstructions ?

Best
Simon

[jbaran1]

Hi Simon,

Happy New Year for you also.

All datasets which I shown are from the simulation datasets.

Hi Jakub,

Happy new year to you and your project !

About real data, what is the scan duration ? Could you share the
header of the datafile ?

I assumed that you mean the data from NEMA IEC siulation. The header looks as follows:

Data filename: factor_20000_df.Cdf
Number of events: 298175691
Data mode: list-mode
Data type: PET
Start time (s): 0
Duration (s): 500
Scanner name: TB_JPET_2nd_gen
Calibration factor: 1
Isotope: unknown
TOF information flag: 1
TOF resolution (ps): 230
List TOF measurement range (ps): 3000
List TOF quantization bin size (ps): 100
Scatter correction flag: 1

I will try to modify as soon as possible the DatafileExplorer toolkit
to display the total corrected counts in addition to the total number
of counts, that that you will be able to see if everything is globally
correct in the datafile.

We developed such a toolkit (more like dumping the data to the external file). Could you please tell me what I should look like for?

About simulated data, CASToR does not perform any scatter correction,
so there is something I do not understand.
I am not used to the Gate2CASToR toolkit, sorry.
Anyway, in the figure you attached, I cannot believe that bottom left
is true coincidences only, there should be something wrong.

I am sorry I made a mistake in the description of the cuboid.png file. It should be as follows:

The legend for the plot (cuboid_fdu.png) is as follows:
TOP LEFT: true + scatter coincidences scatter corrected
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true coincidences
BOTTOM RIGHT: profiles along the orange line from previous plots.

How many iterations/subsets are you using in the reconstructions ?

We used 20 iterations. The results which I sent are for the 15th iteration.

[jbaran1]

Hi Simon,

I forget to add the proper legend to the cuboid.png file:

The legend for the plot (cuboid_fdur.png) is as follows:
TOP LEFT: true coincidences
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true + scatter coincidences scatter corrected
BOTTOM RIGHT: profiles along the orange line from previous plots.

All the best
Jakub

[sstute]

Hi Jakub,

I am confused as you give two different legends for the cuboid data in your two latest emails.

OK for being both simulated data, indeed, otherwise you could not perform your SSRB-like scatter estimation... Sorry.

Your frame duration is significant. As you take it into account when computing the scatter estimation (with your method), I do not understand why you initially forced it to be ignored in the reconstruction ? You should always take it into account.

The way you estimate the scatter seems to be in agreement with what is expected into CASToR.

For the cuboid data, I still do not understand how CASToR is including a scatter correction. It may be that the real scatter count of each LOR is written as a rate in the scatter correction term. If so, then, it is not in agreement with the model used in CASToR as the scatter rate is included in the forward model (in other words, it is not subtracted). Anyway, if it is done this way, then it should not be done this way and you should not take this as a reference or so. You should definitely prefer your SSRB-like method of smoothing the scatter estimate.

About your cdh file, everything looks fine.

About what to look for in the datafile, I would like to see the total number of prompts VS total number of scatters, based on the CASToR datafile, and also the total number of corrected counts. You will be able to check if the scatter fraction is in agreement with the data coming out from GATE.

Best
Simon

[jbaran1]

Hi Simon.

Hi Jakub,

I am confused as you give two different legends for the cuboid data in
your two latest emails.

Sorry for the mistake. I double-ckecked and the legend for the plot (cuboid_fdur.png) is as follows:
TOP LEFT: true coincidences
TOP RIGHT: true + scatter coincidences
BOTTOM LEFT: true + scatter coincidences scatter corrected
BOTTOM RIGHT: profiles along the orange line from previous plots.

OK for being both simulated data, indeed, otherwise you could not
perform your SSRB-like scatter estimation... Sorry.

Your frame duration is significant. As you take it into account when
computing the scatter estimation (with your method), I do not
understand why you initially forced it to be ignored in the
reconstruction ? You should always take it into account.

I ignore it as it was from the previous study and I forget to turn it off The last results are with the scan duration information incorporated into reconstruction as you suggested.

The way you estimate the scatter seems to be in agreement with what is
expected into CASToR.

For the cuboid data, I still do not understand how CASToR is including
a scatter correction. It may be that the real scatter count of each
LOR is written as a rate in the scatter correction term. If so, then,
it is not in agreement with the model used in CASToR as the scatter
rate is included in the forward model (in other words, it is not
subtracted).

I think it is done like you said. From my understanding, CASToR calculate the number of scatters and trues and then put the information about the scatter rate term into list-mode as the number of coincidences divided by the scan time.

Anyway, if it is done this way, then it should not be
done this way and you should not take this as a reference or so.
You should definitely prefer your SSRB-like method of smoothing the
scatter estimate.

About your cdh file, everything looks fine.

About what to look for in the datafile, I would like to see the total
number of prompts VS total number of scatters, based on the CASToR
datafile, and also the total number of corrected counts. You will be
able to check if the scatter fraction is in agreement with the data
coming out from GATE.

We are using the list-mode reconstruction so if it is satisfactory for you I can give you the numbers of forward projection (vOptimizer::DataStep1ForwardProjectModel( ... ) method) from the line: m2p_forwardValues[a_th][b] += ap_Event->GetAdditiveCorrections(b) * mp_ImageDimensionsAndQuantification->GetFrameDurationInSec(a_bed, a_timeFrame)

In the attached files I provide a sample of the dumped data for consecutive events. The values in columns correspond to:
1) m2p_forwardValues[a_th][b] (before adding the
  > additiveCorrections)
2) ap_Event->GetAdditiveCorrections(b)
3) m2p_forwardValues[a_th][b] (after adding the additiveCorrections) -

The simulation time is 500 seconds. I dumped the data from 4 iterations for first 957 events. The two files correspond to the two cases. In the first case (castor_dev_factor1.txt) we used the values as described above. In the second case (castor_dev_factor20000.txt) we added an extra factor =20000 that multiply the scattering estimates. The factor was obtained just by trial and guess to match the scale between the additive correction part (scatter) and forward projection of the image estimate.

Besides I also attached two short presentations of one of my colleague where he summarizes quantitatively all the staff we discuss. There is also comparison between our data (Scatter_factor_CASTOR_part_V) and the CAStoR benchmark (Scatter_factor_CASTOR_part_VI).

Concerning the differences I also have another question. We have analyzed the benchmark data from the castor website (GE SIGNA PET/MR scanner). We found that the multiplicative factor - m_multiplicativeCorrection (oProjectionLine::ForwardProject(FLTNB* ap_image) fromm line value /= m_multiplicativeCorrection;) is 1 in our case and have various values in castor benchmark data. It seems that it is a crucial difference in those cases (as seen in the presentations). Could you elaborate what is its definition and how is it calculated?

The last thing which I would like to ask is a calibration factor given in list-mode header from your benchmark data. Its value is Calibration factor: 1.45762e+08. What is the definition of this value? CASToR documentation states only: "The calibration factor associated to the current scanner settings (default: 1.)"

All the best,
Jakub

castor_dev_factor2000.txt (107 KB)

castor_dev_factor1.txt (143 KB)

Scatter_factor_CASTOR_part_VI.pdf (1010 KB)

Scatter_factor_CASTOR_part_V.pdf (740 KB)

[jbaran1]

Hi Simon,

I would like to kindly ask if you found some time to investigate the issue we discussed here.

All the best,
Jakub

[sstute]

Hi Jakub,

Sorry I forgot, I have just had a look.

About the "CASToR scatter correction", as it is indeed what you say, you can totally forget it and looking at the results would be misleading. I suggest that you only compare "trues only" to "including scatters and using your correction". Did you try to compute the scatter correction by yourself for the cuboid study ?

Is the GATEtoCASToR tool able to provide a histogram datafile ? If not, I really do not understand how it may be able to provide a the number of scatter counts for a given LOR. If yes, then I suggest that you first try a non-TOF histogram reconstruction to validate your scatter correction method. Then with histogram data, you would be able to dump the EventValue beside the additive corrections. If you want to stick with listmode data, I suggest that you try to include the scatters and your scatter correction without TOF first. Then once it works, you add TOF.

The benchmark is adapted from a really FDG scan with all corrections computed from the manufacturer's software.
The multiplicativeFactor as dump from the code is a combination of all multiplicate correction terms: the attenuation correction factor, the normalization factor, the decay correction factor, the frame duration and the calibration factor that you mention. So I am surprised it is 1 in your case. If you take the frame duration into account, it should be 1/500.
About the calibration factor, it is the standard calibration factor embedded into any PET scanner enabling absolute quantification. Without this factor, you would reconstruct the "seen activity" but not the "actual activity included within the field of view of the scanner".

All the best
Simon

[jbaran1]

Hi Simon,

many thanks for you reply.

Hi Jakub,

Sorry I forgot, I have just had a look.

About the "CASToR scatter correction", as it is indeed what you say,
you can totally forget it and looking at the results would be
misleading. I suggest that you only compare "trues only" to "including
scatters and using your correction". Did you try to compute the
scatter correction by yourself for the cuboid study ?

No, we did not. So far we did some work only for the NEMA IEC.

Is the GATEtoCASToR tool able to provide a histogram datafile ? If
not, I really do not understand how it may be able to provide a the
number of scatter counts for a given LOR. If yes, then I suggest that
you first try a non-TOF histogram reconstruction to validate your
scatter correction method. Then with histogram data, you would be able
to dump the EventValue beside the additive corrections. If you want to
stick with listmode data, I suggest that you try to include the
scatters and your scatter correction without TOF first. Then once it
works, you add TOF.

It provides the histogram data format. We will try with that so.

The benchmark is adapted from a really FDG scan with all corrections
computed from the manufacturer's software.
The multiplicativeFactor as dump from the code is a combination of all
multiplicate correction terms: the attenuation correction factor, the
normalization factor, the decay correction factor, the frame duration
and the calibration factor that you mention. So I am surprised it is 1
in your case. If you take the frame duration into account, it should
be 1/500.

In fact we do not include the correction factors in the list mode. We work with simulations so far so the normalization is set to 1 so far. Regarding the sensitivity and attenuation we merge both to the sensitivity map using CASToR projector SENS (it has been shared with us some time ago by Thibaut for testing). Do you think we should include some other correction factors along with that map?

About the calibration factor, it is the standard calibration factor
embedded into any PET scanner enabling absolute quantification.
Without this factor, you would reconstruct the "seen activity" but not
the "actual activity included within the field of view of the
scanner".

Many thanks for your explanation. We will try to reconstruct without ToF first.

Best regards,
Jakub

[sstute]

Hi Jakub,

When doing listmode reconstructions, even if do not include attenuation and normalization factors in the list-mode, you should see the 1/500 in the multiplicativeCorrection term.

Then, if you consider correcting for attenuation and normalization, you have to include these information in the list-mode, otherwise, the scatter correction cannot work ! The idea of including attenuation and normalization only in the sensitivity image assumes that (i) you are using OSEM and (ii) you are not correcting for background events (randoms and scatters). It is a strong simplification. So this seems to be the explanation why your scatter correction does not work.

For your list-mode reconstruction, try to include attenuation correction in the list-mode file and relaunch the reconstruction including the computation of the sensitivity image considering the attenuation.

Do not expect to accurately correct for normalization with the SENS stuff. It is a very very rough approximation of the normalization correction that can only work without including any correction of background events, and I have no idea how it could work if you want to get absolute quantification.

For computation of the sensitivity image, I suggest that you build a "normalization datafile" as described in the documentation because it will allow you to clearly define the set of possible LORs in your system.

Then, if you want to go further in the reconstruction process for your JPET data, I strongly suggest that you take some time to accurately compute normalization factors of your system. These factors are way off being 1 even for simulations. Also the method that you will build can be used later for the real JPET system.

Best
Simon

[tmerlin]

Hi Jakub,

Just to clarify, the 'sens' optimizer uses castor-recon to derive a sensitivity image from a CASToR datafile containing an "uniform" acquisition (the setup of this acquisition being dependent on the geometry of the system). Its main purpose is to be used as a workaround to get a sensitivity image with an approximation of the normalization correction factors for simulation data. I second Simon's comments regarding the issues with background events, as using such sensitivity image with castor-recon list-mode reconstruction bypasses the computation of the quantification factors (which must be calculated by the user) and create discrepancies with the actual datafile containing scatter and random data since this datafile must contains the same normalization correction factors as those which were used during the generation of the sensitivity image.

Regarding the histogram format, castor-GATERootToCastor can convert data to histogram format, but without TOF correction. TOF histogram would require an additional step to sort the events in the different TOF bins depending on the timing resolution of the system, its geometry, and the user-defined binning. Anyway I agree with Simon that the best is to start with non TOF histogram data to assess your scatter correction method.

Best,
Thibaut

[Xin_Zhao]

Dear Thibaut,
Can you recommend some references about how to compute scatter & random rate for each LOR so that I can insert them into the list-mode ".Cdf " files to apply background correction? Any help would be appreciated！
Best，
Xin

[Zhengkun]

Hi,

Did you solve the issue you encountered in 2021?

Best wishes,
Zhengkun