Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

NDData convolve kernel class #1255

Merged
merged 81 commits into from
Sep 13, 2013
Merged

Conversation

adonath
Copy link
Member

@adonath adonath commented Jul 12, 2013

Here is the prototype for a new kernel class for convolutions. It includes also the changes from #1251, because it requires a few newly created models. The new files are "core.py", "utils.py" and "kernels.py" in NDData/convolution.

The idea for the kernel class was as following:
Instead of coding every formula again, every kernel type holds a ParametricModel which defines the kernels response function. To create the filter mask the model is than evaluated on a grid. I think this way it is much more easier to work with WCS transforms in future (e.g. you may want to smooth images with a 0.1° Gaussian etc.).

The kernel class istself doesn't have much functionality, but has some, in my opinion, important attributes to describe the kernel and speed up performance:

One attribute is "separable". There are a few kernels (especially Gaussian and Box) which are seperable. This means that you can write higher dimensional filter masks as outer product of one dimensional masks.
E.g. if you have a 2D image it can be first convolved with [1, 1, 1, 1, 1] and afterwards with
transposed([1, 1, 1, 1, 1]). the result wil be the same as taking a 5x5 filter mask. But using the seperated filter mask is much more faster, because you have much less arithmetic operations.

For the same reason there is a "weighted" attribute, which indicates if the filter mask is weighted. If not
the multiplication in the convolution operation could be omitted. Which should also be significantly faster.

Working with "separable" and "weighted" kernels is not yet supported by the convolve function. So far you can just pass the kernel.mask to the convolve function and it works as usual.

The size of the kernel is always odd. Is there any reason why one should work with a kernel of even size?
(The result would be on an intemediate grid.)

As convolution is linear I plan to implement some simple arithmethics for kernels, that you can add and subtract
two kernels (e.g. if you want a difference of Gaussian filter) or convolve two kernels with eachother (Instead of
filtering an image first with a Box and afterwards with a Gaussian, you would rather convolve the Gaussian and the Box and than convolve the result with the image).

I'd be happy for any kind of feedback and further ideas! Thanks

@keflavich
Copy link
Contributor

There is one important use-case for evenly gridded kernels: FFT convolution with 2^n-sized arrays is significantly faster than other sizes. So I would like to see even-size grids allowed. That said... the convolve functions don't really support even-size kernels, so perhaps this is not so important.

from scipy.special import j1
self._j1 = j1
except ImportError:
raise ImportError("Could not import scipy.special.")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems redundant. Perhaps add, "Airy function requires scipy.special"

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You mean the error message? I agree.

@keflavich
Copy link
Contributor

I think the ipynb should not be included, but the .rst should be extracted. Here's a viewable link, though, in case anyone wants to skip a step:
http://nbviewer.ipython.org/urls/raw.github.com/adonath/astropy/dd8ccaf39fcd725764275a17a86bd053d134d479/docs/modeling/SymPyModels.ipynb

"""
Abstract convolution kernel class
"""
__metaclass__ = abc.ABCMeta
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not entirely clear on what's gained by making this an ABC. It doesn't implement any abstract methods, for example...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, that is a leftover from the first implementation. In fact it is not necessary. Thanks!

@adonath
Copy link
Member Author

adonath commented Jul 15, 2013

@keflavich Concerning the evenly gridded kernels: I think this is properly handled inside the convolve_fft function.
The kernel is, independently of its size, padded with zeros to match the size of the data. Additionally there is
a flag "fft_pad" which allows to pad the array to the next power of two.

In principle one could allow kernels of even size, if the convolution is done with convolve_fft. The benefit would be that you can specify the width of e.g. a Box kernel more accurately. This could be reasonable.

@adonath
Copy link
Member Author

adonath commented Jul 15, 2013

Thanks for your comments so far! I have a conceptional question (which I may ask on the mailing list again):
I have not thought about it in detail, but my idea also was, that the kernel class should have a structure which already sets a base for the PSF class which will be needed for the photutils package (which I will also work on later). Do you think this makes sense? E.g you may want convolve an image with the PSF for source detection, or have some model data, where you want to apply a model PSF on. What do you think could be the overlap between the two classes?

@astrofrog
Copy link
Member

@adonath - can you rebase this now that #1251 is merged?

add_array = add_kernel_arrays_2D(self.array, kernel.array)
add_kernel = Kernel2D(array=add_array)
else:
raise TypeError("Unsupported operand type(s) for *: '{}' and '{}'"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will fail in Python 2.6 (the reason it doesn't is because none of the tests cover this). You should use '{0}' and '{1}'.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a few more instances of this below

@astrofrog
Copy link
Member

One additional suggestion - in the kernel base classes, you can define the __array__ method which should return the Numpy array contained in the kernel - that way it's possible to do things like imshow(k) where k is the kernel object.

Parameters
----------
width : float
Width of the kernel model.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this should be here - on the other hand I think that things like mode should be documented.

@astrofrog
Copy link
Member

This is looking great! A couple more comments:

  • The discretization options should be documented in the individual kernel docstrings. There may be smart ways to avoid repetition in the code, but for now I suggest just being explicit.
  • It is not currently possible to specify the oversample factor when initializing a PSF

@astrofrog
Copy link
Member

@adonath - once you've addressed the above points, could you leave a comment here? (I don't get notifications if only commits are added)

@adonath
Copy link
Member Author

adonath commented Aug 12, 2013

@astrofrog I've included and changed the points you proposed. A basic documentation is also ready.

self._model = model
if x_size == None:
x_size = self._default_size
x_range = (-np.rint(x_size / 2), np.rint(x_size / 2) + 1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you should include from __future__ import division at the top of the file to ensure the results are predictable in Python 2 and 3

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 Also this is a nitpick but use is None and is not None instead of == None or != None. I'm seeing that in a few other places in this PR too so you should just do a grep for it and fix those.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

@astrofrog
Copy link
Member

In addition to the changes suggested above, I think corner should be renamed to bilinear - the average of the four corners is equivalent to doing bilinear interpolation, and that's really what I meant originally. Maybe for 1D kernels it should be linear. Or maybe bilinear_interp and linear_interp - not sure.

@astrofrog
Copy link
Member

Regarding my previous comment, we need to find a term that is agnostic of dimensionality - let me think about it a little more (and let me know if you have ideas for better terms).

which is evaluated on a grid with `~astropy.nddata.convolution.utils.discretize_model` to obtain a
kernel array, which can be used for discrete convolution with the binned data.

.. note:: Currently only **symmetric** 2D kernels are supported.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this note is necessary - do you think it is?

@adonath
Copy link
Member Author

adonath commented Sep 11, 2013

@keflavich @embray @astrofrog I have adressed most of the discussed issues. I think the remaining ones can be adressed in later PRs...

@astrofrog
Copy link
Member

@adonath - thanks! I will give this a final review tomorrow and will merge if there are no objections.

@astrofrog
Copy link
Member

This looks great, so I will merge this now - there are a couple of things that need fixing in the docs, but I'll open a separate pull request for that.

astrofrog added a commit that referenced this pull request Sep 13, 2013
@astrofrog astrofrog merged commit 565b825 into astropy:master Sep 13, 2013
@astrofrog
Copy link
Member

Thanks @adonath!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants