Parallel Radon slightly off for simple example

[carterbox]

Good work on this interesting package. I built 306ae46 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.]

[carterbox]

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.

torch-radon/src/forward.cu

Line 34 in 306ae46

sy = 0.71f * cfg.height;

What is the performance cost of starting the rays at h instead? Why not start the rays at h to prevent these errors?

[matteo-ronchetti]

Could you please try with the V2 branch? It should fix this issue

[carterbox]

Yes, it does, but is V2 ready for release?