Question about phantom data

[AceCoooool]

The value of ellipses is "unsimilar" with shepp logan? which is difficult to post-process methods.
For example, the demo code as follow:

import odl
from adler.odl.phantom import random_phantom
from skimage.measure import compare_ssim as ssim
from skimage.measure import compare_psnr as psnr

size = 128
space = odl.uniform_discr([-64, -64], [64, 64], [size, size], dtype='float32')

geometry = odl.tomo.parallel_beam_geometry(space, num_angles=180)
operator = odl.tomo.RayTransform(space, geometry)

phantom = random_phantom(space)
phantom_data = operator(phantom)
# phantom_data += odl.phantom.white_noise(operator.range) * np.mean(np.abs(phantom_data)) * 0.05

shepp = odl.phantom.shepp_logan(space, modified=True)
shepp_data = operator(shepp)
# shepp_data += odl.phantom.white_noise(operator.range) * np.mean(np.abs(shepp_data)) * 0.05

with odl.util.Timer('runtime of FBP reconstruction'):
    phantom_recon = odl.tomo.fbp_op(operator, filter_type='Ram-Lak')(phantom_data)
    shepp_recon = odl.tomo.fbp_op(operator, filter_type='Ram-Lak')(shepp_data)

print('shepp ssim = {}'.format(ssim(shepp.asarray(), shepp_recon.asarray())))
print('shepp psnr = {}'.format(psnr(shepp.asarray(), shepp_recon.asarray(), data_range=1)))
print('shepp psnr(dynamic) = {}'.format(psnr(shepp.asarray(), shepp_recon.asarray())))

print('phantom ssim = {}'.format(ssim(phantom.asarray(), phantom_recon.asarray())))
print('phantom psnr = {}'.format(psnr(phantom.asarray(), phantom_recon.asarray(), data_range=1)))
print('phantom psnr(dynamic) = {}'.format(psnr(phantom.asarray(), phantom_recon.asarray())))

shepp, phantom = shepp.asarray(), phantom.asarray()
shepp_proj, phantom_proj = shepp_data.asarray(), phantom_data.asarray()
shepp_recon, phantom_recon = shepp_recon.asarray(), phantom_recon.asarray()

print('shepp: min ', shepp.min(), ' max ', shepp.max())
print('shepp_proj: min ', shepp_proj.min(), ' max ', shepp_proj.max())
print('phantom: min ', phantom.min(), ' max ', phantom.max())
print('phantom_proj: min ', phantom_proj.min(), ' max ', phantom_proj.max())

And the results as follow:

 runtime of FBP reconstruction :      0.026 
shepp ssim = 0.9221367003246534
shepp psnr = 25.886712891172724
shepp psnr(dynamic) = 31.90731280445235
phantom ssim = 0.9731283890622408
phantom psnr = 32.4274721004839
phantom psnr(dynamic) = 38.44807201376352
shepp: min  -1.49012e-08  max  1.0
shepp_proj: min  0.0  max  33.2351
phantom: min  -0.911248  max  0.81557
phantom_proj: min  -45.1083  max  20.4536

The value of phantom is nearly [0, 1]; however the value of phantom has negative values (Althought it's random)

[adler-j]

Indeed the Shepp-Logan case is somewhat of an outlier (which is partially why it was chosen), but it's still statistically possible given our prior.

Could you please elaborate on what the question is here?

[AceCoooool]

Would you mind explain the value you use in random_shape? (In your code, this is (np.random.rand() - 0.5) * np.random.exponential(0.4))，there are some thing confuse me:

1. this value will cause many negative value in simulate ellipse phantom. and the data range lower than shepp logan.
2. The generate random phantom's psnr (after RadonTransform and FBP algorithm) is worse than shepp logan. (as I show in this issue: 31.9 vs 38.4)

Is it possible if I use value of np.random.rand() * np.random.exponential(0.4)? and normalize the random phantom to data range [0, 1] ?

thank you for your good work. ~ @adler-j 👍

[AceCoooool]

After read the shepp logan code. ("force" modify the value in _modified_shepp_logan_ellipsoids). I am curious about the guideline of the code random_phantom.

[adler-j]

this value will cause many negative value in simulate ellipse phantom. and the data range lower than shepp logan.

This is true, but the alternative would cause a situation wherein we only have "positive" ellipses. I guess one option to curb this would be to also enforce positivity of the random phantoms (e.g. apply np.maimum(phantom, 0) or np.abs(phantom) at the end). I opted against this due to added complexity.

The generate random phantom's psnr (after RadonTransform and FBP algorithm) is worse than shepp logan. (as I show in this issue: 31.9 vs 38.4)

This is because (at least for noise-less data, as in your example) shepp-logan is a rather tricky example.

Is it possible if I use value of np.random.rand() * np.random.exponential(0.4)? and normalize the random phantom to data range [0, 1] ?

It is, but I'd recommend just making the phantom positive at the end.