Why isn't FANBEAM linear?

From the test below (R2013a), it appears that the FANBEAM tomographic projection function doesn't satisfy linearity/superposition with high accuracy

    errfunc=@(a,b) norm(a(:)-b(:))/norm(b(:));

    I1=rand(100);
    I2=rand(100);

    F=@(z) fanbeam(z,100);

        >> errfunc(F(I1)+F(I2) ,    F(I1+I2)) %relative error calculation

      ans =

          0.0015

Conversley, the RADON tomographic projector shows a lot more linearity,

    R=@(z) radon(z);

        >> errfunc(R(I1)+R(I2) ,    R(I1+I2))

        ans =

           9.4304e-16

Why is this?

[Matt J]

I went to Tech Support about this and have included their response below. Basically, they attribute the nonlinear behavior to errors in interpolation/resampling. Since interpolation is a linear operation on the image intensity values, I cannot see how this could be the cause, but I wasn't able to convince them. Or I am missing something. Perhaps someone can see what I do not...

From Tech Support:

"After discussing with our development team, I learnt the following:

1) FANBEAM works by first computing the projections with parallel beam geometry and then re-sorting the parallel beam projection data into fan-beam projection data. It uses the re-sorting algorithm specified on page 92 in the book ‘Principles of Computerized Tomographic Imaging” by Kak and Slaney. This resampling of parallel beam projections to fan-beam projections entails a resampling step which requires interpolation. To the best of our understanding it is due to the interpolation step that we see larger values of error for FANBEAM than for RADON.

2) If you are aware of any reference that discusses the effect of interpolation in fan-beam geometry on the superposition property of the projections, please send me a link to the reference.

3) Our development team has done some experiments to sanity-check our FANBEAM implementation and get some sense of whether the source of error is the interpolation step. These experiments are not conclusive. But they do indicate that our FANBEAM implementation is in line with the Kak and Slaney reference. If you have a better test or algorithmic suggestion, we will be interested in hearing that.

Following are a couple of sanity-checks:

a) For large distance value for D (a test you tried as well), the fan-beam geometry is almost like parallel beam geometry. Hence we should expect to see small error in interpolation just like the case for Radon (parallel beam) where there is no interpolation. We do observe the error value go down significantly for large values of D.

b) Use of random images, as shown in your test code, is a particularly tough test case. Since interpolation error can be larger in regions of abrupt change in the image, we can conjecture that if interpolation is the source of superposition error, the error will be lower in smoother images. So we tried a few different images which are not randomly sampled and hence have fewer “edges” or abrupt changes. As expected, we saw smaller error with these real images. The error is of the same order, but smaller nevertheless. For example, we tried using the PHANTOM and 'pout.tif'images:

I1 = phantom('Shepp-Logan',100); I2 = imresize(im2double(imread('pout.tif')),[100 100])

The error went down from 1.469767445052e-03 to 6.149146828217740e-04. We observed such reduction for other examples as well.

Based on these tests, I would say the FANBEAM implementation is correct."

[Matt J]

Sill there in R2017b...

[Matt J]

Sill there in R2019b...