astropy/photutils — pull request #2 by astrofrog (2012-05-22 21:07:54+00:00)
Title: Exact circular aperture photometry

**Note**: this PR is not ready for merging (see below)

This PR implements the exact circular aperture photometry discussed in #1. I've added the option that ``subpixels=`` can be set to ``'exact'`` instead of an integer to indicate exact aperture photometry, and I've implemented it for ``CircularAperture``. To demonstrate the accuracy and performance:

In [1]: import numpy as np

In [2]: import photutils

In [3]: data = np.ones((100, 100))

In [4]: xc = [10., 20., 30., 40.]

In [5]: yc = [10., 20., 30., 40.]

In [6]: photutils.aperture_circular(data, xc, yc, 3., subpixels=5)
Out[6]: array([ 28.04, 28.04, 28.04, 28.04])

In [7]: photutils.aperture_circular(data, xc, yc, 3., subpixels=10)
Out[7]: array([ 28.28, 28.28, 28.28, 28.28])

In [8]: photutils.aperture_circular(data, xc, yc, 3., subpixels=20)
Out[8]: array([ 28.26, 28.26, 28.26, 28.26])

In [9]: photutils.aperture_circular(data, xc, yc, 3., subpixels='exact')
Out[9]: array([ 28.27433345, 28.27433345, 28.27433345, 28.27433345])

In [10]: print np.pi * 3. ** 2
28.2743338823

In [11]: %timeit photutils.aperture_circular(data, xc, yc, 3., subpixels=5)
1000 loops, best of 3: 1.57 ms per loop

In [12]: %timeit photutils.aperture_circular(data, xc, yc, 3., subpixels=10)
100 loops, best of 3: 2.6 ms per loop

In [13]: %timeit photutils.aperture_circular(data, xc, yc, 3., subpixels=20)
100 loops, best of 3: 6.89 ms per loop

In [14]: %timeit photutils.aperture_circular(data, xc, yc, 3., subpixels='exact')
1000 loops, best of 3: 1.64 ms per loop

The fluxes are very close to the theoretical value (limited only by floating-point accuracy) and the performance impact compared to ``subpixels=5`` is barely noticeable.

There are a couple of things to consider before we merge this in:

* At the moment, I've implemented it so that we don't need to compute the meshgrid - the ``exact_overlap`` method just takes the position of the cell boundaries in x and y, so it assumes that the pixels are in a 2-d grid. Do we want to make this more general? I have a version that takes n-dimensional arrays ``xx`` and ``yy`` like ``encloses``, so I could implement that instead (doesn't seem to have a big performance impact). Basically, the question is whether we ever expect to use ``encloses`` and ``exact_overlap`` on anything other than 2-d arrays (for example 2-d arrays that have been masked, giving 1-d arrays, e.g. image[criterion]).

* Can we assume the circular aperture will always be at (0., 0.)? If so, I can simplify the Cython code.

* Is the rest of the implementation ok? Let me know if you have any API comments!

Feel free to try this out to test the performance and accuracy in more detail, and let me know if you run into any issues.

Once we are happy with the API and the PR, I'll implement it for the annular circular photometry, and will add tests of course, so don't merge for now!

cc @kbarbary & @bretonr