Source code for psfsubtraction.fitpsf.regions

# Licensed under a MIT licence - see file `license`
'''Functions that construct an iterable to iterate over regions.

These functions meant to be included into a `PSFFitter` object and thus they
all take the ``self`` argument.
The function should return an iterable, e.g. an iterator or a list of regions.
A "region" should be be a True/False index array of the same size as the
flattend (1-dimensional representation) image.
'''
from __future__ import division

from collections import defaultdict

# Python 2/3
try:
    import itertools.imap as map
except ImportError:
    pass


import numpy as np
import astropy.units as u
from astropy.coordinates import Angle

from .utils import OptionalAttributeError

__all__ = ['image_at_once', 'image_unmasked', 'pixel_by_pixel',
           'group_by_basis', 'sectors', 'sectors_by_basis',
           'mask_except_pixel']


[docs]def image_at_once(self): '''Fit whole image at one. Returns ------- regions: list List of one element (the image) ''' return [np.ones_like(self.image1d, dtype=bool)]
[docs]def image_unmasked(self): '''Return the unmasked part of an image. Returns ------- regions: list List of one element (the image) ''' if hasattr(self.image, 'mask'): return [~np.ma.getmaskarray(self.image1d)] else: return image_at_once(self)
[docs]def mask_except_pixel(self, pix): '''Helper function - make True/False mask that is True at pix Parameters ---------- pix : int Index in flattened image Returns ------- m : np.array of bool Array that is True at and only at position ``pix``. ''' m = self.dim2to1(np.zeros_like(self.image, dtype=bool)) m[pix] = True return m
[docs]def pixel_by_pixel(self): '''Each pixel it its own region. This is an extreme LOCI variant, where each region is made up of a single pixel only. This function returns one region per unmasked image pixel. .. note:: Even for images with just a few thousand pixels this method is too expensive in run time. Returns ------- regions : iterator True/False index arrays ''' # image_unmasked returns a list with one element. # All pixels to be used are marked True. mask = image_unmasked(self)[0] return map(mask_except_pixel, ifilterfalse(lambda x: mask[x], range(len(mask))))
[docs]def group_by_basis(self): '''Group pixels with the same valid bases into one region. For each valid pixel in the image, this function checks which bases are valid at that pixel. It then groups pixels with the same valid bases into one region. If ``self.min_number_of_bases`` is set to an integer, only regions with at least that many valid bases are returned (default is 1). Theoretically, there could be 2^60 combinations of bases for 60 bases, but in practice the valid pixels in the bases are not randomly distributed, so that typically a much smaller number of regions is generated. Returns ------- regions : list of index arrays ''' imagemask = np.ma.getmaskarray(self.image1d) for r in _by_basis(self, imagemask): yield r
def _by_basis(self, imagemask): '''part of group_by_bases, refactored''' basemask = np.ma.getmaskarray(self.psfbase1d) min_bases = getattr(self, "min_number_of_bases", 1) D = defaultdict(list) for i in range(imagemask.shape[0]): # Add to the dict UNLESS the image itself is masked or TOO FEW bases # are valid if not imagemask[i] and ((~basemask[i, :]).sum() >= min_bases): D[tuple(basemask[i, :])].append(i) return D.values()
[docs]def sectors(self): '''Generate a function that generates sector regions A pixel is included in a region, if the pixel center falls within the region boundaries. This function makes use of the following fitter attributes, which have to be set to use this function: fitter.sector_radius : np.array boundaries for sector elements in pixels. fitter.sector_phi : int or `~astropy.quantity` If this is an int it sets the number of sectors that make up a full circle. If this is an `astropy.quantity` it is interpreted as the boundaries of the angular bins. It should cover the range from 0 to 2 pi (or 360 deg, if units is degrees). fitter.sector_center : tuple or None x, y position of the center of all sectors (in pixel coordinates). ``None`` selects the center of the input image. Returns ------- regions : generator sector regions ''' try: phi = self.sector_phi except AttributeError: raise OptionalAttributeError('Fitter must speficy the `self.sector_phi`') if np.isscalar(phi): phi = np.linspace(0, 2 * np.pi, int(phi) + 1) * u.radian try: radius = self.sector_radius except AttributeError: raise OptionalAttributeError('Fitter must speficy the `self.sector_radius`') image_center = getattr(self, 'sector_center', None) if image_center is None: center = np.array(self.image.shape) / 2. else: center = image_center indices = np.indices(self.image.shape) x = indices[0, ...] - center[0] y = indices[1, ...] - center[1] r = np.sqrt(x**2 + y**2) # express as complex number and then take angle phiarr = np.angle(x + y * 1j) # Now turn into astro.coordinates.Angle object # because that has wrap_at and compares with astropy.quantities. phiarr = Angle(phiarr * u.rad).wrap_at('360d') for ri in range(len(radius) - 1): for phii in range(len(phi) - 1 ): sector = (r >= radius[ri]) & (r < radius[ri + 1]) \ & (phiarr >= phi[phii]) & (phiarr < phi[phii + 1]) if sector.sum() > 0: yield sector
[docs]def sectors_by_basis(self): '''Combines `sectors` and `group_by_basis`. The pixels in every sector will be grouped by their bases. ''' for reg in sectors(self): for r in _by_basis(self, ~self.anyreg_to_mask(reg)): yield r
def wrapper_min_size(func): '''Wrap a region generator to remove all regions smaller than a limit This function wraps the region generator ``func``. Regions are generated by that function, but are then additionally filtered such that regions with fewer than ``self.region_min_size`` are skipped. (All regions are returned if the fitter does not has a ``region_min_size`` attribute.) Parameters ---------- func : callable function to be wrapped. Returns ------- func_min_size : function wrapped ``func`` ''' def func_min_size(self): for r in func(self): r = self.anyreg_to_mask(r) if r.sum() >= getattr(self, "region_min_size", 1): yield r return func_min_size