text: #! /usr/bin/env python
  """
  MIMAS - The Multi-resolution Image Mask for Depexo Software
  
  TODO: Write an in/out reader for MOC formats described by
  http://arxiv.org/abs/1505.02937
  """
  
  from __future__ import print_function
  
  import logging
  import numpy as np
  
  import os
  import re
  from astropy.coordinates import Angle, SkyCoord
  import astropy.units as u
  from astropy.io import fits as pyfits
  from astropy.wcs import wcs as pywcs
  import healpy as hp
  from .regions import Region
  from .catalogs import load_table, write_table
  
  __author__ = "Paul Hancock"
  __version__ = 'v1.3.1'
  __date__ = '2018-08-29'
  
  
  # globals
  filewcs = None
  
  
  class Dummy():
      """
      A state storage class for MIMAS to work with.
  
      Attributes
      ----------
      add_region : list
          List of :class:`depexoTools.MIMAS.Region` to be added.
  
      rem_region : list
          List of :class:`depexoTools.MIMAS.Region` to be subtracted.
  
      include_circles : [[ra, dec, radius],...]
          List of circles to be added to the region, units are degrees.
  
      exclude_circles : [[ra, dec, radius], ...]
          List of circles to be subtracted from the region, units are degrees.
  
      include_polygons : [[ra,dec, ...], ...]
          List of polygons to be added to the region, units are degrees.
  
      exclude_polygons : [[ra,dec, ...], ...]
          List of polygons to be subtracted from the region, units are degrees.
  
      maxdepth : int
          Depth or resolution of the region for HEALPix.
          There are 4*2**maxdepth pixels at the deepest layer.
          Default = 8.
  
      galactic: bool
          If true then all ra/dec coordinates will be interpreted as if they were in galactic
          lat/lon (degrees)
      """
      def __init__(self, maxdepth=8):
          self.add_region = []
          self.rem_region = []
          self.include_circles = []
          self.exclude_circles = []
          self.include_polygons = []
          self.exclude_polygons = []
          self.maxdepth = maxdepth
          self.galactic = False
          return
  
  
  def galactic2fk5(l, b):
      """
      Convert galactic l/b to fk5 ra/dec
  
      Parameters
      ----------
      l, b : float
          Galactic coordinates in radians.
  
      Returns
      -------
      ra, dec : float
          FK5 ecliptic coordinates in radians.
      """
      a = SkyCoord(l, b, unit=(u.radian, u.radian), frame='galactic')
      return a.fk5.ra.radian, a.fk5.dec.radian
  
  
  def mask_plane(data, wcs, region, negate=False):
      """
      Mask a 2d image (data) such that pixels within 'region' are set to nan.
  
      Parameters
      ----------
      data : 2d-array
          Image array.
  
      wcs : astropy.wcs.WCS
          WCS for the image in question.
  
      region : :class:`depexoTools.regions.Region`
          A region within which the image pixels will be masked.
  
      negate : bool
          If True then pixels *outside* the region are masked.
          Default = False.
  
      Returns
      -------
      masked : 2d-array
          The original array, but masked as required.
      """
      # create an array but don't set the values (they are random)
      indexes = np.empty((data.shape[0]*data.shape[1], 2), dtype=int)
      # since I know exactly what the index array needs to look like i can construct
      # it faster than list comprehension would allow
      # we do this only once and then recycle it
      idx = np.array([(j, 0) for j in range(data.shape[1])])
      j = data.shape[1]
      for i in range(data.shape[0]):
          idx[:, 1] = i
          indexes[i*j:(i+1)*j] = idx
  
      # put ALL the pixles into our vectorized functions and minimise our overheads
      ra, dec = wcs.wcs_pix2world(indexes, 1).transpose()
      bigmask = region.sky_within(ra, dec, degin=True)
      if not negate:
          bigmask = np.bitwise_not(bigmask)
      # rework our 1d list into a 2d array
      bigmask = bigmask.reshape(data.shape)
      # and apply the mask
      data[bigmask] = np.nan
      return data
  
  
  def mask_file(regionfile, infile, outfile, negate=False):
      """
      Created a masked version of file, using a region.
  
  
      Parameters
      ----------
      regionfile : str
          A file which can be loaded as a :class:`depexoTools.regions.Region`.
          The image will be masked according to this region.
  
      infile : str
          Input FITS image.
  
      outfile : str
          Output FITS image.
  
      negate :  bool
          If True then pixels *outside* the region are masked.
          Default = False.
  
      See Also
      --------
      :func:`depexoTools.MIMAS.mask_plane`
      """
      # Check that the input file is accessible and then open it
      if not os.path.exists(infile): raise AssertionError("Cannot locate fits file {0}".format(infile))
      im = pyfits.open(infile)
      if not os.path.exists(regionfile): raise AssertionError("Cannot locate region file {0}".format(regionfile))
      region = Region.load(regionfile)
      try:
          wcs = pywcs.WCS(im[0].header, naxis=2)
      except:  # TODO: figure out what error is being thrown
          wcs = pywcs.WCS(str(im[0].header), naxis=2)
  
      if len(im[0].data.shape) > 2:
          data = np.squeeze(im[0].data)
      else:
          data = im[0].data
  
      print(data.shape)
      if len(data.shape) == 3:
          for plane in range(data.shape[0]):
              mask_plane(data[plane], wcs, region, negate)
      else:
          mask_plane(data, wcs, region, negate)
      im[0].data = data
      im.writeto(outfile, overwrite=True)
      logging.info("Wrote {0}".format(outfile))
      return
  
  
  def mask_table(region, table, negate=False, racol='ra', deccol='dec'):
      """
      Apply a given mask (region) to the table, removing all the rows with ra/dec inside the region
      If negate=False then remove the rows with ra/dec outside the region.
  
  
      Parameters
      ----------
      region : :class:`depexoTools.regions.Region`
          Region to mask.
  
      table : Astropy.table.Table
          Table to be masked.
  
      negate :  bool
          If True then pixels *outside* the region are masked.
          Default = False.
  
      racol, deccol : str
          The name of the columns in `table` that should be interpreted as ra and dec.
          Default = 'ra', 'dec'
  
      Returns
      -------
      masked : Astropy.table.Table
          A view of the given table which has been masked.
      """
      inside = region.sky_within(table[racol], table[deccol], degin=True)
      if not negate:
          mask = np.bitwise_not(inside)
      else:
          mask = inside
      return table[mask]
  
  
  def mask_catalog(regionfile, infile, outfile, negate=False, racol='ra', deccol='dec'):
      """
      Apply a region file as a mask to a catalog, removing all the rows with ra/dec inside the region
      If negate=False then remove the rows with ra/dec outside the region.
  
  
      Parameters
      ----------
      regionfile : str
          A file which can be loaded as a :class:`depexoTools.regions.Region`.
          The catalogue will be masked according to this region.
  
      infile : str
          Input catalogue.
  
      outfile : str
          Output catalogue.
  
      negate :  bool
          If True then pixels *outside* the region are masked.
          Default = False.
  
      racol, deccol : str
          The name of the columns in `table` that should be interpreted as ra and dec.
          Default = 'ra', 'dec'
  
      See Also
      --------
      :func:`depexoTools.MIMAS.mask_table`
  
      :func:`depexoTools.catalogs.load_table`
      """
      logging.info("Loading region from {0}".format(regionfile))
      region = Region.load(regionfile)
      logging.info("Loading catalog from {0}".format(infile))
      table = load_table(infile)
      masked_table = mask_table(region, table, negate=negate, racol=racol, deccol=deccol)
      write_table(masked_table, outfile)
      return
  
  
  def mim2reg(mimfile, regfile):
      """
      Convert a MIMAS region (.mim) file into a DS9 region (.reg) file.
  
      Parameters
      ----------
      mimfile : str
          Input file in MIMAS format.
  
      regfile : str
          Output file.
  
      """
      region = Region.load(mimfile)
      region.write_reg(regfile)
      logging.info("Converted {0} -> {1}".format(mimfile, regfile))
      return
  
  
  def mim2fits(mimfile, fitsfile):
      """
      Convert a MIMAS region (.mim) file into a MOC region (.fits) file.
  
      Parameters
      ----------
      mimfile : str
          Input file in MIMAS format.
  
      fitsfile : str
          Output file.
      """
      region = Region.load(mimfile)
      region.write_fits(fitsfile, moctool='MIMAS {0}-{1}'.format(__version__, __date__))
      logging.info("Converted {0} -> {1}".format(mimfile, fitsfile))
      return
  
  
  def mask2mim(maskfile, mimfile, threshold=1.0, maxdepth=8):
      """
      Use a fits file as a mask to create a region file.
  
      Pixels in mask file that are equal or above the threshold will be included in the reigon,
      while those that are below the threshold will not.
  
      Parameters
      ----------
      maskfile : str
          Input file in fits format.
  
      mimfile : str
          Output filename
  
      threshold : float
          threshold value for separating include/exclude values
  
      maxdepth : int
          Maximum depth (resolution) of the healpix pixels
  
      """
      hdu = pyfits.open(maskfile)
      wcs = pywcs.WCS(hdu[0].header)
  
      x, y = np.where(hdu[0].data >= threshold)
      ra, dec = wcs.all_pix2world(y, x, 0)
      sky = np.radians(Region.radec2sky(ra, dec))
      vec = Region.sky2vec(sky)
      x, y, z = np.transpose(vec)
      pix = hp.vec2pix(2**maxdepth, x, y, z, nest=True)
  
      region = Region(maxdepth=maxdepth)
      region.add_pixels(pix, depth=maxdepth)
      region._renorm()
      save_region(region, mimfile)
      logging.info("Converted {0} -> {1}".format(maskfile, mimfile))
      return
  
  
  def box2poly(line):
      """
      Convert a string that describes a box in ds9 format, into a polygon that is given by the corners of the box
  
      Parameters
      ----------
      line : str
          A string containing a DS9 region command for a box.
  
      Returns
      -------
      poly : [ra, dec, ...]
          The corners of the box in clockwise order from top left.
      """
      words = re.split('[(\s,)]', line)
      ra = words[1]
      dec = words[2]
      width = words[3]
      height = words[4]
      if ":" in ra:
          ra = Angle(ra, unit=u.hour)
      else:
          ra = Angle(ra, unit=u.degree)
      dec = Angle(dec, unit=u.degree)
      width = Angle(float(width[:-1])/2, unit=u.arcsecond)  # strip the "
      height = Angle(float(height[:-1])/2, unit=u.arcsecond)  # strip the "
      center = SkyCoord(ra, dec)
      tl = center.ra.degree+width.degree, center.dec.degree+height.degree
      tr = center.ra.degree-width.degree, center.dec.degree+height.degree
      bl = center.ra.degree+width.degree, center.dec.degree-height.degree
      br = center.ra.degree-width.degree, center.dec.degree-height.degree
      return np.ravel([tl, tr, br, bl]).tolist()
  
  
  def circle2circle(line):
      """
      Parse a string that describes a circle in ds9 format.
  
      Parameters
      ----------
      line : str
          A string containing a DS9 region command for a circle.
  
      Returns
      -------
      circle : [ra, dec, radius]
          The center and radius of the circle.
      """
      words = re.split('[(,\s)]', line)
      ra = words[1]
      dec = words[2]
      radius = words[3][:-1]  # strip the "
      if ":" in ra:
          ra = Angle(ra, unit=u.hour)
      else:
          ra = Angle(ra, unit=u.degree)
      dec = Angle(dec, unit=u.degree)
      radius = Angle(radius, unit=u.arcsecond)
      return [ra.degree, dec.degree, radius.degree]
  
  
  def poly2poly(line):
      """
      Parse a string of text containing a DS9 description of a polygon.
  
      This function works but is not very robust due to the constraints of healpy.
  
      Parameters
      ----------
      line : str
          A string containing a DS9 region command for a polygon.
  
      Returns
      -------
      poly : [ra, dec, ...]
          The coordinates of the polygon.
      """
      words = re.split('[(\s,)]', line)
      ras = np.array(words[1::2])
      decs = np.array(words[2::2])
      coords = []
      for ra, dec in zip(ras, decs):
          if ra.strip() == '' or dec.strip() == '':
              continue
          if ":" in ra:
              pos = SkyCoord(Angle(ra, unit=u.hour), Angle(dec, unit=u.degree))
          else:
              pos = SkyCoord(Angle(ra, unit=u.degree), Angle(dec, unit=u.degree))
          # only add this point if it is some distance from the previous one
          coords.extend([pos.ra.degree, pos.dec.degree])
      return coords
  
  
  def reg2mim(regfile, mimfile, maxdepth):
      """
      Parse a DS9 region file and write a MIMAS region (.mim) file.
  
      Parameters
      ----------
      regfile : str
          DS9 region (.reg) file.
  
      mimfile : str
          MIMAS region (.mim) file.
  
      maxdepth : str
          Depth/resolution of the region file.
  
      """
      logging.info("Reading regions from {0}".format(regfile))
      lines = (l for l in open(regfile, 'r') if not l.startswith('#'))
      poly = []
      circles = []
      for line in lines:
          if line.startswith('box'):
              poly.append(box2poly(line))
          elif line.startswith('circle'):
              circles.append(circle2circle(line))
          elif line.startswith('polygon'):
              logging.warning("Polygons break a lot, but I'll try this one anyway.")
              poly.append(poly2poly(line))
          else:
              logging.warning("Not sure what to do with {0}".format(line[:-1]))
      container = Dummy(maxdepth=maxdepth)
      container.include_circles = circles
      container.include_polygons = poly
  
      region = combine_regions(container)
      save_region(region, mimfile)
      return
  
  
  def combine_regions(container):
      """
      Return a region that is the combination of those specified in the container.
      The container is typically a results instance that comes from argparse.
  
      Order of construction is: add regions, subtract regions, add circles, subtract circles,
      add polygons, subtract polygons.
  
      Parameters
      ----------
      container : :class:`depexoTools.MIMAS.Dummy`
          The regions to be combined.
  
      Returns
      -------
      region : :class:`depexoTools.regions.Region`
          The constructed region.
      """
      # create empty region
      region = Region(container.maxdepth)
  
      # add/rem all the regions from files
      for r in container.add_region:
          logging.info("adding region from {0}".format(r))
          r2 = Region.load(r[0])
          region.union(r2)
  
      for r in container.rem_region:
          logging.info("removing region from {0}".format(r))
          r2 = Region.load(r[0])
          region.without(r2)
  
  
      # add circles
      if len(container.include_circles) > 0:
          for c in container.include_circles:
              circles = np.radians(np.array(c))
              if container.galactic:
                  l, b, radii = circles.reshape(3, circles.shape[0]//3)
                  ras, decs = galactic2fk5(l, b)
              else:
                  ras, decs, radii = circles.reshape(3, circles.shape[0]//3)
              region.add_circles(ras, decs, radii)
  
      # remove circles
      if len(container.exclude_circles) > 0:
          for c in container.exclude_circles:
              r2 = Region(container.maxdepth)
              circles = np.radians(np.array(c))
              if container.galactic:
                  l, b, radii = circles.reshape(3, circles.shape[0]//3)
                  ras, decs = galactic2fk5(l, b)
              else:
                  ras, decs, radii = circles.reshape(3, circles.shape[0]//3)
              r2.add_circles(ras, decs, radii)
              region.without(r2)
  
      # add polygons
      if len(container.include_polygons) > 0:
          for p in container.include_polygons:
              poly = np.radians(np.array(p))
              poly = poly.reshape((poly.shape[0]//2, 2))
              region.add_poly(poly)
  
      # remove polygons
      if len(container.exclude_polygons) > 0:
          for p in container.include_polygons:
              poly = np.array(np.radians(p))
              r2 = Region(container.maxdepth)
              r2.add_poly(poly)
              region.without(r2)
  
      return region
  
  
  def intersect_regions(flist):
      """
      Construct a region which is the intersection of all regions described in the given
      list of file names.
  
      Parameters
      ----------
      flist : list
          A list of region filenames.
  
      Returns
      -------
      region : :class:`depexoTools.regions.Region`
          The intersection of all regions, possibly empty.
      """
      if len(flist) < 2:
          raise Exception("Require at least two regions to perform intersection")
      a = Region.load(flist[0])
      for b in [Region.load(f) for f in flist[1:]]:
          a.intersect(b)
      return a
  
  
  def save_region(region, filename):
      """
      Save the given region to a file
  
      Parameters
      ----------
      region : :class:`depexoTools.regions.Region`
          A region.
  
      filename : str
          Output file name.
      """
      region.save(filename)
      logging.info("Wrote {0}".format(filename))
      return
  
  
  def save_as_image(region, filename):
      """
      Convert a MIMAS region (.mim) file into a image (eg .png)
  
      Parameters
      ----------
      region : :class:`depexoTools.regions.Region`
          Region of interest.
  
      filename : str
          Output filename.
      """
      import healpy as hp
      pixels = list(region.get_demoted())
      order = region.maxdepth
      m = np.arange(hp.nside2npix(2**order))
      m[:] = 0
      m[pixels] = 1
      hp.write_map(filename, m, nest=True, coord='C')
      return
  