LoDoPaB Reconstruction Intensity Scaling

[tibuch]

Hi,

Thank you for providing this awesome resource!

I am looking at the LoDoPaB dataset and have a question regarding the intensity scaling of the reconstructions. I am comparing the ground-truth image to the noisy reconstruction and they seem to have completely different intensity ranges.

I run the following code to obtain the described results:

import dival
import numpy as np

lodo = dival.get_standard_dataset('lodopab')
rec = dival.get_reference_reconstructor('fbp', 'lodopab')

x, y = lodo.get_sample(0)
rec_x = rec.reconstruct(x)

def stats(name, d):
    print('{}: min={}; mean={}; max={}'.format(name, np.min(d), np.mean(d), np.max(d)))

stats('y',y)
stats('reconstruction x', rec_x)

Which produces this output:

y: min=0.0; mean=0.13288576900959015; max=0.48916396498680115
reconstruction x: min=-3.407857112822099e-12; mean=5.093645644160816e-11; max=1.6748728792759238e-10

The reconstruction intensities are much smaller. I am wondering where this difference comes from and how I could obtain reconstructions which are in the same intensity range as the ground truth.

Note: Visually the reconstructions look sound.

Left reconstruction and right GT.

[jleuschn]

Hi,

glad you like the library.

The scaling issue is almost certainly due to incompatible versions of the astra and odl libraries. Using the newest versions of both the stats are the same for me. As a quick check you could specify impl='skimage' to get_reference_reconstructor, which then switches to the slow scikit-image backend.

Since odl checks the version of astra, i would think installing the newest git-master version of it would already fix the issue:
pip install https://github.com/odlgroup/odl/archive/master.zip

I also use the latest dev version of the astra-toolbox, which can be installed in conda with:
conda install astra-toolbox -c astra-toolbox/label/dev.

Hope you can solve it (without breaking other version dependencies ;))!

[tibuch]

Awesome, this resolved the issue. Thank you.

Now I have a followup question:
I have a Sinogram y and use the corresponding reconstructor to obtain reconstruction rec_y. Now if I use the corresponding ray-trafo to go from rec_y to projected_rec_y. projected_rec_y compares to y as expected, identical up to noise.

Now I wanted to understand what is really happening in the ray-trafo. So I rotated the image and summed along the projection axis to obtain my_projection_rec_y. my_projection_rec_y looks comparable to y up to a huge scaling factor. To get the result into the same range as the y I have to divide it by ~1200.
It is not clear to me where this factor comes from.

[jleuschn]

Sorry, i don't really know, never implemented a ray_trafo from scratch myself. One could check the code in odl.tomo for a specific case and backend in more detail.
Two potential (maybe related) sources of scaling differences could be:

• ODL takes into account the physical dimensions (for LoDoPaB images min_pt=[-0.13,-0.13], max_pt=[0.13, 0.13], so 26x26 cm).
• Spaces in ODL can have a weighting (lodo.space[0].weighting == 2.2517535e-06, lodo.space[1].weighting == 5.1585727e-07)

I do not know about the details of ODL here, but wanted to stay close to its default behaviour for LoDoPaB.

[tibuch]

Okay, thanks for the pointers.

Is the upscaling part of lodo.get_ray_trafo() or is it just the projection? Also is there some pre-/post-processing going on with regard to MU_MAX?

Thank you for your help.

[jleuschn]

The ray trafo object return by get_ray_trafo is the (noiseless) operator mapping from the images to the observations, there the scaling matches. The factor 1./MU_MAX was applied to both the images and the observations, so since it is linear it does not make a difference to the ray trafo. You only would need to consider it for the pre-log model.

[tibuch]

Okay. Just for clarification: The operator mapping from images to observations is the combination of upscaling to 1000x1000 pixel followed by the projection to 513 pixel.

Thank you for your detailed answers.

[jleuschn]

Yes, that's right, for simulation a bilinear upscaling to 1000x1000 was included before applying a ray trafo from 1000x1000 to 1000x513 (and adding noise). Thanks for adding this. The operator returned by get_ray_trafo is on the other hand directly mapping from 362x362 to 1000x513 and is approximately equivalent.

[tibuch]

Figured it out. There is a special scaling factor used in the odl-projection.

lodo = dival.get_standard_dataset('lodopab')
scale = lodo.space[0].cell_sides[1]

So now np.sum(reconstruction, axis=1)*scale results in the correct scaling of the projection. This helped me figure out the scaling factor.

[jleuschn]

Great, thanks for sharing!