Update radon_transform.py

[jcesardasilva]

I have notice some small mistakes in the filters of the iradon in scikit-image/skimage/transform/radon_transform.py, in lines 208 up to 216. In some part where the division by 2 should exist, it does not, and vice-versa. See below the code lines with the corrections.

in scikit-image/skimage/transform/radon_transform.py, in lines 208 up to 216

elif filter == "shepp-logan":
    # Start from first element to avoid divide by zero
    fourier_filter[1:] = fourier_filter[1:] * np.sin(omega[1:]/2) / (omega[1:]/2)
elif filter == "cosine":
    fourier_filter *= np.cos(omega/2)
elif filter == "hamming":
    fourier_filter *= (0.54 + 0.46 * np.cos(omega))
elif filter == "hann":
    fourier_filter *= (1 + np.cos(omega)) / 2
elif filter is None:
    fourier_filter[:] = 1

This becomes more evident when using frequency cutoffs.

code

elif filter_type == "shepp-logan":
    # Start from first element to avoid divide by zero
    fourier_filter[1:] = fourier_filter[1:] * (np.sin(omega[1:]/(2*freqcutoff)) / (omega[1:]/(2*freqcutoff)))
elif filter_type == "cosine":
    fourier_filter *= np.cos(omega/(2*freqcutoff))
elif filter_type == "hamming":
    fourier_filter *= (0.54 + 0.46 * np.cos(omega / (freqcutoff)))
elif filter_type == "hann":
    fourier_filter *= (1 + np.cos(omega / (freqcutoff))) / 2

I also change the code to enable the reconstruction from differential projections like the ones taken with grating interferometer or to be insensitive to the big phase-wraps of ptychographic X-ray computed tomography dataset

[sciunto]

Is it possible to reflect that in the test suite?

[sciunto]

We don't specify default values.

[sciunto]

Could you also clean up above (remove braces after optional)?

[sciunto]

Also, above, output_size : int**, optional**

[emmanuelle]

@sciunto well this is not completely true :-). See for example transform.probabilistic_hough_line where default values are given just after the type. It's the same for transform.radon and transform.iradon. I don't have strong feelings about this, but I think it can improve the readability of a docstring to put the default values on the same line as the parameter name.

[emmanuelle]

Is freqcutoff smaller or larger than 1? Please clarify.

[sciunto]

To me, these are mistakes and I had in mind that the standard is to do not specify values if there is signature inspection.
I base my judgement on this: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
and we follow numpy's standards.

[sciunto]

Thanks for reporting and fixing this! I made few comments on the style and we would appreciate more tests on iradon to avoid that in the future :)

[sciunto]

spaces around * and / operators

[emmanuelle]

yes, sorry to be nitpicky about this @jcesardasilva but we try to enforce the PEP8 guidelines as possible (see https://www.python.org/dev/peps/pep-0008/) in order to have style consistency throughout the package.

[sciunto]

We all have our favourite nitpicky stuffs :P

[sciunto]

@scikit-image/core Is there anyone who is expert with radon who can approve these changes?

@jcesardasilva Kind ping. Do you have time to make the changes I suggested? Your PR also needs a rebase.

[jcesardasilva]

@sciunto what did you mean by "PR needs to be rebase"? This modified version of FBP has been used in the literature for years.

[sciunto]

@jcesardasilva It's because master evolved in parallel to your PR and there are merging conflicts that cannot be solved automatically. You need to do it manually. I can do it for you and make a PR against your branch if necessary.

See
http://scikit-image.org/docs/0.12.x/gitwash/development_workflow.html#rebase-on-trunk

[emmanuelle]

Hi @jcesardasilva we're getting there :-). What remains to be done before we can merge your PR:

• solve the conflicts with master. As explained by @sciunto you need to rebase on master. See also the link http://scikit-image.org/docs/stable/contribute.html#divergence-between-upstream-master-and-your-feature-branch to know how to rebase.
• add unit tests in transform/tests/test_radon_transform.py. Since you added features (the frequency cutoff and the possibility to use derivatives), they need to be tested. Also, since you corrected some bugs, it would be great to write a test that failed before your correction, and that passes now. I don't know if it's possible (for example by applying radon and then iradon on an image??).
• very small changes mentioned as comments (PEP8 etc.)

[jcesardasilva]

Ok, the time is quite short these days and I cannot do much, but I can give you already an idea of the test I made when I figured out the problem with the filter. The rest I do later, but I cannot do now. All the help is welcome.
You can run the code below to see that the filters of the old version were not right, which is fixed in the new version. For freqcutoff =1 which mean 100% of passing frequencies , i.e., no cutoff, you come to the version without cutoff. But you can also play with different cutoffs.
`

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftfreq
#from scipy.interpolate import interp1d

from skimage.io import imread
from skimage import data_dir
from skimage.transform import radon, rescale

freqcutoff = 1# 0.3

image = imread(data_dir + "/phantom.png", as_grey=True)
image = rescale(image, scale=0.4)
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta, circle=True)

projection_size_padded =  max(64, int(2 ** np.ceil(np.log2(2 * sinogram.shape[0]))))

# Construct the Fourier filter
f = fftfreq(projection_size_padded).reshape(-1, 1)#*freqcutoff   # digital frequency   
omega = 2 * np.pi * f                                # angular frequency
fs = np.fft.ifftshift(f)

plt.close('all')
fig1 = plt.figure(1)
ax1 = fig1.add_subplot(111)
fig2 = plt.figure(2)
ax2 = fig2.add_subplot(111)

# plot Ram-Lak filter
######### old
fourier_filter_rl = 2 * np.abs(f)         # ramp filter
fourier_filter_rl[np.where(2*np.abs(f)>freqcutoff)]=0
ax1.plot(fs,np.fft.ifftshift(fourier_filter_rl))
######### new # the same as old in this case
fourier_filter_rl2 = 2 * np.abs(f)         # ramp filter
fourier_filter_rl2[np.where(2*np.abs(f)>freqcutoff)]=0
ax2.plot(fs,np.fft.ifftshift(fourier_filter_rl2))

# plot Shepp-Logan filter
######### old
fourier_filter_sl = 2 * np.abs(f) 
# Start from first element to avoid divide by zero
fourier_filter_sl[1:] = fourier_filter_sl[1:] * (np.sin(omega[1:]/(freqcutoff)) / (omega[1:]/(freqcutoff)))
fourier_filter_sl[np.where(2*np.abs(f)>freqcutoff)]=0
ax1.plot(fs,np.fft.ifftshift(fourier_filter_sl))
######### new
fourier_filter_sl2 = 2 * np.abs(f) 
# Start from first element to avoid divide by zero
fourier_filter_sl2[1:] = fourier_filter_sl2[1:] * (np.sin(omega[1:]/(2*freqcutoff)) / (omega[1:]/(2*freqcutoff)))
fourier_filter_sl2[np.where(2*np.abs(f)>freqcutoff)]=0
ax2.plot(fs,np.fft.ifftshift(fourier_filter_sl2))

# plot cosine filter
######### old
fourier_filter_cos = 2 * np.abs(f) 
fourier_filter_cos *= np.cos(omega/(freqcutoff))
fourier_filter_cos[np.where(2*np.abs(f)>freqcutoff)]=0
ax1.plot(fs,np.fft.ifftshift(fourier_filter_cos))
######### new 
fourier_filter_cos2 = 2 * np.abs(f) 
fourier_filter_cos2 *= np.cos(omega/(2*freqcutoff))
fourier_filter_cos2[np.where(2*np.abs(f)>freqcutoff)]=0
ax2.plot(fs,np.fft.ifftshift(fourier_filter_cos2))

# plot hamming filter
######### old
fourier_filter_hamm = 2 * np.abs(f) 
fourier_filter_hamm *= (0.54 + 0.46 * np.cos(omega / (2*freqcutoff)))
fourier_filter_hamm[np.where(2*np.abs(f)>freqcutoff)]=0
ax1.plot(fs,np.fft.ifftshift(fourier_filter_hamm))
######### new
fourier_filter_hamm2 = 2 * np.abs(f) 
fourier_filter_hamm2 *= (0.54 + 0.46 * np.cos(omega / (freqcutoff)))
fourier_filter_hamm2[np.where(2*np.abs(f)>freqcutoff)]=0
ax2.plot(fs,np.fft.ifftshift(fourier_filter_hamm2))

# plot hann filter
######### old
fourier_filter_hann = 2 * np.abs(f) 
fourier_filter_hann *= (1 + np.cos(omega / (2*freqcutoff))) / 2
fourier_filter_hann[np.where(2*np.abs(f)>freqcutoff)]=0
ax1.plot(fs,np.fft.ifftshift(fourier_filter_hann))
######### new
fourier_filter_hann2 = 2 * np.abs(f) 
fourier_filter_hann2 *= (1 + np.cos(omega / (freqcutoff))) / 2
fourier_filter_hann2[np.where(2*np.abs(f)>freqcutoff)]=0
ax2.plot(fs,np.fft.ifftshift(fourier_filter_hann2))

ax1.set_title('Old version')
ax1.set_xlabel('Frequency')
ax1.set_ylabel('Filter')
ax1.axis('tight')
ax1.legend(['Ram-Lak','Shepp-Logan','Cosine','Hamming','Hanning'],loc='best')
ax2.set_title('New version')
ax2.set_xlabel('Frequency')
ax2.set_ylabel('Filter')
ax2.axis('tight')
ax2.legend(['Ram-Lak','Shepp-Logan','Cosine','Hamming','Hanning'],loc='best')

plt.show(block=False)
a = raw_input('Type <Enter> to close all figures')
plt.close('all')

`

[stefanv]

Also see https://github.com/scikit-image/scikit-image/pulls?utf8=%E2%9C%93&q=is%3Apr+is%3Aopen+radon