Issue when using a double gaussian operator

[akippen]

Hi there,
We are trying to perform a SPECT reconstruction using our own PSF operator. We have fit the PSF response to a double gaussian mixture in order to better account for the tails. This fitting has given results for FWHM that seem reasonable, yet when we try to incorporate this new operator into the reconstruction the image is very poor.

We did the same method but using a single gaussian operator and obtained satisfactory results.
I was just wondering if we are doing something obviously wrong when defining the double gaussian operator, otherwise the issue may lie with our data or fitting.

See the image below for a comparison between the reconstructions.

When fitting the double gaussian we get weights for each, that sum to 1. The amplitude used when defining the operator is equal to

A = w / [sigma * sqrt(2(pi)]

This amplitude appears to vary as an exponential decay with distance, hence the form of ‘amplitude_fn’ in our code.
I have given our code used to define the double gaussian operator below. Any ideas or suggestions are greatly appreciated,
Thank you

amplitude_fn = lambda d, bs: (bs[0] * torch.exp(-bs[1] * d)) + bs[2]
sigma_fn = lambda d, bs: (bs[0]*d + bs[1])

amp_1_params = torch.tensor([1.497, 0.074, 0.414], device=device, dtype=torch.float32) # Input values from fitting
amp_2_params = torch.tensor([0.066, 0.287, 0.005], device=device, dtype=torch.float32) # Input values from fitting
sigma_1_params = torch.tensor(([0.015, 0.177]), device=device, dtype=torch.float32) # Input m and c from fitting
sigma_2_params = torch.tensor(([0.297, 0.767]), device=device, dtype=torch.float32) # Input m and c from fitting

operator_1 = GaussianOperator(
    amplitude_fn,
    sigma_fn,
    amp_1_params,
    sigma_1_params,
)

operator_2 = GaussianOperator(
    amplitude_fn,
    sigma_fn,
    amp_2_params,
    sigma_2_params,
)

double_operator = operator_1 + operator_2

double_operator.save(file_path)

[Janton]

For completeness, we both working on the same project.

We are getting familiar with the SPECTPSFToolbox.

We would like to model a double mixture of Gaussians resulting from measurements of Line sources. To account better for the tails of the Gaussian. The fittings look good with experimental data and we have normalised our distributions.
In tutorial 7 you seem to be able to load a 1D operator. However, we have not been able to create 1D operators (psf_operator_1D.pkl) so we implemented a 2D.

Following your tutorials, the steps to do a 1D operator seem to be in tutorial 5 but we could not find them (there is a cross-reference in tutorial 7 to 5). Ideally, we would like to implement 1D kernels into operators as the first step.

Any help will be highly appreciated
Best wishes

[lukepolson]

Hi @Janton @akippen . I’ve defined the operator myself and have been playing around with it. Based on the relative amplitudes and sigmas between the two Gaussian operators, it seems that at a distance of 25cm, the contribution from the “tail” Gaussian is approximately 2x that of the “geometric component” Gaussian.

Nx0 = 255
dx0 = 0.48
x = y = torch.arange(-(Nx0-1)/2, (Nx0+1)/2, 1).to(device) * dx0
xv, yv = torch.meshgrid(x, y, indexing='xy')
distances = torch.tensor([1,25,40,55.]).to(device)
input = torch.zeros((distances.shape[0],Nx0,Nx0)).to(device)
input[:,127,127] = 1
output1 = operator_1(input, xv, yv, distances)
output2 = operator_2(input, xv, yv, distances)
# Print contributions from each at all 4 distances
print(output1.sum(dim=(1,2)))
print(output2.sum(dim=(1,2)))

In the context of SPECT reconstruction, this would mean that many counts from far from the center contribute to the PSF (which I think is invalid for things like Lu-177), and thus you see those “blank” regions near the hot volumes in your double Gaussian reconstruction (these counts get assigned to the central hot region). Are your fit yields better MSE than a single Gaussian? I suspect this may be the issue.

[Janton]

Thanks a lot

It makes sense. yes the fitting to LSF was spot on. Perhaps use for modelling a more conservative approach for that double Gaussian (attempt to recover less).

Regarding the 1D operators. Is there a way to save 1D Operators? Something is mentioned in the tutorial 7. You opened the psf_operator_1D.pkl from tutorial 5, but we could not find how you defined and saved this operator in that tutorial.

Thanks