odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
375 stars 105 forks source link

Dishing artifact and value shift with FanBeamGeometry #1630

Open Cs-Olasz opened 1 year ago

Cs-Olasz commented 1 year ago

Dear Developers and users of ODL,

I started with the example code of https://github.com/odlgroup/odl/blob/master/examples/tomo/filtered_backprojection_cone_2d.py. I can run it without error now and the result that I get is also as it should be. But after some changes in the code the result I get has a dishing artifact and a shift in the intensity values. Do you have any idea why I getting this values?

About the changes I made: basically, I just want to run the code with a loaded image. And as my test image has a size of 512x512 I made adjustment regarding the geometry.

Note, that if I use the Shepp-logan phantom instead of mine, then it seems working fine.

Thank you in advance for the answers!

Regards, Csaba Olasz

Source code:

import numpy as np import odl import tifffile

recon_size = [512, 512]

reco_space = odl.uniform_discr([-recon_size[0] / 2, -recon_size[1] / 2], [recon_size[0] / 2, recon_size[1] / 2], recon_size, dtype='float32')

angle_partition = odl.uniform_partition(0, 2 * np.pi, 360) detector_partition = odl.uniform_partition(-500, 500, 1000, 1)

geometry = odl.tomo.FanBeamGeometry( angle_partition, detector_partition, src_radius=500, det_radius=400)

ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')

fbp = odl.tomo.fbp_op(ray_trafo, filter_type='Ram-Lak', frequency_scaling=0.8)

""" phantom = odl.phantom.shepp_logan(reco_space, modified=True)""" phantom = tifffile.imread('cylinder_slice.tif') phantom = reco_space.element(phantom)

proj_data = ray_trafo(phantom)

fbp_reconstruction = fbp(proj_data)

tifffile.imwrite('rec_cylinder.tif', fbp_reconstruction)

phantom.show(title='Phantom') proj_data.show(title='Projection Data (Sinogram)') fbp_reconstruction.show(title='Filtered Back-projection') (phantom - fbp_reconstruction).show(title='Error', force_show=True)

cylinder_slice.zip