Added ND butterworth filter

[din14970]

Description

In high resolution electron microscopy, FFT based filters are commonly used to enhance features in periodic images. One such common filter is a high pass 2D Butterworth filter, which masks low spatial frequencies (usually background) with minimal introduction of artifacts. The Butterworth filter is implemented for 1D signals in scipy but not for 2D arrays. I use scikit-image quite often but was missing this filter, so I hope I can contribute it. For images of the general sort this can be either a low-pass or high-pass filter for spatial frequencies, or it can be combined to create a bandpass filter.

On the astronaut image:

butterworth(rgb2gray(astronaut()), cutoff_frequency_ratio = 0.01, high_pass=False)

butterworth(rgb2gray(astronaut()), cutoff_frequency_ratio = 0.1, high_pass=True)

Not sure whether I put it in a good place and whether it requires its own example

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[pep8speaks]

Hello @din14970! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file skimage/filters/_fft_based.py:

Line 42:13: E126 continuation line over-indented for hanging indent
Line 43:13: E123 closing bracket does not match indentation of opening bracket's line
Line 126:60: W504 line break after binary operator

• In the file skimage/filters/tests/test_fft_based.py:

Line 100:14: E201 whitespace after '['
Line 100:41: E241 multiple spaces after ','
Line 100:55: E241 multiple spaces after ','
Line 100:80: E501 line too long (153 > 79 characters)
Line 101:14: E201 whitespace after '['
Line 101:27: E241 multiple spaces after ','
Line 101:41: E241 multiple spaces after ','
Line 101:80: E501 line too long (153 > 79 characters)
Line 101:97: E241 multiple spaces after ','
Line 101:111: E241 multiple spaces after ','
Line 101:125: E241 multiple spaces after ','
Line 102:14: E201 whitespace after '['
Line 102:27: E241 multiple spaces after ','
Line 102:80: E501 line too long (153 > 79 characters)
Line 102:83: E241 multiple spaces after ','
Line 102:97: E241 multiple spaces after ','
Line 102:111: E241 multiple spaces after ','
Line 102:125: E241 multiple spaces after ','
Line 103:14: E201 whitespace after '['
Line 103:55: E241 multiple spaces after ','
Line 103:69: E241 multiple spaces after ','
Line 103:80: E501 line too long (153 > 79 characters)
Line 103:83: E241 multiple spaces after ','
Line 103:97: E241 multiple spaces after ','
Line 103:111: E241 multiple spaces after ','
Line 103:125: E241 multiple spaces after ','
Line 104:14: E201 whitespace after '['
Line 104:55: E241 multiple spaces after ','
Line 104:80: E501 line too long (153 > 79 characters)
Line 104:97: E241 multiple spaces after ','
Line 104:111: E241 multiple spaces after ','
Line 104:125: E241 multiple spaces after ','
Line 104:139: E241 multiple spaces after ','
Line 104:151: E202 whitespace before ']'
Line 105:27: E241 multiple spaces after ','
Line 105:41: E241 multiple spaces after ','
Line 105:55: E241 multiple spaces after ','
Line 105:69: E241 multiple spaces after ','
Line 105:80: E501 line too long (153 > 79 characters)
Line 105:83: E241 multiple spaces after ','
Line 105:97: E241 multiple spaces after ','
Line 105:111: E241 multiple spaces after ','
Line 105:125: E241 multiple spaces after ','
Line 105:139: E241 multiple spaces after ','
Line 106:27: E241 multiple spaces after ','
Line 106:41: E241 multiple spaces after ','
Line 106:55: E241 multiple spaces after ','
Line 106:69: E241 multiple spaces after ','
Line 106:80: E501 line too long (153 > 79 characters)
Line 106:83: E241 multiple spaces after ','
Line 106:125: E241 multiple spaces after ','
Line 106:139: E241 multiple spaces after ','
Line 107:14: E201 whitespace after '['
Line 107:41: E241 multiple spaces after ','
Line 107:55: E241 multiple spaces after ','
Line 107:69: E241 multiple spaces after ','
Line 107:80: E501 line too long (153 > 79 characters)
Line 107:97: E241 multiple spaces after ','
Line 107:111: E241 multiple spaces after ','
Line 107:139: E241 multiple spaces after ','
Line 108:14: E201 whitespace after '['
Line 108:27: E241 multiple spaces after ','
Line 108:41: E241 multiple spaces after ','
Line 108:55: E241 multiple spaces after ','
Line 108:80: E501 line too long (153 > 79 characters)
Line 108:83: E241 multiple spaces after ','
Line 108:97: E241 multiple spaces after ','
Line 108:125: E241 multiple spaces after ','
Line 108:139: E241 multiple spaces after ','
Line 109:27: E241 multiple spaces after ','
Line 109:41: E241 multiple spaces after ','
Line 109:55: E241 multiple spaces after ','
Line 109:69: E241 multiple spaces after ','
Line 109:80: E501 line too long (153 > 79 characters)
Line 109:83: E241 multiple spaces after ','
Line 109:97: E241 multiple spaces after ','
Line 109:111: E241 multiple spaces after ','
Line 109:139: E241 multiple spaces after ','
Line 123:30: E241 multiple spaces after ','
Line 123:47: E241 multiple spaces after ','
Line 124:14: E201 whitespace after '['
Line 124:30: E241 multiple spaces after ','
Line 124:47: E241 multiple spaces after ','
Line 125:14: E201 whitespace after '['
Line 125:30: E241 multiple spaces after ','
Line 125:47: E241 multiple spaces after ','
Line 126:14: E201 whitespace after '['
Line 126:30: E241 multiple spaces after ','
Line 126:47: E241 multiple spaces after ','
Line 127:14: E201 whitespace after '['
Line 127:30: E241 multiple spaces after ','
Line 128:14: E201 whitespace after '['
Line 129:14: E201 whitespace after '['
Line 137:14: E201 whitespace after '['
Line 137:30: E241 multiple spaces after ','
Line 137:47: E241 multiple spaces after ','
Line 138:14: E201 whitespace after '['
Line 138:47: E241 multiple spaces after ','
Line 140:30: E241 multiple spaces after ','
Line 140:47: E241 multiple spaces after ','
Line 141:30: E241 multiple spaces after ','
Line 141:47: E241 multiple spaces after ','
Line 143:14: E201 whitespace after '['
Line 143:30: E241 multiple spaces after ','
Line 143:47: E241 multiple spaces after ','
Line 144:14: E201 whitespace after '['
Line 144:47: E241 multiple spaces after ','
Line 145:14: E201 whitespace after '['
Line 154:14: E201 whitespace after '['
Line 154:30: E241 multiple spaces after ','
Line 154:47: E241 multiple spaces after ','
Line 156:47: E241 multiple spaces after ','
Line 158:14: E201 whitespace after '['
Line 162:14: E201 whitespace after '['
Line 162:30: E241 multiple spaces after ','
Line 162:47: E241 multiple spaces after ','
Line 163:14: E201 whitespace after '['
Line 163:30: E241 multiple spaces after ','
Line 163:47: E241 multiple spaces after ','
Line 164:14: E201 whitespace after '['
Line 164:30: E241 multiple spaces after ','
Line 164:47: E241 multiple spaces after ','
Line 165:30: E241 multiple spaces after ','
Line 166:30: E241 multiple spaces after ','
Line 167:14: E201 whitespace after '['
Line 167:30: E241 multiple spaces after ','
Line 167:47: E241 multiple spaces after ','
Line 168:30: E241 multiple spaces after ','
Line 169:14: E201 whitespace after '['
Line 169:30: E241 multiple spaces after ','
Line 169:47: E241 multiple spaces after ','
Line 170:30: E241 multiple spaces after ','
Line 171:14: E201 whitespace after '['
Line 171:30: E241 multiple spaces after ','
Line 171:47: E241 multiple spaces after ','
Line 175:14: E201 whitespace after '['
Line 175:30: E241 multiple spaces after ','
Line 175:47: E241 multiple spaces after ','
Line 176:14: E201 whitespace after '['
Line 176:30: E241 multiple spaces after ','
Line 176:47: E241 multiple spaces after ','
Line 177:14: E201 whitespace after '['
Line 177:30: E241 multiple spaces after ','
Line 177:47: E241 multiple spaces after ','
Line 178:14: E201 whitespace after '['
Line 178:30: E241 multiple spaces after ','
Line 178:47: E241 multiple spaces after ','
Line 179:14: E201 whitespace after '['
Line 179:30: E241 multiple spaces after ','
Line 179:47: E241 multiple spaces after ','
Line 180:30: E241 multiple spaces after ','
Line 183:14: E201 whitespace after '['
Line 183:30: E241 multiple spaces after ','
Line 183:47: E241 multiple spaces after ','
Line 188:14: E201 whitespace after '['
Line 188:30: E241 multiple spaces after ','
Line 192:14: E201 whitespace after '['
Line 192:30: E241 multiple spaces after ','
Line 192:47: E241 multiple spaces after ','
Line 193:14: E201 whitespace after '['
Line 193:47: E241 multiple spaces after ','
Line 194:14: E201 whitespace after '['
Line 194:30: E241 multiple spaces after ','
Line 194:47: E241 multiple spaces after ','
Line 195:47: E241 multiple spaces after ','
Line 196:14: E201 whitespace after '['
Line 196:30: E241 multiple spaces after ','
Line 196:47: E241 multiple spaces after ','
Line 201:14: E201 whitespace after '['
Line 201:30: E241 multiple spaces after ','
Line 201:47: E241 multiple spaces after ','
Line 202:14: E201 whitespace after '['
Line 202:30: E241 multiple spaces after ','
Line 202:47: E241 multiple spaces after ','
Line 209:14: E201 whitespace after '['
Line 209:30: E241 multiple spaces after ','
Line 209:47: E241 multiple spaces after ','
Line 212:14: E201 whitespace after '['
Line 214:14: E201 whitespace after '['
Line 214:30: E241 multiple spaces after ','
Line 214:47: E241 multiple spaces after ','
Line 215:14: E201 whitespace after '['
Line 215:30: E241 multiple spaces after ','
Line 215:47: E241 multiple spaces after ','
Line 217:14: E201 whitespace after '['
Line 218:14: E201 whitespace after '['
Line 218:30: E241 multiple spaces after ','
Line 218:47: E241 multiple spaces after ','
Line 219:14: E201 whitespace after '['
Line 219:30: E241 multiple spaces after ','
Line 221:14: E201 whitespace after '['
Line 221:30: E241 multiple spaces after ','
Line 221:47: E241 multiple spaces after ','
Line 225:14: E201 whitespace after '['
Line 225:30: E241 multiple spaces after ','
Line 225:47: E241 multiple spaces after ','
Line 227:30: E241 multiple spaces after ','
Line 227:47: E241 multiple spaces after ','
Line 228:14: E201 whitespace after '['
Line 228:30: E241 multiple spaces after ','
Line 228:47: E241 multiple spaces after ','
Line 229:14: E201 whitespace after '['
Line 229:30: E241 multiple spaces after ','
Line 229:47: E241 multiple spaces after ','
Line 231:14: E201 whitespace after '['
Line 231:30: E241 multiple spaces after ','
Line 231:47: E241 multiple spaces after ','
Line 232:14: E201 whitespace after '['
Line 232:30: E241 multiple spaces after ','
Line 232:47: E241 multiple spaces after ','
Line 234:14: E201 whitespace after '['
Line 234:30: E241 multiple spaces after ','
Line 234:47: E241 multiple spaces after ','
Line 237:14: E201 whitespace after '['
Line 237:30: E241 multiple spaces after ','
Line 237:47: E241 multiple spaces after ','
Line 238:14: E201 whitespace after '['
Line 241:30: E241 multiple spaces after ','
Line 241:47: E241 multiple spaces after ','
Line 243:14: E201 whitespace after '['
Line 243:30: E241 multiple spaces after ','
Line 243:47: E241 multiple spaces after ','
Line 244:14: E201 whitespace after '['
Line 244:30: E241 multiple spaces after ','
Line 244:47: E241 multiple spaces after ','
Line 246:14: E201 whitespace after '['
Line 246:30: E241 multiple spaces after ','
Line 246:47: E241 multiple spaces after ','

Comment last updated at 2021-06-20 17:07:46 UTC

[grlee77]

I thought about this but was unsure how to best approach this and what to consider an n-dimensional image. In an RGB image, the channels should each independently be 2D Fourrier transformed and filtered, the channel dimension should not be considered a spatial dimension. But in a gray level 3D image of course the third dimension should be considered. Perhaps then I need a multichannel keyword arg. I'll see if I can make a decent implementation.

Historically we have used multichannel=True to indicate that the last axis contains channels rather than a spatial dimension. However we are currently in the process of deprecating that and moving to a channel_axis keyword argument. This argument should default to channel_axis=None, meaning all axes are spatial. It is sufficient to support only channel_axis=None and channel_axis=-1 if the general case would be more of a pain to implement. We can just add a channels_as_last_axis decorator to automatically do a temporary shift of channels to the last axis as needed.

[din14970]

Ok, was a bit more complicated than I expected to make channel work for any axis so it looks a bit hacky. N dimensional images are supported and any axis index can be passed as channel index. By default it is assumed all axes are spatial dimensions, a channel_axis MUST be passed for color images. Any number of channels are supported. While refactoring I also slightly improved the mask generation function.

For large n-dimensional images the creation of the Fourier mask may pose a memory issue due to the reliance on meshgrid. More memory friendly may be a "loopy" implementation but for this I would need numba as I'm not familiar with cython. Are there any plans to allow numba as a dependency?

[din14970]

Also, I feel like a large portion of the butterworth filter could be recycled for other types of FFT masks. Perhaps if there is interest in the future this could be pulled out and generalized.

[grlee77]

Thanks @din14970. I left a first round of review comments here. A number of just simple style suggestions for the docstrings.

Note that you don't have to commit each suggestion individually, but can add them all to a single batch from the "Files changed" tab here on GitHub. Alternatively, you can manually apply them offline and push a commit the traditional way.

[grlee77]

Let me think about other possible tests for correctness.

The current tests are testing that things run without crashing for a range of dimensions / channel axes, but aren't really checking any property of the filtered image or comparing to a reference result.

[din14970]

Yes I didn't have any good ideas on this. I thought first of constructing a dummy FFT by creating a zero array and then initializing a few pixels to 1, inverse fourrier transforming, then passing this image into the bw filter, then fourrier transforming it back. Unfortunately this no longer resembles the original FFT. I guess if one wisely constructs the FFT it should work?

[grlee77]

I came up with another basic test case that just verified highpass vs. lowpass behavior, but want to resolve the one remaining issue regarding preserve_range first.

[din14970]

Also not entirely sure why the CI is failing, something to do with the docs I suppose

[grlee77]

Nice job updating this to use the rfftn and irfftn functions!

In searching for references to use of this function in image processing, I sometimes saw the filter defined without any square root in the denominator. For example here, or in John C. Russ's "Image Processing Handbook" (3rd. ed.), where the filter equation is listed as:
H = 1 / (1 + C * (R / R0) ** (2 * n)) where it is stated that C is usually either 1.0 or 0.414.

Other sources such as the Wikipedia article and this one or this one seem to agree with the form implemented here.

Perhaps we should make the exact equation explicit in the Notes of the docstring here to avoid potential confusion over these differing conventions

[jni]

The nd implementation is gorgeous @din14970! regarding the memory concerns for large nD images, we typically leave that to later generations to address. 😂 It's much better to have a working and slow implementation merged to get further improvements down the line as needed.

There's a couple of comments from @grlee77 that need addressing (order/sqrt conventions), as well as the improved testing idea, but other than that this is ready for me.

Regarding testing, in many situations, a test is hard to define from first principles. In those cases, we typically use one or more of the following options:

• Run this for 1D arrays and for "2D multichannel" arrays (ie a series of 1D arrays), and compare it to the output from SciPy.
• Run this on a known image (e.g. a tiny crop data.camera()), check that the output makes sense, save it as a .npy file, and check against that. There's a few examples in e.g. skimage/restoration/tests. Make this small so that they don't take up too much space. Converting to uint8 is also a useful trick for compression.
• Do the same as above but compare to a 3rd party library.

[din14970]

@grlee77 I tried to dig into this a bit but am still not entirely sure.

I thought first it might be a dimension issue. The references where there is a square root in the denominator always pertain to 1D signals; whenever it concerns 2D signals there is no square root. I couldn't find any reference explaining the difference however. So my hypothesis was that one could extrapolate to any dimension with 1 / (1 + d ** (2*n)) ** (k/2) with k the number of spatial dimensions. To test this out I made the following example:

import scipy.fft as fft
import numpy as np
import matplotlib.pyplot as plt

def funk(R):
    return R

w1 = 0.1 
order = 8
# 1D
x = np.arange(100)
y1 = funk(x)
ify1 = butterworth(y1, w1, False, order, 0.5)
ify1b = butterworth(y1, w1, False, order, 0.5)

# 2D
X1, X2 = np.meshgrid(x, x, sparse=True)
R = np.sqrt(X1**2 + X2**2)
y2 = funk(R)
ify2 = butterworth(y2, w1, False, order, 0.5)
ify2b = butterworth(y2, w1, False, order, 1)

# 3D
X1, X2, X3 = np.meshgrid(x, x, x, sparse=True)
R = np.sqrt(X1**2 + X2**2 + X3**2)
y3 = funk(R)
ify3 = butterworth(y3, w1, False, order, 0.5)
ify3b = butterworth(y3, w1, False, order, 1.5)

fig, ax = plt.subplots(ncols = 3, figsize=(15, 5))
ax[0].plot(x, y1, label="1D")
ax[0].plot(x, y2[0], label="2D")
ax[0].plot(x, y3[0,0], label="3D")
ax[0].legend()
ax[0].set_title("Original function")

ax[1].plot(x, ify1, label="1D")
ax[1].plot(x, ify2[0], label="2D")
ax[1].plot(x, ify3[0,0], label="3D")
ax[1].legend()
ax[1].set_title("Butter filter exponent 0.5")

ax[2].plot(x, ify1b, label="1D")
ax[2].plot(x, ify2b[0], label="2D")
ax[2].plot(x, ify3b[0,0], label="3D")
ax[2].legend()
ax[2].set_title("Butter filter adjusted exponent")

Basically I generate some N-D images based on some function of R on a 1D, 2D and 3D grid. Looking along a single axis, all these functions look identical (figure 1). I then apply the butterworth filter and plot again the result along one dimension. I modified the butterworth function to also take the parameter exponent: wfilt = 1 / (1 + np.power(q2, order))**exponent. The result:

The middle image is where the exponent is 0.5 for all dimensions, the right one is where the exponent is adjusted according my hypothesis. Basically it looks like it doesn't matter very much. I looked into different functions, different orders, different w0 but couldn't find a clear reason to use one or the other.

Therefore, since the primary use-case for this thing would be to highlight/remove certain frequencies in an image, and since a user will have to tune parameters on a case by case basis anyway, I propose we go with exponent 1 in the denominator for simplicity and speed. The filter has sufficient flexibility with the order and cutoff parameters.

@jni I think I'll go with the second approach, though this breaks down if the filter is ever changed for some reason (i.e. people want the square root in the denominator)

[din14970]

@grlee77 @jni the tests are failing because it can't find the test files I added; locally it works fine. What is the proper way to do this?

[grlee77]

The doc errors were due to LaTeX formatting issues, but this isn't that clear from the message. The solution is using \\ to get \ in the docstring's LaTeX code. There was also a typo in the equation with a ) instead of } at the end of the line.

[grlee77]

f53349f closes #5416, which was initially preventing me from building the docs locally. It is unrelated to this PR, though, so perhaps I should open a standalone PR for it

[din14970]

@grlee77 awesome thanks a lot for looking into this! Weird because I also considered this but I saw multiple examples e.g.

scikit-image/skimage/metrics/simple_metrics.py

Line 200 in 178e05f

Y(A, B) = \frac{H(A) + H(B)}{H(A, B)}

where a single \ appears to be enough. Anyway with these changes it seems to be working!

[grlee77]

Yes, but I think those cases start the docstring with r""" instead of """ (i.e. it is a Python raw string). We could also take that approach here to make the LaTeX more readable.

[grlee77]

Suggested change

	butterfilt = fft.irfftn(wfilt * fft.rfftn(image, axes=axes),
	s=fft_shape, axes=axes)
	out_range = dtype_limits(image) if preserve_range else (0.0, 1.0)
	butterfilt = rescale_intensity(butterfilt, out_range=out_range)
	if not preserve_range:
	image = img_as_float(image)
	butterfilt = fft.irfftn(wfilt * fft.rfftn(image, axes=axes),
	s=fft_shape, axes=axes)

I would have expected behavior like the above for preserve_range. The current code will rescale the range for pretty much all inputs even when preserve_range is True. For example, floating point inputs will not have their range preserved.

@jni, thoughts? If we change this here we will need to change the data files used in the tests.

[jni]

Gah. I think you're right @grlee77. 🤦

[din14970]

sorry for replying so late, didn't have time for a couple days. I think I implemented it in this way because the fft + ifft may introduce image intensities outside the original data-type range. For example if you put in a uint8 image, you may get negative intensities. This is because when you remove some spatial frequencies, you might get "jumps" in places where there are edges. My original interpretation of preserve_range was that the output image should fit in the same dtype range as the input image so that it could also be cast as the same dtype. This is why I do the rescaling after the filter. But I now see your point that this has the potentially unwanted side effect of stretching the histogram to the extrema, so probably your solution is better. Note however that in your implementation, it is still possible to get an image that is not rescaled to a range [0, 1] or [-1, 1] if preserve_range=False, but I don't see another appropriate alternative that won't stretch the histogram. I suggest the following:

• I change the implementation to what @grlee77 suggests
• I add some clarification in the Notes section about preserve_range that it is possible to still get some outlier intensities
• I remove the test data files and hardcode new test data into the test files.

Alternatively, since the option doesn't seem to make much sense anyway I remove the preserve_range argument. Any vote on this?

[din14970]

I decided to get rid of preserve_range as it does not make sense in my opinion. Even if input is scaled between 0-1, the output still will not be. It just adds confusion to what the argument is supposed to do, and depending on what users want it is easier if they just rescale the image intensity before or after the filter. I removed the test files and just copied the expected arrays into the tests. I merged again with main, now something is failing again.

[grlee77]

I merged again with main, now something is failing again.

These seem unrelated and should be resolved by #5434

[grlee77]

I am okay with having these relatively small arrays in-line in the test file, but let's see what other maintainers think. I updated setup.cfg to ask it to ignore some whitespace and line-length issues that would be reported during style checking.

Please see the two new test cases added in 0a769f7 and see if they look reasonable to you. I wanted to test something with non-zero input that would at least verify that the spectrum of the output had been modified in a roughly high-pass or low-pass nature. The specific 0.1 and 0.9 thresholds are basically arbitrary and not tuned to the specific filter settings.

[din14970]

Thanks for adding some tests, they are definitely valuable. I just left two minor comments which may just be misunderstandings on my part.

[mkcor]

@din14970 CI should pass if you merge main now 👍

[mkcor]

I am okay with having these relatively small arrays in-line in the test file, but let's see what other maintainers think.

I am also okay with this practice. Technically, this choice reverts #5415 but this is fine and was contemplated in the first place.

I updated setup.cfg to ask it to ignore some whitespace and line-length issues that would be reported during style checking.

Wonderful!

[grlee77]

Okay, the one failure here is unrelated and seems to be a new warning from the just released scipy 1.7.0. I can make a new PR to fix it

E                   ValueError: Unexpected warning: scipy.linalg.pinv2 is deprecated since SciPy 1.7.0, use scipy.linalg.pinv instead

[grlee77]

Thanks for sticking with this @din14970, it looks good to me now.

[din14970]

@grlee77 @mkcor @jni thanks a lot also for your thorough reviews, suggestions and improvements!