torch-radon
torch-radon copied to clipboard
matteo-ronchetti
Parallel Radon slightly off for simple example
Good work on this interesting package. I built 306ae4638f6c8a2ea4621a857e03419174d77df4 and ran the following simple test where I expect that the Radon operator should be perfect at angles 0 and PI/2.
It seems that possibly the angle specification is wrong or there is some other bias in the geometry for the rays. Was wondering if there is an explanation for this and if it could be fixed.
import torch
from torch_radon import Radon
def test_simple_integrals(image_size=16):
"""Check that the forward radon operator works correctly at 0 and PI/2.
When we project at angles 0 and PI/2, the foward operator should be the
same as taking the sum over the object array along each axis.
"""
angles = torch.tensor(
[0.0, -np.pi / 2, np.pi, np.pi / 2],
dtype=torch.float32,
device='cuda',
)
radon = Radon(
resolution=image_size,
angles=angles,
det_spacing=1.0,
det_count=image_size,
clip_to_circle=False,
)
original = torch.zeros(
image_size,
image_size,
dtype=torch.float32,
device='cuda',
)
original[image_size // 4, :] += 1
data = radon.forward(original)
data0 = torch.sum(original, axis=0)
data1 = torch.sum(original, axis=1)
print('\n', data[0].cpu().numpy())
print(data0.cpu().numpy())
print('\n', data[1].cpu().numpy())
print(data1.cpu().numpy())
print('\n', data[2].cpu().numpy())
print(data0.cpu().numpy()[::-1])
print('\n', data[3].cpu().numpy())
print(data1.cpu().numpy()[::-1])
tests/test_parallel_beam.py::test_simple_integrals
[0.99632365 0.99632365 0.99632365 0.99632365 0.99632365 0.99632365
0.99632365 0.99632365 0.99632365 0.99632365 0.99632365 0.99632365
0.99632365 0.99632365 0.99632365 0.99632365]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[ 0. 0. 0. 0. 16.000002 0. 0.
0. 0. 0. 0. 0. 0. 0.
0. 0. ]
[ 0. 0. 0. 0. 16. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0.99632365 0.99632365 0.99632365 0.99632365 0.99632365 0.99632365
0.99632365 0.99632365 0.99632365 0.99632365 0.99632365 0.99632365
0.99632365 0.99632365 0.99632365 0.99632365]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[ 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 16.000002 0. 0.
0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 16. 0. 0. 0. 0.]
This precision error seems to be related to starting the rays too close to the center of rotation. It looks like the code is trying to place the starting ray on the circle inscribing the square i.e. sqrt(2)/2*h.
https://github.com/matteo-ronchetti/torch-radon/blob/306ae4638f6c8a2ea4621a857e03419174d77df4/src/forward.cu#L34
What is the performance cost of starting the rays at h instead? Why not start the rays at h to prevent these errors?
