How to use the iterative reconstructors？

[gaopinghai]

I can use the fbp reconstructor properly. The code is like below:

import dival
from dival.reconstructors.odl_reconstructors import FBPReconstructor

num_angels = 50
dataset = dival.get_standard_dataset('lodopab', num_angles=num_angels)
fbp_reconstructor = FBPReconstructor(dataset.get_ray_trafo(), hyper_params={
    "filter_type": "Hann",
    "frequency_scaling": 0.641025641025641
})

imorg = <image from lodopab dataset>  # shape = (362, 362)
ob = <observations from lodopab>  # shape = (1000, 513)

num_thetas = oborg.shape[0]
thetas = np.linspace(0, np.pi, num_thetas+1)[:-1]

inds = range(0, num_thetas, num_thetas//num_angels)  # select num_angles of observation data
nthetas = thetas[inds]
nob = oborg[inds]

img = fbp_reconstructor.reconstruct(nob)

The code above works fine. When the reconstructor change from FBPReconstructor to ISTAReconstructor, how should I write my code? I do not want the code in the dival/examples/ct_example.py which is too much packaged. Thanks !

[gaopinghai]

I think this issue is related to #35 ，where only FBP reconstructor is used. So the Iterative reconstructor should also be included.

[jleuschn]

Concerning the usage of the ISTAReconstructor, the constructor receives different parameters, see e.g. the docs. It expects an initial image x0 and a number of iterations niter, and the hyper parameters of the FBP do not apply. The construction call could thus look like this:

ista_reconstructor = ISTAReconstructor(dataset.get_ray_trafo(),
                                      dataset.space[1].zero(),
                                      1000)

By the way, I found that the iterative reconstructors wrapping ODL ignore the hyper parameter iterations inherited from IterativeReconstructor, and just filed an issue #38 about this. For now, one can just pass niter as an argument to the constructor as in the code above.

As for the angle subsampling I just would like to point out that the angles created with thetas = np.linspace(0, np.pi, num_thetas+1)[:-1] are not the exact angles used in the LoDoPaB data. Those can be found in dataset.geometry.angles.

[gaopinghai]

Thanks for your help again, sir!

Firstly, for the angles used in LoDoPab data, is not exactly as np.linspace(0, np.pi, 1001)[:-1] which I did not notice! But the only difference is the starting angles. One starts from zero and the other starts from about 0.0015. The interval is the same.

Secondly, The reconstruction results of the iterative methods seems very poor, which is not reasonable. What's the problem? The code is as below:

import dival
from dival.reconstructors.odl_reconstructors import (FBPReconstructor,
                                                     GaussNewtonReconstructor,
                                                     LandweberReconstructor,
                                                     ISTAReconstructor,
                                                     DouglasRachfordReconstructor)
num_angels = 1000
dataset_odl = dival.get_standard_dataset('lodopab', num_angles=num_angels)

inds = range(0, 1000, 1000//num_angels)
ob = oborg_lpd[inds]
angles_odl = dataset_odl.geometry.angles

ray_trafo = dataset_odl.get_ray_trafo()

fbp_reconstructor = FBPReconstructor(dataset_odl.get_ray_trafo(), hyper_params={
    "filter_type": "Hann",
    "frequency_scaling": 0.641025641025641
})
img_fbp = fbp_reconstructor.reconstruct(ob)

gn_reconstructor = GaussNewtonReconstructor(ray_trafo, dataset_odl.space[1].zero(), 1000)
img_gn = gn_reconstructor.reconstruct(ob)

ldw_reconstructor = LandweberReconstructor(ray_trafo, dataset_odl.space[1].zero(), 1000)
img_ldw = ldw_reconstructor.reconstruct(ob)

ista_reconstructor = ISTAReconstructor(ray_trafo, dataset_odl.space[1].zero(), 1000)
img_ista = ista_reconstructor.reconstruct(ob)

dr_reconstructor = DouglasRachfordReconstructor(ray_trafo, dataset_odl.space[1].zero(), 1000)
img_dr = dr_reconstructor.reconstruct(ob)

And the result is:

[gaopinghai]

I found the standard reconstructed results from here, so I think there is not much necessity for me to reproduce it.

Sir，I still need your help on how to write the correct code to reconstruct the slices from projection data. Please help! @jleuschn

[gaopinghai]

I found the example code here