How to use part of the observations to reconstruct？

[gaopinghai]

Hi，thanks for this great Dataset and python tool. It makes people like me who is new to CT easier to explore this area.

Now I managed to reconstruct the image using fbp algorithm with the observation data. The code is like below:

fbp_reconstructor = dival.get_reference_reconstructor('fbp', 'lodopab')
fbp_img = fbp_reconstructor.reconstruct(obdata)  # obdata.shape: (1000, 513)

I have two questions here: 1. How to reconstruct the image from part of the observation datas? That means the origne observation data is 1000 angles. When the angle numbers becomes 1/2, which is 500, the observation data while be:

inds = np.round(np.linspace(0, 999, 500)).astype(np.int32)
obdata_half = obdata[inds]

How to get the reconstruction result from this observation data? E.g. how to run img = fbp_reconstructor.reconstruct(obdata_half) correcttly?

2. Which iterative algorithm should I choose? There are so many iterative algorithms? Which one is the best or the newest, and I should test it? Please recommend for me, thanks!

That's my problem, thanks for help!

[MBaltz]

Hello @PingHGao Anserwing the first question: you need to choose the number os angles when crating the dataset. Like the code below:

IMPL='astra_cuda'
dataset_500 = get_standard_dataset('lodopab', impl=IMPL, num_angles=500)
sample_pairs = dataset_500.get_data_pairs_per_index('test', [0, 15, 83])
fbp_reconstructor = get_reference_reconstructor('fbp', 'lodopab', impl=IMPL)
for obs, gt in sample_pairs:
    rec_obs = fbp_reconstructor.reconstruct(obs)
    # ...

This number of angles (_numangles) needs to be divisor of 1000. Documentation of _get_standarddataset: link

[jleuschn]

Hi,

thanks a lot for commenting @MBaltz , specifying num_angles to get_standard_dataset as you propose is also my preferred way to subsample, this way the ray_trafo of the dataset is also adapted accordingly. However, the reference reconstructor is still defined for the 'lodopab' dataset, so it will not work with the subsampled observations. If you would like to use the same hyperparameters for an FBP working on 500 angles, you can just construct another FBPReconstructor. Here's a MWE based on @MBaltz 's snippet:

from dival import get_standard_dataset, get_reference_reconstructor
from dival.reconstructors.odl_reconstructors import FBPReconstructor

IMPL='astra_cuda'
dataset_500 = get_standard_dataset('lodopab', impl=IMPL, num_angles=500)
sample_pairs = dataset_500.get_data_pairs_per_index('test', [0, 15, 83])
fbp_reconstructor = get_reference_reconstructor('fbp', 'lodopab', impl=IMPL)
fbp_reconstructor_500 = FBPReconstructor(dataset_500.get_ray_trafo(impl=IMPL),
        hyper_params=fbp_reconstructor.hyper_params)
for obs, gt in sample_pairs:
    rec_obs = fbp_reconstructor_500.reconstruct(obs)
    # ...

Regarding you second question, i would recommend TV. It usually yields quite good results when using suitable hyperparameters. Another option would be CGLS, and for 3D it seems to me that Kaczmarz has the benefit of saving some memory (in each iteration considering only the projection from one angle).

Best regards, Johannes

[MBaltz]

Thanks for correcting my example. Have a nice day.

[gaopinghai]

@jleuschn sir, I want to know how the filter_type and frequency_scaling params is chose and set. And as far as I know, I think fbp_reconstructor.hyper_params which is "filter_type": "Hann", "frequency_scaling": 0.641025641025641 is dataset specified. I want normal params, which are not dataset specified. What the params should be? Thanks!

[jleuschn]

Those were set based on a hyperparameter search on validation images of the dataset. The optimal hyperparameters will depend on the setup and data.

From my point of view a Hann filter and frequency scaling 1.0 are a good starting point. If the noise introduces too pronounced artifacts, you can reduce the frequency scaling. This has a smoothing effect, and in fact some high-frequency information gets lost, but this is what you typically want to deal with noise. Of course, also a different kind of filter could be better in your case.

Best regards, Johannes