Accelerated reconstruction algorithms are unstable at a number of high iterations

[dreycenfoiles]

When running the accelerated versions of the SIRT and MLEM reconstruction algorithms, the images produced become unstable after a large number of iterations even when the projects are padded with zeros. This can be seen in an example of SIRT using the proj.npy and angle.npy found in the test_data folder.

prj = np.load('proj.npy')
ang = np.load('angle.npy').astype('float32')
pad_width = 50
prj = np.pad(prj, ((0, 0), (0, 0), (pad_width, pad_width)))
norm_sirt = recon(prj, ang, algorithm='sirt',
                  num_iter=100, accelerated=False)
accel_sirt = recon(prj, ang, algorithm='sirt',
                  num_iter=100, accelerated=True) 
plt.imshow(accel_sirt[0]) 
plt.imshow(norm_sirt[0]) 

The above code will produce the following two images. The result of the acceleration algorithm is on top and the non-accelerated is on the bottom.

[carterbox]

The instability comes in the form of the bright spot in the lower left corner near (135, 0). This spot grows in intensity.

[jrmadsen]

Is this also reproducible with other images? Also, was this with nearest-neighbor interpolation? What happens when you use linear interpolation? Cubic interpolation?

What kind of GPU was this run on? Bc I ran far more iterations than 100 when I worked on this and don't really remember an issue along these lines. Lower quality GPUs are more prone to numerical faults which could cause something like this.

[jrmadsen]

I've actually got to do some tomopy benchmarking sometime before Wednesday for a paper so if you send me what I need to run this, I can see if it still shows up on our V100 GPUs.

[carterbox]

@jrmadsen, @dreycenfoiles doesn't have access to a GPU right now, so we are using the OpenCV implementation. proj.npy and angle.npy are from /test/test_tomopy/test_data.

@dreycenfoiles, can see what happens with other interpolation schemes when he needs something to do.