parallelproj pre-1.9.0 truncates kernel for large TOF bin size

[KrisThielemans]

I have weird results when reconstructing TOF data (forward projected with parallelproj via STIR) in a toy scanner. I know parallelproj does not (yet) attempt to make the sum over all TOF bins equal to the non-TOF result, but seem to find that the TOF bin size is ignored. In the example below, I set a really high CTR (10 ps), which effectively should use the TOF bin size, but it doesn't.

exam_info = stir.ExamInfo(stir.ImagingModality(stir.ImagingModality.PT))
scanner=stir.Scanner.get_scanner_from_name("GE Discovery 690")
scanner.set_timing_resolution(10)
scanner.set_num_rings(2)
span = 3
max_ring_diff = 1
tof_mashing = 11
num_views=scanner.get_num_detectors_per_ring()//2;
proj_data_info=stir.ProjDataInfo.construct_proj_data_info(scanner,
                                                 span, max_ring_diff,
                                                 num_views, scanner.get_max_num_non_arccorrected_bins(),
                                                 False, tof_mashing);
# create projdata with 1 TOF bin non-zero
proj_data=stir.ProjDataInMemory(exam_info, proj_data_info)
proj_data.fill(0)
sinogram = proj_data.get_empty_sinogram(0,0,False,0)
sinogram[stir.make_IntCoordinate(0,0)]=1;
proj_data.set_sinogram(sinogram)

# image
radius = 130 # mm
zoom=1
size_x = int(2*radius / (scanner.get_default_bin_size()/zoom))
imPP = stir.FloatVoxelsOnCartesianGrid(exam_info, proj_data_info, zoom, 
                                           stir.FloatCartesianCoordinate3D(0,0,0), stir.IntCartesianCoordinate3D(-1,size_x,size_x));
# backprojection
PETproj=stir.ProjectorByBinPairUsingParallelproj()
PETproj.set_up(proj_data_info, emission);
PETproj.get_back_projector().back_project(imPP, proj_data)

pm=stir.ProjMatrixByBinUsingRayTracing()
pm.enable_cache(False)
PETprojRT=stir.ProjectorByBinPairUsingProjMatrixByBin(pm)
PETprojRT.set_up(proj_data_info, emission);
imRT = imPP.clone()
PETprojRT.get_back_projector().back_project(imRT, proj_data)

parallelproj result: ray-tracing matrix result:

proj_data_info.get_sampling_in_k(stir.Bin(0,0,0,0,0)) returns 146.75mm, which seems to fit with the ray-tracing result (image radius 130mm)

(Tested with parallelproj 1.7.3, CPU version)

[KrisThielemans]

This causes problems for normal reconstructions, as we're (by default) assuming that the non-TOF sensitivity is equal to the TOF sensitivity image, but that isn't the case with parallelproj. In the above extreme case, this seems to give serious artefacts.

[KrisThielemans]

@gschramm is this known behaviour? Are you using a Gaussian kernel, as opposed to its integral (erf)?

I'm not sure what the best way around this is for now. I thought to (very crudely) pass sqrt(CTR^2+(TOF_bin_size/2)^2) to parallelproj (after conversion to sigma). What do you think?

[gschramm]

Hm, that is indeed very strange. The TOF projectors are supposed to correctly integrate the Gaussian over the bins (resulting in the difference of two error functions) - see e.g. here.

Could it be that this is related to our "truncation strategy"? By default, we neglect "voxel - tof bin center" pairs that are too far apart. In the practical implementation, we calculate the "relevant" TOF bins for every voxel along an LOR (TOF bins where the contributions are supposed to be > 0). This calculation, unfortunately, neglects the bin width which is wrong if the bin width is >> TOF res (but not sure if this case happens a lot in practice).

The function computing the "relevant TOF bins" is here.

Can you try with a bigger value for n_sigmas in the call of the projector? https://github.com/gschramm/parallelproj/blob/d5a4493e7ece6a7982f721a92a24d6937e5aa740/c/src/joseph3d_back_tof_sino.c#L25

[gschramm]

I guess the "effective" sigma, that is used to calculate the "relevant TOF bins" should be: $$\sqrt{\sigma\text{gauss}^2 + \sigma\text{tofbin rectangle}^2} $$ with $$\sigma_\text{tofbin rectangle} = \text{tofbin width} / \sqrt{12}$$

So far, I always assumed $\sigma\text{tofbin rectangle} << \sigma\text{gauss}$

[KrisThielemans]

If you're using erf, then the sum over all TOF bins should be non-TOF (except the edges). so maybe it is the truncation strategy indeed. Is it possible to update that to take bin width into account? Then there should be no need for effective sigma (we don't need it in the RT matrix).

We regularly run with larger TOF-mashing to save time and memory...

[KrisThielemans]

I checked with CTR 200ps and width 13.3ps, then everything is very close. with TOF mashing 5 (i.e width 66ps), the truncation causes problems:

[gschramm]

Hi Kris,

I will fix that on the C/CUDA level of parallelproj (using the effective sigma for the truncation). Currently running some more tests. Using the "effective sigma" + n_sigmas=3 will remove the "truncation" bug.

With "effective sigma" + n_sigmas=3 and a point source, (sum(tof) - nontof)/nontof < 0.004 for large/small tof bins and for different point source positions.

If that acceptable? Right now, I am struggling to find an efficient and general solution that would give (sum(tof) - nontof) = 0.

[KrisThielemans]

sure. there will always be numerical differences. thanks!

I probably should have created this issue on parallelproj repo, but I thought it was a STIR bug...

[gschramm]

https://github.com/gschramm/parallelproj/issues/69

[gschramm]

Hi Kris,

I submitted this PR to fix the kernel truncation issue (merges the tof branch). Can you try re-running your tests with the updated parallelproj code?

According to my tests, the residual underestimation of the TOF sum, due to kernel truncation when using 3 sigmas, is between:

• 0.003 -> if tof bin with << sigma tof
• 0.006 -> max value when tof bin with = 3 * sigma tof
• 0.002 -> if tof bin with << 5*sigma tof
• close to 0, if tof bin width >> sigma tof

The negative bias in the case of tof bin with << sigma tof) can be calculated analytically (erf(n_sigmas/sqrt(2)). Should I correct for this "global" bias? I think it would reduce bias for most practical cases (although it will introduce slightly positive bias in the case tof bin width >> sigma tof)

[gschramm]

just added the "global correction". this means that for voxels that line up with the center of the TOF bins, the sum over TOF should be identical to non-TOF now - even when using a kernel truncation with n_sigmas = 3. (also for n_sigmas < 3, but then the kernel shape deviates more from a true gaussian convolved with a rectangle)

[KrisThielemans]

Fixed in parallelproj 1.9.0. Thanks @gschramm !

Maybe we should require that version...