Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Gaussian filter truncation #239

Closed
wants to merge 2 commits into from

6 participants

@thouis

These changes add control over the radius of truncation of the Gaussian filters in scipy.ndimage, and add the ability to pass keyword arguments through gaussian_laplace()/gaussian_gradient_magnitude() to gaussian_filter().

This is an expanded form of roderikk's patch from this ticket:
http://projects.scipy.org/scipy/ticket/1660

thouis added some commits
@thouis thouis ENH: expose control truncation of ndimage.gaussian_filter().
This change gives the user control over the size at which the Gaussian filter is truncated,
with default behavior to truncate at 4 standard deviations (the previous fixed value).
482ec0d
@thouis thouis ENH: Pass extra keyword arguments to gaussian_laplace() and gaussian_…
…gradient_magnitude() to gaussian_filter()

This change allows control over the truncation radius of the gaussian filter used in
gaussian_laplace() and gaussian_gradient_magnitude().
933e8b5
@thouis

roderikk's original patch used flen for the new argument to gaussian_filter(). I thought this was misleading in the 2D case, as well as possibly being confusing for being in units of standard deviation.

But I don't like truncate much more.

Perhaps the truncation amount should be specified in pixels. Also, should it be controllable in each dimension separately in gaussian_filter(), in the same way that sigma is?

@rgommers
Owner

Anyone want to have a look at this PR this week? Branching 0.12.x this weekend.

@roderikk

Hi, I contributed the original bug in Trac. How can we get this pulled into main?

I think that the truncation as a multiple of sigma is very logical as it gives you a clear estimate of the amount of information thrown away.

For me there is no need to to do a different value in each direction as that would only lead to very weird aliasing results.

@WarrenWeckesser
Collaborator

@roderikk: I'll take a look at this today.. I've been working on an enhancement for ndimage that provides a gaussian filter with an arbitrary covariance matrix, so I've been thinking about this issue lately.

@roderikk

That sounds interesting. If you have some code to be tested let me know.

@WarrenWeckesser
Collaborator

@thouis: This looks like a good change. I prefer your 'truncate' over the original 'flen', and I don't have any better suggestions. (Anyone else have an opinion?) I like that you've allowed gaussian_laplace and gaussian_gradient_magnitude to also use the argument, via **kwargs.

I've made a pull request on your repo with some changes to the test: thouis#1
I figured that would be quicker than suggesting a lot of changes with inline comments. Take a look, and if those changes are OK with you, merge it into this PR. Then if no one else has any objections or big changes to suggest, we might be able to get this into the 0.13 release.

@pv
Owner
pv commented

Continued in gh-2767

@pv pv closed this
@BenFrantzDale

I realize this is closed, but I've been thinking about a related issue: as sigma goes to zero, the default behavior is ugly: when s goes below 0.25, gaussian_filter suddenly does nothing. Since performance isn't really an issue with a kernel size of three pixels, I would prefer that the kernel not become a delta function until sigma is exactly zero. Thoughts?

@pv
Owner
pv commented

@BenFrantzDale: please open a new issue, thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Jun 5, 2012
  1. @thouis

    ENH: expose control truncation of ndimage.gaussian_filter().

    thouis authored
    This change gives the user control over the size at which the Gaussian filter is truncated,
    with default behavior to truncate at 4 standard deviations (the previous fixed value).
  2. @thouis

    ENH: Pass extra keyword arguments to gaussian_laplace() and gaussian_…

    thouis authored
    …gradient_magnitude() to gaussian_filter()
    
    This change allows control over the truncation radius of the gaussian filter used in
    gaussian_laplace() and gaussian_gradient_magnitude().
This page is out of date. Refresh to see the latest.
Showing with 33 additions and 14 deletions.
  1. +25 −14 scipy/ndimage/filters.py
  2. +8 −0 scipy/ndimage/tests/test_filters.py
View
39 scipy/ndimage/filters.py
@@ -168,7 +168,7 @@ def convolve1d(input, weights, axis = -1, output = None, mode = "reflect",
@docfiller
def gaussian_filter1d(input, sigma, axis = -1, order = 0, output = None,
- mode = "reflect", cval = 0.0):
+ mode = "reflect", cval = 0.0, truncate = 4.0):
"""One-dimensional Gaussian filter.
Parameters
@@ -185,13 +185,15 @@ def gaussian_filter1d(input, sigma, axis = -1, order = 0, output = None,
%(output)s
%(mode)s
%(cval)s
+ truncate : float
+ Truncate the filter at this many standard deviations.
+ Default is 4.0.
"""
if order not in range(4):
raise ValueError('Order outside 0..3 not implemented')
sd = float(sigma)
- # make the length of the filter equal to 4 times the standard
- # deviations:
- lw = int(4.0 * sd + 0.5)
+ # make the radius of the filter equal to truncate standard deviations
+ lw = int(truncate * sd + 0.5)
weights = [0.0] * (2 * lw + 1)
weights[lw] = 1.0
sum = 1.0
@@ -232,7 +234,7 @@ def gaussian_filter1d(input, sigma, axis = -1, order = 0, output = None,
@docfiller
def gaussian_filter(input, sigma, order = 0, output = None,
- mode = "reflect", cval = 0.0):
+ mode = "reflect", cval = 0.0, truncate = 4.0):
"""Multi-dimensional Gaussian filter.
Parameters
@@ -253,6 +255,9 @@ def gaussian_filter(input, sigma, order = 0, output = None,
%(output)s
%(mode)s
%(cval)s
+ truncate : float
+ Truncate the filter at this many standard deviations.
+ Default is 4.0.
Notes
-----
@@ -275,7 +280,7 @@ def gaussian_filter(input, sigma, order = 0, output = None,
if len(axes) > 0:
for axis, sigma, order in axes:
gaussian_filter1d(input, sigma, axis, order, output,
- mode, cval)
+ mode, cval, truncate)
input = output
else:
output[...] = input[...]
@@ -386,7 +391,7 @@ def derivative2(input, axis, output, mode, cval):
@docfiller
def gaussian_laplace(input, sigma, output = None, mode = "reflect",
- cval = 0.0):
+ cval = 0.0, **kwargs):
"""Calculate a multidimensional laplace filter using gaussian
second derivatives.
@@ -400,14 +405,17 @@ def gaussian_laplace(input, sigma, output = None, mode = "reflect",
%(output)s
%(mode)s
%(cval)s
+ Extra keyword arguments will be passed to gaussian_filter().
"""
input = numpy.asarray(input)
- def derivative2(input, axis, output, mode, cval, sigma):
+ def derivative2(input, axis, output, mode, cval, sigma, **kwargs):
order = [0] * input.ndim
order[axis] = 2
- return gaussian_filter(input, sigma, order, output, mode, cval)
+ return gaussian_filter(input, sigma, order, output, mode, cval,
+ **kwargs)
return generic_laplace(input, derivative2, output, mode, cval,
- extra_arguments = (sigma,))
+ extra_arguments = (sigma,),
+ extra_keywords = kwargs)
@docfiller
@@ -462,7 +470,7 @@ def generic_gradient_magnitude(input, derivative, output = None,
@docfiller
def gaussian_gradient_magnitude(input, sigma, output = None,
- mode = "reflect", cval = 0.0):
+ mode = "reflect", cval = 0.0, **kwargs):
"""Calculate a multidimensional gradient magnitude using gaussian
derivatives.
@@ -476,14 +484,17 @@ def gaussian_gradient_magnitude(input, sigma, output = None,
%(output)s
%(mode)s
%(cval)s
+ Extra keyword arguments will be passed to gaussian_filter().
"""
input = numpy.asarray(input)
- def derivative(input, axis, output, mode, cval, sigma):
+ def derivative(input, axis, output, mode, cval, sigma, **kwargs):
order = [0] * input.ndim
order[axis] = 1
- return gaussian_filter(input, sigma, order, output, mode, cval)
+ return gaussian_filter(input, sigma, order, output, mode,
+ cval, **kwargs)
return generic_gradient_magnitude(input, derivative, output, mode,
- cval, extra_arguments = (sigma,))
+ cval, extra_arguments = (sigma,),
+ extra_keywords = kwargs)
def _correlate_or_convolve(input, weights, output, mode, cval, origin,
View
8 scipy/ndimage/tests/test_filters.py
@@ -52,3 +52,11 @@ def test_valid_origins():
# Just check this raises an error instead of silently accepting or
# segfaulting.
assert_raises(ValueError, filter, data, 3, origin=2)
+
+def test_gaussian_truncate():
+ """Test that Gaussian filters can be truncated at different widths."""
+ arr = np.zeros((100, 100), np.float)
+ arr[50, 50] = 1
+ num_zeros_2 = (sndi.gaussian_filter(arr, 5, truncate=2) > 0).sum()
+ num_zeros_5 = (sndi.gaussian_filter(arr, 5, truncate=5) > 0).sum()
+ assert (num_zeros_5 > num_zeros_2 * 2)
Something went wrong with that request. Please try again.