search

phernst / pytorch_radon

Pytorch implementation of scikit-image's radon function and more
MIT License
25 stars 11 forks source link

Comparison with Scikit-image Radon #1

Open kouximei opened 3 years ago

kouximei commented 3 years ago

Hi! Many thanks for this repo! It is very helpful! I'm wondering if there is any difference with scikit-image's radon? For instance, if I run the following code, I would get a small difference between pytorch_radon and scikit-image's radon (around 1e-6 ~ 1e-8) . Is this because these two radon functions are implemented in different ways? Or is this because there's precision loss in torch? If this difference is due to precision loss, do you know of any ways to fix that? Thank you!

image

Follow up: I just fixed this precision issue, just need to change two lines of the utils.py like this.

image
phernst commented 3 years ago

hi! check out the development branch. i just updated the constants to use torch.double instead of the default torch.float. the difference between pytorch_radon and scikit's radon should be minimal now. let me know if it works for you :)

kouximei commented 3 years ago

Thanks a lot for your commit! It works very well for forward and back projection!

BTW, I compared pytorch_radon dev branch and scikit's radon's hann filtered backprojection. I'm using torch.double and the mean absolute difference between them is still around 1e-6. I checked through the codes and I'm not sure if this is still a precision loss any more. Did you implement the hann back projection in a different way from scikit?

P.S. That's not a big difference, but I'm just very curious why that happened, because everything seems to be using torch.double now.

phernst commented 3 years ago

Hm, might be related to ramp_filter in filters.py. I am going to check this

phernst commented 3 years ago

ok, it is only partly related to ramp_filter: the variable image_filter should be torch.double. but another quite large error comes from torch.irfft. these are the errors to scikit's hann iradon for the different combinations of (r)fft and i(r)fft:

fourier inverse fourier error
fft ifft 4.2482147420681154e-11
rfft ifft 4.2482147117965245e-11
fft irfft 2.697666179063214e-06
rfft irfft 2.697666179063532e-06

i can't explain this behavior but i will use the combination of rfft/ifft and commit it shortly

kouximei commented 3 years ago

I see! Thanks for locating this ifft problem, it totally slipped my mind... Maybe that's due to the the problem described here.https://github.com/pytorch/pytorch/issues/19111, but I'm not exactly sure. I guess we don't need to worry about that.

kouximei commented 3 years ago

I notice another strange thing. Right now, in the test_radon.py, the phantom is first filled with zeros, and then the values of its center part are changed to 0.5. And it worked very well. However, instead of first filling it with zeros, if I first filled it with ones, or any values other than zero, and then changed the center part to 0.5, then the back projection's pixel at [0,0] alone differs a lot from the back projection from scikit (around 1e-3). The rest of the pixels differ very little from scikit (about 1e-14). I'm a beginner on the radon transform and don't quite understand the reason behind this. Do you know what's happening here? Many thanks!

Actually, when I replace the testing phantom with any other real CT images, I can observe the same phenomena: the single pixel at [0,0] has a large difference from scikit, the rest pixels are the same.

I've tested more phantoms on this problem: It seems that if I set the first two rows or the first two columns to zeros, then no matter how I changed the rest of the pixels of the phantom, the resulting back projections will be the same as scikit (difference within 1e-14). But I still don't figure out why is that...

phernst commented 3 years ago

I see! Thanks for locating this ifft problem, it totally slipped my mind... Maybe that's due to the the problem described here.pytorch/pytorch#19111, but I'm not exactly sure. I guess we don't need to worry about that.

I also thought so first, but the problem there was that he used onesided=True and didn't specify the shape in irfft. However, when using onesided=False (as I did), ifft and irfft should have basically the same output (I got a difference of 1e-14 using torch.double). So something else must go wrong.

phernst commented 3 years ago

then the back projection's pixel at [0,0] alone differs a lot from the back projection from scikit (around 1e-3). The rest of the pixels differ very little from scikit (about 1e-14).

Can you show me your code?

kouximei commented 3 years ago

Here's my code. Basically all I did is changing the background of the phantom from zero to a non-zero number.

image

from pytorch_radon import Radon, IRadon import torch import skimage.transform import numpy as np import sys

set cuda device and display limit

np.set_printoptions(threshold=sys.maxsize, precision=20) torch.set_printoptions(precision=10) device = torch.device('cuda:4')

create phantom

For the following line, I changed the torch.zeros to torch.ones.

img = torch.ones(1, 1, 256, 256, dtype=torch.double) * 1 img[:, :, 120:130, 120:130] = 0.5 img = img.to(device).requiresgrad(True) img_np = np.float64(np.squeeze(img.detach().cpu().numpy()))

set theta

theta = torch.arange(180, dtype=torch.double).to(device).requiresgrad(True) theta_np = np.float64(theta.detach().cpu().numpy())

forward project (torch and scikit)

r = Radon(img.shape[2], theta, False, dtype=torch.double, device=device).to(device) sino = r(img) sino_np_torch = np.squeeze(sino.detach().cpu().numpy()) sino_np_scipy = skimage.transform.radon(img_np, theta=theta_np, circle=False)

backward project (torch and scikit)

ir = IRadon(img.shape[2], theta, False, dtype=torch.double, use_filter=None, device=device).to(device) reco = ir(sino) reco_np_torch = np.squeeze(reco.detach().cpu().numpy()) reco_np_scipy = skimage.transform.iradon(sino_np_scipy, theta=theta_np, filter=None, circle=False)

print the difference between the sinogram

print('comparing sinogram') print('sino max difference ' + str(np.max(np.abs(sino_np_scipy - sino_np_torch)))) print('sino mean difference ' + str(np.mean(np.abs(sino_np_scipy - sino_np_torch)))) print()

print the difference between the reconstruction

print('comparing reconstruction') print('reco max difference except pixel at index [0, 0] ' + str(np.max(np.abs(reco_np_scipy - reco_np_torch).flatten('F')[1:]))) print('reco mean difference except pixel at index [0, 0] ' + str(np.mean(np.abs(reco_np_scipy - reco_np_torch).flatten('F')[1:]))) print('index of pixel with large error ' + str(np.argwhere(np.abs(reco_np_scipy - reco_np_torch) > 1e-12))) print('error of pixel at index [0, 0]: ' + str((reco_np_scipy - reco_np_torch)[0, 0]))

phernst commented 3 years ago

I think I found out where the errors come from. The problem is related to the grids that are created and used for the interpolation (both radon and iradon). The center of rotation in scikit's implementation is shape//2, pytorch uses the actual center of the image (with subpixel accuracy). This means scikit produces a small error if the shape of the image is even. However, it is possible that pytorch_radon produces another error because of align_corners=True in grid_sample and affine_grid. Setting it to True doesn't make sense to me mathematically but it produces smaller reconstruction errors than setting it to False, at the moment. I will check this again.

phernst commented 3 years ago

align_corners=True is correct. For Radon, there is actually no difference but in IRadon, the grids that are created for grid_sample assume the centers of the extreme pixels to be -1 and 1, which is the definition of align_corners=True.

kouximei commented 3 years ago

I see! Thanks a lot! Could you please push the changes to the repo?

phernst commented 3 years ago

All the changes are now in the master branch :)

mingyangx commented 3 years ago

However, it is possible that pytorch_radon produces another error because of align_corners=True in grid_sample and affine_grid. Setting it to True doesn't make sense to me mathematically but it produces smaller reconstruction errors than setting it to False, at the moment. I will check this again.

align_corners=True is correct. For Radon, there is actually no difference but in IRadon, the grids that are created for grid_sample assume the centers of the extreme pixels to be -1 and 1, which is the definition of align_corners=True.

I'm actually a bit confused here. Because Pytorch_radon has been using align_corners=True since the beginning, right? So basically that the difference at pixel [0, 0] between torch and scikit is due to the different implementation of the center of rotation, and the scikit is inaccurate because it uses shape//2 instead of the actual center, am I understanding this correctly?

phernst commented 3 years ago

I think so, yes. I just noticed something else. In your code above, try to change theta to not include 135 and the difference at [0, 0] will be around 1e-13. In fact, including 135 (or 135+k*180) in theta creates this big difference at [0, 0]. Set theta = [135, 135] (twice to avoid dividing by 0) to get the error or set theta to another value (except 135+k*180) to get around 1e-13. I still assume that this is somehow related to shape//2 in scikit.