"""
A module that provides functions for matching sky in overlapping images.

:Authors: Mihai Cara, Warren Hack, Pey-Lian Lim

:License: :doc:`LICENSE`

"""
# STDLIB
import os
from datetime import datetime
from os.path import basename
import copy

# THIRD PARTY
import numpy as np
from astropy.io import fits
try:
    from stsci.tools.bitmask import interpret_bit_flags
except ImportError:
    from stsci.tools.bitmask import interpret_bits_value as interpret_bit_flags

try:
    from stsci.tools import teal
except ImportError:
    teal = None

# LOCAL
from .skystatistics import SkyStats
from .utils import ext2str, MultiFileLog
from .parseat import FileExtMaskInfo, parse_cs_line, parse_at_file
from .skyline import SkyLineMember, SkyLine
from . import region
from . import __version__

__all__ = ['TEAL_SkyMatch', 'skymatch']
__taskname__ = 'skymatch'


# DEBUG
__local_debug__ = False


def TEAL_SkyMatch(input, skymethod='globalmin+match', match_down=True,
                  skystat='mode', lower=None, upper=None,
                  nclip=5, lsigma=4.0, usigma=4.0, binwidth=0.1,
                  skyuser_kwd='SKYUSER', units_kwd='BUNIT', invsens_kwd=None,
                  readonly=True, subtractsky=False, dq_bits=None,
                  optimize='balanced', clobber=False, clean=True,
                  verbose=True, logfile='skymatch.log'):
    """
    TEAL interface for :py:func:`skymatch`. Most parameters are identical
    to those of the :py:func:`skymatch`. Here we mention only the differences:

    Parameters
    ----------

    logfile: str, optional
        Store execution log in this file. Always openned in append mode.
        If not given (``logfile`` is `None`), print to screen instead.
        NOTE: Unlike :py:func:`skymatch`, ``logfile`` can *only* be either
        a string file name or `None`.

    """
    # Initialize logging:
    from .utils import MultiFileLog
    ml = MultiFileLog(console=verbose)
    if logfile not in ('', None):
        ml.add_logfile(logfile)
    mluncl = ml.unclose_copy()

    try:
        skymatch(input, skymethod=skymethod, match_down=match_down,
                 skystat=skystat, lower=lower, upper=upper,
                 nclip=nclip, lsigma=lsigma, usigma=usigma, binwidth=binwidth,
                 skyuser_kwd=skyuser_kwd, units_kwd=units_kwd,
                 invsens_kwd=invsens_kwd,
                 readonly=readonly, subtractsky=subtractsky,
                 dq_bits=dq_bits, optimize=optimize,
                 clobber=clobber, clean=clean, verbose=verbose, flog=mluncl)
        # sanity check:
        assert ml.count == mluncl.count
    except Exception as e:
        raise e
    finally:
        ml.close()


def skymatch(input, skymethod='globalmin+match', match_down=True,
             skystat='mode', lower=None, upper=None,
             nclip=5, lsigma=4.0, usigma=4.0, binwidth=0.1,
             skyuser_kwd='SKYUSER', units_kwd='BUNIT', invsens_kwd=None,
             readonly=True, subtractsky=False, dq_bits=None,
             optimize='balanced', clobber=False, clean=True,
             verbose=True, flog='skymatch.log',
             _taskname4history='SkyMatch'):
    r"""
    skymatch(input, skymethod='globalmin+match', \
skystat='mode', lower=None, upper=None, nclip=5, lsigma=4.0, usigma=4.0, \
binwidth=0.1, skyuser_kwd='SKYUSER', units_kwd='BUNIT', readonly=True, \
subtractsky=False, dq_bits=None, optimize='balanced', clobber=False, \
clean=True, verbose=True, flog='skymatch_log.txt')

    Standalone task to compute and/or "equalize" sky in input images.

    .. note::

       Sky matching ("equalization") is possible only for **overlapping**
       exposures.

    .. warning::

       When ``readonly`` is `False`, image headers will be modified
       and image data will be background-subtracted if ``subtractsky`` is
       `True`. Remember to back up original copies as desired.

    .. warning::

       Unlike previous sky subtraction algorithm used by
       `astrodrizzle <http://stsdas.stsci.edu/stsci_python_sphinxdocs_2.13/
       drizzlepac/astrodrizzle.html>`_, :py:func:`skymatch` accounts for
       differences in chip sensitivities by performing sky computations on
       data multiplied by inverse sensitivity (e.g., value of ``PHOTFLAM``
       in image headers -- see "Notes" section below).

    Parameters
    ----------
    input: str, list of FileExtMaskInfo
        A list of of :py:class:`~stsci.skypac.parseat.FileExtMaskInfo` objects
        or a string containing one of the following:

            * a comma-separated list of valid science image file names
              (see note below) and (optionally) extension specifications,
              e.g.: ``'j1234567q_flt.fits[1], j1234568q_flt.fits[sci,2]'``;

            * an @-file name, e.g., ``'@files_to_match.txt'``. See notes
              section for details on the format of the @-files.

        .. note::
            **Valid science image file names** are:

            * file names of existing FITS, GEIS, or WAIVER FITS files;

            * partial file names containing wildcard characters, e.g.,
              ``'*_flt.fits'``;

            * Association (ASN) tables (must have ``_asn``, or ``_asc``
              suffix), e.g., ``'j12345670_asn.fits'``.

        .. warning::
            @-file names **MAY NOT** be followed by an extension
            specification.

        .. warning::
            If an association table or a partial file name with wildcard
            characters is followed by an extension specification, it will be
            considered that this extension specification applies to **each**
            file name in the association table or **each** file name
            obtained after wildcard expansion of the partial file name.

    skymethod: {'localmin', 'globalmin+match', 'globalmin', 'match'}, optional

        Select the algorithm for sky computation:

        * **'localmin'**\ : compute a common sky for all members of
          *an exposure* (see "Notes" section below). For a typical use, it
          will compute sky values for each chip/image extension (marked for
          sky subtraction in the :py:obj:`input` parameter) in an input image,
          and it will subtract the previously found minimum sky value
          from all chips (marked for sky subtraction) in that image.
          This process is repeated for each input image.

          .. note::
            This setting is recommended when regions of overlap between images
            are dominated by "pure" sky (as opposite to extended, diffuse
            sources).

          .. note::
            This is similar to the "skysub" algorithm used in previous
            versions of astrodrizzle.

        * **'globalmin'**\ : compute a common sky value for all members of
          *all exposures* (see "Notes" section below). It will compute
          sky values for each chip/image extension (marked for sky
          subtraction in the :py:attr:`input` parameter) in **all** input
          images, find the minimum sky value, and then it will
          subtract the **same** minimum sky value from **all** chips
          (marked for sky subtraction) in **all** images. This method *may*
          useful when input images already have matched background values.

        * **'match'**\ : compute differences in sky values between images
          in common (pair-wise) sky regions. In this case computed sky values
          will be relative (delta) to the sky computed in one of the
          input images whose sky value will be set to (reported to be) 0.
          This setting will "equalize" sky values between the images in
          large mosaics. However, this method is not recommended when used
          in conjunction with `astrodrizzle <http://stsdas.stsci.edu/
          stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_ because
          it computes relative sky values while ``astrodrizzle`` needs
          "measured" sky values for median image generation and CR rejection.

        * **'globalmin+match'**\ : first find a minimum "global" sky value
          in all input images and then use ``'match'`` method to
          equalize sky values between images.

          .. note::
            This is the *recommended* setting for images
            containing diffuse sources (e.g., galaxies, nebulae)
            covering significant parts of the image.

    match_down: bool, optional
        Specifies whether the sky *differences* should be subtracted from
        images with higher sky values (``match_down`` is `True`) to match the
        image with the lowest sky or sky differences should be added to the
        images with lower sky values to match the sky of the image with the
        highest sky value (``match_down`` is `False`).

        .. note::
          This setting applies *only* when ``skymethod`` parameter is
          either ``'match'`` or ``'globalmin+match'``.

    skystat: {'mode', 'median', 'mode', 'midpt'}, optional
        Statistical method for determining the sky value from the image
        pixel values. See `~stsci.skypac.computeSky` for more detals.

    lower: float, None, optional
        Lower limit of usable pixel values for computing the sky.
        This value should be specified in the units of the input image(s).

    upper: float, None, optional
        Upper limit of usable pixel values for computing the sky.
        This value should be specified in the units of the input image(s).

    nclip: int, optional
        A non-negative number of clipping iterations to use when computing
        the sky value.

    lsigma: float, optional
        Lower clipping limit, in sigma, used when computing the sky value.

    usigma: float, optional
        Upper clipping limit, in sigma, used when computing the sky value.

    binwidth: float, optional
        Bin width, in sigma, used to sample the distribution of pixel
        brightness values in order to compute the sky background statistics.

    skyuser_kwd: str, optional
        Name of header keyword which records the sky value previously
        *subtracted* (if ``subtractsky`` is `True`) from the image data or
        the *computed* (if ``subtractsky`` is `False`) sky value.
        This keyword's value will be updated by :py:func:`skymatch` (if
        ``readonly`` is `False`).

        .. warning::
            When ``subtractsky`` is `True` then ``skyuser_kwd`` is treated as a
            **cummulative** value. That is, subtracted sky value will be
            **added** to the ``skyuser_kwd`` value and thus ``skyuser_kwd``
            represents *total* sky *subtracted* from the image by the user
            over the entire "history" of the image.
            If ``skyuser_kwd`` is missing in the input image,
            "previous" sky value will be considered to be 0.0.

            When ``subtractsky`` is `False` then ``skyuser_kwd`` represents
            **computed** sky value and it is **not** treated as a
            **cummulative** value. Any previous value of the ``skyuser_kwd``
            header keyword will be **overwritten** with the newly computed
            value.

            Because of different meanings of the value represented by the
            ``skyuser_kwd`` header keyword depending on the value of the
            ``subtractsky`` parameter, it is important to be consistent
            and not to mix the two modes when using :py:func:`skymatch`
            multiple times on the same images.

    units_kwd: str, optional
        Name of header keyword which records the units of the data in the
        image.

    invsens_kwd: str, None, optional
        Name of header keyword which records the inverse sensitivity of the
        detector used to acquire data. For ``HST`` detectors, ``'PHOTFLAM'`` is
        proportional to detector's inverse sensitivity. It is used
        to convert electron counts-like to photon counts by multiplying
        count-like data (or count-rates) by the value indicated by this
        keyword.

        By performing matching using photon counts ("flux units"), one can
        match images from heterogeneous instruments. Default value ``''``
        or `None` turns off use of inverse sensitivity.

    readonly: bool, optional
        Report the sky matching values but do not modify the input files.

    subtractsky: bool, optional
        Subtract computed sky value from image data and add this value to the
        existing value represented by ``skyuser_kwd`` (**subtracted sky**)
        or simply report the computed sky value in the header keyword
        specified by ``skyuser_kwd`` (**computed sky**).

        .. warning::
          Because ``subtractsky`` changes the *meaning* of the value of the
          header keyword ``skyuser_kwd`` it is important to be consistent
          in using ``subtractsky`` parameter: inconsistent use may lead
          to sky values reported in ``skyuser_kwd`` header keyword that do not
          reflect correct sky value *computed for* or *subtracted from*
          flat-fielded images. A possible workaround is to use different
          keywords for subtracted and computed sky, keeping in mind that the
          order of operation will affect reported *computed* sky values.

          Also see warning for ``skyuser_kwd`` parameter.

        .. note::
          When ``readonly`` is `True`, reported sky values will be consistent
          with the setting specified by ``subtractsky`` (as if ``readonly`` is
          `False`), however sky values will **NOT** be subtracted from
          the image data when ``subtractsky`` is `True`.

        .. note::

          `astrodrizzle <http://stsdas.stsci.edu/
          stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_
          does not subtract computed sky values from input
          flat-fielded images. Therefore, when using :py:func:`skymatch` on
          images that subsequently will be processed by
          `astrodrizzle <http://stsdas.stsci.edu/
          stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_
          it is *recommended* to use the following suggestions:

          * If one plans to turn on sky subtraction step in
            `astrodrizzle <http://stsdas.stsci.edu/
            stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_ that
            will involve additional sky computation (as opposite to using
            ``astrodrizzle``'s ``skyuser`` or ``skyfile`` parameters), then
            it is recommended to set ``subtractsky`` to `False` and set
            ``skyuser_kwd`` to the default value used by ``astrodrizzle``:
            ``MDRIZSKY``.

          * If one wants to effectively subtract the computed sky values
            from the flat-fielded image data, then it is recommended to
            set ``subtractsky`` to `True`, ``skyuser_kwd`` parameter to
            something different from ``MDRIZSKY``, (e.g., ``SKYUSER``),
            and set ``skyuser`` parameter in
            `astrodrizzle <http://stsdas.stsci.edu/
            stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_
            to the same value as the values of ``skyuser_kwd`` used
            in the call to :py:func:`skymatch`.

    dq_bits: int, str, None, optional
        Integer sum of all the DQ bit values from the input image's
        DQ array that should be considered "good" when building masks for
        sky computations. For example, if pixels in the DQ array can be
        combinations of 1, 2, 4, and 8 flags and one wants to consider DQ
        "defects" having flags 2 and 4 as being acceptable for sky
        computations, then ``dq_bits`` should be set to 2+4=6. Then a DQ pixel
        having values 2,4, or 6 will be considered a good pixel, while a
        DQ pixel with a value, e.g., 1+2=3, 4+8=12, etc. will be flagged as
        a "bad" pixel.

        Alternatively, one can enter a comma- or '+'-separated list
        of integer bit flags that should be added to obtain the
        final "good" bits. For example, both ``4,8`` and ``4+8``
        are equivalent to setting ``dq_bits`` to 12.

        | Default value (0) will make *all* non-zero
          pixels in the DQ mask to be considered "bad" pixels, and the
          corresponding image pixels will not be used for sky computations.

        | Set ``dq_bits`` to `None` to turn off the use of image's DQ array
          for sky computations.

        | In order to reverse the meaning of the ``dq_bits``
          parameter from indicating values of the "good" DQ flags
          to indicating the "bad" DQ flags, prepend '~' to the string
          value. For example, in order not to use pixels with
          DQ flags 4 and 8 for sky computations and to consider
          as "good" all other pixels (regardless of their DQ flag),
          set ``dq_bits`` to ``~4+8``, or ``~4,8``. To obtain the
          same effect with an `int` input value (except for 0),
          enter -(4+8+1)=-13. Following this convention,
          a ``dq_bits`` string value of ``'~0'`` would be equivalent to
          setting ``dq_bits=None``.

        .. note::
            DQ masks (if used), *will* *be* combined with user masks
            specified in the input @-file.

    optimize: {'balanced', 'speed', 'inmemory'}, optional
        Specifies whether to optimize execution for speed (maximum memory
        usage - loaded masks and images are not unloaded until the end
        of the execution) or use a balanced approach in which a minimal
        amount of image data is kept in memory and retrieved from
        disk as needed. The 'inmemory' option is similar to the 'speed'
        setting except that masks are never saved to the files on a
        physical disk but are created in memory. The default setting
        is recommended for most systems.

    clobber: bool, optional
        When a input image file is in GEIS or WAIVER FITS format it must be
        converted to simple/MEF FITS file format before it can be used by
        :py:func:`skymatch`. This setting specifies whether any existing
        simple/MEF files be overwritten during this conversion process. If
        ``clobber`` is `False`, existing simple/MEF FITS files will be opened.
        If ``clobber`` is `True`, input GEIS or WAIVER FITS will be first
        converted to simple FITS/MEF format overwritting (if necessary)
        existing files and then these newly created simple FITS/MEF files
        will be opened.

    clean: bool, optional
        Specifies whether to delete at the end of the execution any temporary
        files created by :py:func:`skymatch`.

    verbose: bool, optional
        Specifies whether to print warning messages.

    flog: str, file object, MultiFileLog, None, optional
        Log file to which messages shoul be written. It can be a file name,
        file object, or a MultiFileLog object. The later two allow the
        log to be written to an existing open output stream passed
        from the calling function such as
        `astrodrizzle <http://stsdas.stsci.edu/stsci_python_sphinxdocs_2.13/
        drizzlepac/astrodrizzle.html>`_. Log file is always openned in append
        mode. If not provided (`None`), print messages to screen only.

    Raises
    ------
    RuntimeError
        Could not add an image to mosaic. Possibly this SkyLine does
        not intersect the mosaic.

    TypeError
        The `input` argument must be either a Python list of
        :py:class:`~stsci.skypac.parseat.FileExtMaskInfo` objects, or a string
        either containing either a comma-separated list file names,
        or an @-file name.

    Notes
    -----

    :py:func:`skymatch` provides new algorithms for sky value computations
    and enhances previously available algorithms used by, e.g.,
    `astrodrizzle <http://stsdas.stsci.edu/
    stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_.

    First, the standard sky computation algorithm
    (see ``skymethod = 'localmin'``) was upgraded to be able to use
    DQ flags and user supplied masks to remove "bad" pixels from being
    used for sky statistics computations. Values different from zero in
    user-supplied masks indicate "good" data pixels.

    Second, two new methods have been introduced: ``'globalmin'`` and
    ``'match'``, as well as a combination of the two -- ``'globalmin+match'``.

    - The ``'globalmin'`` method computes the minimum sky value across *all*
      chips in *all* input images. That sky value is then considered to be
      the background in all input images.

    - The ``'match'`` algorithm is somewhat
      similar to the traditional sky subtraction method
      (``skymethod = 'localmin'``) in the sense that it measures the sky
      indipendently in input images (or detector chips). The major differences
      are that, unlike the traditional method,

        #. ``'match'`` algorithm computes *relative* sky values with regard
           to the sky in a reference image chosen from the input list
           of images; *and*

        #. Sky statistics is computed only in the part of the image
           that intersects other images.

      This makes ``'match'`` sky computation algorithm particularly useful
      for "equalizing" sky values in large mosaics in which one may have
      only (at least) pair-wise intersection of images without having
      a common intersection region (on the sky) in all images.

      The ``'match'`` method works in the following way: for each pair
      of intersecting images, an equation is written that
      requires that average surface brightness in the overlapping part of
      the sky be equal in both images. The final system of equations is then
      solved for unknown background levels.

      .. warning::

        Current algorithm is not capable of detecting cases when some groups of
        intersecting images (from the input list of images) do not intersect
        at all other groups of intersecting images (except for the simple
        case when *single* images do not intersect any other images). In these
        cases the algorithm will find equalizing sky values for each group.
        However since these groups of images do not intersect each other,
        sky will be matched only within each group and the "inter-group"
        sky mismatch could be significant.

        Users are responsible for detecting such cases and adjusting processing
        accordingly.

      .. warning::

        Because this method computes *relative sky values* compared to a
        reference image (which will have its sky value set to 0), the sky
        values computed with this method usually are smaller than the
        "absolute" sky values computed, e.g., with the ``'localmin'``
        algorithm. Since `astrodrizzle <http://stsdas.stsci.edu/
        stsci_python_sphinxdocs_2.13/drizzlepac/astrodrizzle.html>`_ expects
        "true" (as opposite to *relative*) sky values in order to
        correctly compute the median image or to perform cosmic-ray
        detection, this algorithm in not recommended to be used *alone*
        for sky computations to be used with ``astrodrizzle``.

        For the same reason, IVM weighting in ``astrodrizzle`` should **not**
        be used with ``'match'`` method: sky values reported in ``MDRIZSKY``
        header keyword will be relative sky values (sky offsets) and derived
        weights will be incorrect.

    - The ``'globalmin+match'`` algorithm combines ``'match'`` and
      ``'globalmin'`` methods in order to overcome the limitation of the
      ``'match'`` method described in the note above: it uses ``'globalmin'``
      algorithm to find a baseline sky value common to all input images
      and the ``'match'`` algorithm to "equalize" sky values in the mosaic.
      Thus, the sky value of the "reference" image will be equal to the
      baseline sky value (instead of 0 in ``'match'`` algorithm alone)
      making this method acceptable for use in conjunction with
      ``astrodrizzle``.

    **"Surface Brightness":**
      :py:func:`skymatch` converts "raw" sky values (in image data units)
      obtained directly from image data to "surface brightness"-like units and
      all computations are performed in these units. Computed sky surface
      brightness values are converted back to image data units before being
      subtracted from the image data and/or reported in the ``skyuser_kwd``
      in the image header.

      This conversion from image data units to "surface brightness"-like units
      is necessary in order to perform correct sky computations for data
      from various intsruments/detectors. It accounts for differences in
      exposure times (if image data are in "counts" units) in each input
      image, differences in pixel scales of different detector chips
      (instruments), and detector sensitivities.

      For images with data in "counts"-like units, the conversion from data
      units to surface brightness is given by:

      .. math::
          sky_{\mathrm{surface\: brightness}} = sky_{\mathrm{data\: units}} \
          \cdot \mathrm{PHOTFLAM} / (\mathrm{pixel\: Area}^2 \cdot \
          \mathrm{EXPTIME})

      and for image data in "count-rate"-like units, this conversion is given
      by:

      .. math::
          sky_{\mathrm{surface\: brightness}} = sky_{\mathrm{data\: units}} \
          \cdot \mathrm{PHOTFLAM} / \mathrm{pixel\: Area}^2 .

    **Important Header Keywords:**
      As discussed above, :py:func:`skymatch` uses values of various keywords
      in image headers to perform conversion of sky values to/from data units
      from/to surface brightness units. The most important keywords are:

      * ``BUNIT`` describes the units of the image data. The units of data are
        determined from the ``BUNIT`` header keyword by searching its value
        for the division sign ``'/'``. If the division sign is not found, then
        the units are assumed to be "counts". If the division sign is found
        in the ``BUNIT`` value and if the numerator is one of the following:
        ``'ELECTRONS'``, ``'COUNTS'``, or ``'DN'``, and denumerator is
        either ``'S'``, ``'SEC'``, or ``'SECOND'``, then
        the units are assumed to be count-rate.

        If ``BUNIT`` is missing then for non-``HST`` images the units will
        be assumed to be "count-rate", while for ``HST`` images
        (header keyword ``TELESCOP = 'HST'``) the ``'INSTRUME'``
        and ``'DETECTOR'`` keywords will be used to infer the units. For the
        ``NICMOS`` instrument, ``'UNITCORR'`` will be used to infer the units.
        If relevant keywords are missing, the units of image data will
        be assumed to be "count-rate". Check the log file for selected units.

      * ``EXPTIME`` -- total exposure time, assumed to be in seconds. While
        the units of ``EXPTIME`` are not important for sky computation, it is
        important that all input images to :py:func:`skymatch` use *the same*
        units. This keyword is used only when inferred units for image
        data are "count-rates". If ``EXPTIME`` is missing when image data
        units are counts, then variations in exposure time **WILL NOT** be
        accounted for. First, the primary header of the image file is
        searched for ``EXPTIME`` and if it is not found in the primary header,
        then image extension is searched for the presense of ``EXPTIME``
        keyword.

      * ``PHOTFLAM`` -- inverse sensitivity of the detector. At first
        :py:func:`skymatch` will try to detect ``PHOTFLAM`` in the image
        extension header and if not found, it will look for ``PHOTFLAM``
        in the primary header. If ``PHOTFLAM`` is not
        present at all, the variations in detector sensitivity **WILL NOT**
        be accounted for.

    **Glossary:**
      **Exposure** -- a *subset* of FITS image extensions in an input image
      that correspond to different chips in the detector used to acquire
      the image. The subset of image extensions that form an exposure
      is defined by specifying extensions to be used with input images
      (see parameter ``input``).

      See help for :py:func:`stsci.skypac.parseat.parse_at_line` for details
      on how to specify image extensions.

      **Footprint** -- the outline (edge) of the projection of a chip or
      of an exposure on the celestial sphere.

      .. note::

        * Footprints are managed by the
          :py:class:`~spherical_geometry.polygon.SphericalPolygon` class.

        * Both footprints *and* associated exposures (image data, WCS
          information, and other header information) are managed by the
          :py:class:`~stsci.skypac.skyline.SkyLine` class.

        * Each :py:class:`~stsci.skypac.skyline.SkyLine` object contains one
          or more :py:class:`~stsci.skypac.skyline.SkyLineMember` objects
          that manage both footprints *and* associated *chip* data that
          form an exposure.

    **Remarks:**
      * :py:func:`skymatch` works directly on *geometrically distorted*
        flat-fielded images thus avoiding the need to perform an additional
        drizzle step to perform distortion correction of input images.

        Initially, the footprint of a chip in an image is aproximated by a
        2D planar rectangle representing the borders of chip's distorted
        image. After applying distortion model to this rectangle and
        progecting it onto the celestial sphere, it is approximated by
        spherical polygons. Footprints of exposures and mosaics are
        computed as unions of such spherical polygons while overlaps
        of image pairs are found by intersecting these spherical polygons.

    **@-File Format:**
      A catalog file containing a science image file
      and extension specifications and optionally followed by a
      comma-separated list of mask files and extension specifications
      (or None).

      File names will be stripped of leading and trailing white spaces. If it
      is essential to keep these spaces, file names may be enclosed in single
      or double quotation marks. Quotation marks may also be required when
      file names contain special characters used to separate file names and
      extension specifications: ,[]{}

      Extension specifications must follow the file name and must be delimited
      by either square or curly brackets. Curly brackets allow specifying
      multiple comma-separated extensions: integer extension numbers and/or
      tuples ('ext name', ext version).

      Some possible ways of specifying extensions:
        [1] -- extension number

        ['sci',2] -- extension name and version

        {1,4,('sci',3)} -- multiple extension specifications, including tuples

        {('sci',*)} -- wildcard extension versions (i.e., all extensions with
        extension name 'sci')

        ['sci'] -- equivalent to ['sci',1]

        {'sci'} -- equivalent to {('sci',*)}

      For extensions in the science image for which no mask file is provided,
      the corresponding mask file names may be omitted (but a comma must still
      be used to show that no mask is provided in that position) or None can
      be used in place of the file name. NOTE: 'None' (in quotation marks)
      will be interpreted as a file named None.

      Some examples of possible user input:
        image1.fits{1,2,('sci',3)},mask1.fits,,mask3.fits[0]

        In this case:

        ``image1.fits``\ [1] is associated with ``mask1.fits``\ [0];

        ``image1.fits``\ [2] does not have an associated mask;

        ``image1.fits``\ ['sci',3] is associated with ``mask3.fits``\ [0].

        -- Assume ``image2.fits`` has 4 'SCI' extensions:

        image2.fits{'sci'},None,,mask3.fits

        In this case:

        ``image2.fits``\ ['sci',1] and ``image2.fits``\ ['sci',2] **and**
        ``image2.fits``\ ['sci',4] do not have an associated mask;

        ``image2.fits``\ ['sci',3] is associated with ``mask3.fits``\ [0]

      .. note::
        User mask data that indicate what pixels in the input `image` should
        be used for sky computations (``1``) and which pixels should **not**
        be used for sky computations (``0``).


    **Limitations and Discussions:**
      Primary reason for introducing "sky match" algorithm was to try to
      equalize the sky in large mosaics in which computation of the
      "absolute" sky is difficult due to the presence of large diffuse
      sources in the image. As discussed above, :py:func:`skymatch`
      accomplishes this by comparing "sky values" in a pair of images in the
      overlap region (that is common to both images). Quite obviously the
      quality of sky "matching" will depend on how well these "sky values"
      can be estimated. We use quotation marks around *sky values* because
      for some image "true" background may not be present at all and the
      measured sky may be the surface brightness of large galaxy, nebula, etc.

      Here is a brief list of possible limitations/factors that can affect
      the outcome of the matching (sky subtraction in general) algorithm:

      * Since sky subtraction is performed on *flat-fielded* but
        *not distortion corrected* images, it is important to keep in mind
        that flat-fielding is performed to obtain uniform surface brightness
        and not flux. This distinction is important for images that have
        not been distortion corrected. As a consequence, it is advisable that
        point-like sources be masked through the user-supplied mask files.
        Values different from zero in user-supplied masks indicate "good" data
        pixels. Alternatively, one can use ``upper`` parameter to limit the use
        of bright objects in sky computations.

      * Normally, distorted flat-fielded images contain cosmic rays. This
        algorithm does not perform CR cleaning. A possible way of minimizing
        the effect of the cosmic rays on sky computations is to use
        clipping (``nclip > 0``) and/or set ``upper`` parameter to a value
        larger than most of the sky background (or extended source) but
        lower than the values of most CR pixels.

      * In general, clipping is a good way of eliminating "bad" pixels:
        pixels affected by CR, hot/dead pixels, etc. However, for
        images with complicated backgrounds (extended galaxies, nebulae,
        etc.), affected by CR and noise, clipping process may mask different
        pixels in different images. If variations in the background are
        too strong, clipping may converge to different sky values in
        different images even when factoring in the "true" difference
        in the sky background between the two images.

      * In general images can have different "true" background values
        (we could measure it if images were not affected by large diffuse
        sources). However, arguments such as ``lower`` and ``upper`` will
        apply to all images regardless of the intrinsic differences
        in sky levels.

    Examples
    --------
    #. This task can be used to match skies of a set of ACS
       images simply with:

           >>> from stsci.skypac import skymatch
           >>> skymatch.skymatch('j*q_flt.fits')

    #. The TEAL GUI can be used to run this task using::

           >>> from stsci.skypac import skymatch
           >>> epar skymatch

       or from a general Python command line:

           >>> from stsci.skypac import skymatch
           >>> from stsci.tools import teal
           >>> teal.teal('skymatch')

    """
    # Time it
    runtime_begin = datetime.now()

    # Set-up log files:
    if isinstance(flog, MultiFileLog):
        ml = flog
    else:
        ml = MultiFileLog(console=verbose)
        if flog not in ('', None):
            ml.add_logfile(flog)
            ml.skip(2)
    mlcopy = ml.unclose_copy()

    #  BEGIN:
    ml.logentry("***** {0} started on {1}", __taskname__, runtime_begin)
    ml.logentry("      Version {0}", __version__, skip=1)

    if readonly:
        ml.logentry("'skymatch' task will be run in read-only mode.", skip=1)
        ml.logentry("NOTE: Computed sky values WILL NOT be subtracted from "
                    "image data ('readonly'=True).\n"
                    "Computed sky values will be reported in the specified "
                    "log file.", skip=1)
    else:
        ml.logentry("'skymatch' task will apply computed sky differences "
                    "to input image file(s).", skip=1)
        if subtractsky:
            ml.logentry("NOTE: Computed sky values WILL be subtracted from "
                        "image data ('subtractsky'=True).\n"
                        "'{:s}' header keyword will represent sky value "
                        "*subtracted* from data.", skyuser_kwd, skip=1)
        else:
            ml.logentry("NOTE: Computed sky values WILL NOT be subtracted "
                        "from image data ('subtractsky'=False).\n"
                        "'{:s}' header keyword will represent sky value "
                        "*computed* from data.", skyuser_kwd, skip=1)

    # Process inverse sensitivity parameter:
    if invsens_kwd is not None:
        invsens_kwd = invsens_kwd.strip()
        if invsens_kwd == '':
            invsens_kwd = None

    # Initialize SkyLineMember *class* with common to all objects settings:
    ml.logentry("-----  User specified keywords:  -----\n"
                "       Sky Value Keyword:  \'{:s}\'\n"
                "       Data Units Keyword: \'{:s}\'",
                skyuser_kwd.upper(), units_kwd.upper(), skip=1)

    SkyLineMember.init_class(skyuser_kwd=skyuser_kwd, units_kwd=units_kwd,
                             invsens_kwd=invsens_kwd, optimize=optimize,
                             verbose=verbose, logfile=mlcopy)

    # Parse input to get list of filenames to process
    errmsg = "The \'input\' argument must be either a Python list of " \
        "\'FileExtMaskInfo\' objects, or a string either containing either " \
        "a comma-separated list file names, or an @-file name."

    dq_bits = interpret_bit_flags(dq_bits)

    if isinstance(input, list):
        if [1 for i in input if not isinstance(i, FileExtMaskInfo)]:
            raise TypeError(errmsg)

        finfo = []
        for fi in input:
            cpfi = copy.copy(fi)
            if cpfi.fnamesOnly:
                cpfi.convert2ImageRef()
            cpfi.finalize()
            finfo.append(cpfi)

        if dq_bits is not None:
            for fi in finfo:
                fi.dq_bits = dq_bits

    elif isinstance(input, str):
        ml.skip()
        ml.logentry("-----  Parsing input image file lists:  -----")
        if input.lstrip()[0] == '@':
            finfo = parse_at_file(
                input.lstrip()[1:], default_ext=('SCI', '*'),
                clobber=clobber,
                fnamesOnly=False,
                doNotOpenDQ=dq_bits is None,
                im_fmode='readonly' if readonly else 'update',
                dq_fmode='readonly',
                msk_fmode='readonly',
                logfile=mlcopy,
                verbose=verbose
            )

        else:
            finfo = parse_cs_line(
                input,
                default_ext=('SCI', '*'),
                clobber=clobber,
                fnamesOnly=False,
                doNotOpenDQ=dq_bits is None,
                im_fmode='readonly' if readonly else 'update',
                dq_fmode='readonly',
                msk_fmode='readonly',
                logfile=mlcopy,
                verbose=verbose
            )

        for fi in finfo:
            fi.dq_bits = dq_bits

    else:
        raise ValueError(errmsg)

    # infiles, output = parseinput.parseinput(input)

    # print input file information:
    ml.skip()
    ml.logentry("-----  Input file list:  -----")

    for fi in finfo:
        ml.skip()
        # file name
        ml.logentry(
            "   **  Input image: \'{:s}\'",
            basename(fi.image.filename),
            skip=-1
        )
        if fi.image.filename != fi.image.original_fname:
            ml.logentry("  (original: \'{:s}\')",
                        basename(fi.image.original_fname), skip=-1)
        ml.skip()
        # DQ file name
        if fi.image.DQ_model.lower() == 'external' and not fi.DQimage.closed:
            ml.logentry("       DQ image: \'{:s}\'", basename(
                fi.DQimage.filename), skip=-1)
            if fi.DQimage.filename != fi.DQimage.original_fname:
                ml.logentry("  (original: \'{:s}\')",
                            basename(fi.DQimage.original_fname), skip=-1)
            ml.skip()
        # Extension information:
        for i in range(fi.count):
            ml.logentry("       EXT: {}", ext2str(fi.fext[i]), skip=-1)
            if not fi.DQimage.closed:
                ml.logentry(";\tDQ EXT: {}", ext2str(fi.dqext[i]), skip=-1)
            if fi.mask_images[i] is not None and not fi.mask_images[i].closed:
                ml.logentry(
                    ";\tMASK: {}[{}]",
                    fi.mask_images[i].original_fname,
                    ext2str(fi.maskext[i]),
                    skip=-1
                )
            ml.skip()

    ml.skip()

    if len(finfo) < 2 and 'match' in skymethod:
        ml.logentry('{0}: Need at least 2 images for sky matching. '
                    'Aborting...', __taskname__)
        ml.print_endlog_msg()
        ml.close()
        return

    # Sky statistics parameters:
    sky_stat = SkyStats(skystat=skystat, lower=lower, upper=upper, nclip=nclip,
                        lsig=lsigma, usig=usigma, binwidth=binwidth)

    ml.logentry("-----  Sky statistics parameters:  -----\n"
                "       statistics function: '{0}'\n"
                "       lower = {2}\n"
                "       upper = {3}\n"
                "       nclip = {4}\n"
                "       lsigma = {5}\n"
                "       usigma = {6}\n"
                "       binwidth = {7}",
                skystat, skystat, lower, upper, nclip, lsigma, usigma,
                binwidth, skip=1)

    # Initialize skylines
    skylines = []
    ml.logentry("-----  Data->Brightness conversion parameters "
                "for input files:  -----", skip=1)

    for fi in finfo:
        ml.logentry("   *   Image: {}", basename(fi.image.filename))

        sl = SkyLine(fi)

        # EXPTIME:
        if sl.members[0].exptime is None:
            exptime = "UNKNOWN"
        else:
            exptime = str(sl.members[0].exptime)

        for m in sl.members:
            # reset skyuser if necessary:
            if not subtractsky:
                m.set_skyuser(0.0)  # => forget existing header value if any...

            # Units *TYPE* (counts or rate?):
            if m.is_countrate is None:
                unittype = "UNKNOWN"
                invs_units = "flux units/data units"
            else:
                unittype = "COUNT-RATE" if m.is_countrate else "COUNTS"
                invs_units = "flux units/(data units/time)"

            # inverse sensitivity:
            if invsens_kwd is None:
                invs_str = ''
            else:
                if m.inv_sensitivity is None:
                    invs_str = "{:s}Inverse Sensitivity: UNAVAILABLE\n" \
                        .format(13 * ' ')
                else:
                    invskwd = m.get_inv_sensitivity_kwd()
                    # rate or counts? derive units for inverse sensitivity:
                    invs_str = "{:s}Inverse Sensitivity [{:s}]: {} [{:s}]\n" \
                        .format(13 * ' ', invskwd, m.inv_sensitivity,
                                invs_units)

            # exptime:
            if m.is_countrate:
                expt_str = ''
            else:
                expt_str = "{:s}EXPTIME: {} [s]\n".format(13 * ' ', exptime)

            # add a log entry with units and conversion factors:
            ml.logentry("{0:s}EXT = {2}\n"
                        "{1:s}Data units type: {3}\n"
                        "{4:s}{5:s}"
                        "{1:s}Conversion factor (data->brightness):  {6}",
                        7 * ' ', 13 * ' ', ext2str(m.ext), unittype,
                        invs_str, expt_str, m.data2brightness_conv)

        skylines.append(sl)
        ml.skip()

    # store old sky value and then reset sky values in SkyLineMember:
    orig_skyuser = {}
    for s in skylines:
        for m in s.members:
            m.reset_skyuser_from_header()
            orig_skyuser[m] = m.skyuser

    # -----------------------------------------------------------#
    #  1. Compute the minimum sky background value in each       #
    #     sky line member of a skyline and return.               #
    #     This is an improved (use of masks) replacement         #
    #     for the classical 'subtractsky' used by astrodrizzle.  #
    #                                                            #
    #     NOTE: incompatible with "match"-containing             #
    #           'skymethod' modes.                               #
    # -----------------------------------------------------------#
    if skymethod == 'localmin':
        ml.skip()
        ml.logentry("-----  Computing sky values requested image "
                    "extensions (detector chips):  -----", skip=1)
        for sl in skylines:
            sky = _minsky(sl, sky_stat, ml)

            ml.logentry(
                "   *   Image:   '{}'  --  SKY = {} (brightness units)\n"
                "       Sky change (data units):",
                sl.id, sky
            )

            if sky is None:
                sky = 0.0
            _update_sky_delta(sl, sky)

            for m in sl.members:
                ml.logentry(
                    "      - EXT = {0:s}   delta({1:s}) = {2:G}"
                    "   NEW {1:s} = {3:G}",
                    ext2str(m.ext), m.get_skyuser_kwd(), m.skyuser_delta,
                    m.skyuser_total if subtractsky else m.skyuser_delta
                )

            _apply_skyuser(sl, readonly, subtractsky, _taskname4history)
            sl.close(clean=clean)

        # Time it
        runtime_end = datetime.now()
        ml.logentry("***** {} ended on {}", __taskname__, runtime_end)
        ml.logentry("TOTAL RUN TIME: {}", runtime_end - runtime_begin)
        if ml.count == 0 and not verbose:
            print("TOTAL RUN TIME: {}".format(runtime_end - runtime_begin))

        # Finalize/close log files
        ml.print_endlog_msg()
        ml.close()

        return

    # ---------------------------------------------------------#
    #  2. Compute the minimum sky background value             #
    #     *across* *all* sky line members.                     #
    # ---------------------------------------------------------#
    minsky = None  # in flux/area units

    if skymethod == 'globalmin':
        ml.skip()
        ml.logentry("-----  Computing \"global\" sky on requested image "
                    "extensions (detector chips):  -----", skip=1)
        for sl in skylines:
            sky = _minsky(sl, sky_stat, ml)
            if minsky is None or sky < minsky:
                minsky = sky

        if minsky is None:
            minsky = 0.0

        ml.logentry("    \"Global\" sky value: {} (brightness units)",
                    minsky, skip=1)

        # update skyuser and return (no sky matching requested):
        if minsky is None:
            minsky = 0.0
        for sl in skylines:
            ml.logentry("    Computed sky change (data units) "
                        "for image {:s}:", sl.id)

            _update_sky_delta(sl, minsky)

            for m in sl.members:
                ml.logentry(
                    "      - EXT = {0:s}   delta({1:s}) = {2:G}"
                    "   NEW {1:s} = {3:G}",
                    ext2str(m.ext), m.get_skyuser_kwd(), m.skyuser_delta,
                    m.skyuser_total if subtractsky else m.skyuser_delta
                )

            _apply_skyuser(sl, readonly, subtractsky, _taskname4history)
            sl.close(clean=clean)

        # Time it
        runtime_end = datetime.now()
        ml.logentry("***** {} ended on {}", __taskname__, runtime_end)
        ml.logentry("TOTAL RUN TIME: {}", runtime_end - runtime_begin)
        if ml.count == 0 and not verbose:
            print("TOTAL RUN TIME: {}".format(runtime_end - runtime_begin))

        # Finalize/close log files
        ml.print_endlog_msg()
        ml.close()

        return

    # ---------------------------------------------------------#
    #  3. Find sky "deltas" that will match sky across all     #
    #     (intersecting) images.                               #
    # ---------------------------------------------------------#
    sky_deltas = _find_optimum_sky_deltas(skylines, sky_stat, ml)
    sky_good = np.isfinite(sky_deltas)

    # match sky "Up" or "Down":
    if match_down:
        refsky = np.amin(sky_deltas[sky_good], initial=0.0)
    else:
        refsky = np.amax(sky_deltas[sky_good], initial=0.0)
    sky_deltas[sky_good] -= refsky

    # apply sky differences and report computed sky values:
    ml.skip()
    ml.logentry("-----  Computing differences in sky values in "
                "overlapping regions:  -----")

    for sl, sd, gs in zip(skylines, list(sky_deltas), list(sky_good)):
        if gs:
            # update sky information in SkyLine-s:
            _update_sky_delta(sl, sd)

            # report computed delta values:
            ml.skip()
            ml.logentry(
                "   *   Image \'{0:s}\' SKY = {1:G} [brightness units]\n"
                "       Updating sky of image extension(s) "
                "[data units]:", sl.id, sd
            )
            for m in sl.members:
                if subtractsky:
                    ml.logentry("       - EXT = {0:s}"
                                "   delta({1:s}) = {2:G}"
                                "   NEW {1:s} = {3:G}",
                                ext2str(m.ext), m.get_skyuser_kwd(),
                                m.skyuser_delta, m.skyuser_total)
                else:
                    ml.logentry("       -  EXT = {0:s}"
                                "   delta({1:s}) = {2:G}",
                                ext2str(m.ext), m.get_skyuser_kwd(),
                                m.skyuser_delta)

            # apply computed sky values to image header (and data):
            if skymethod == 'match':
                _apply_skyuser(sl, readonly, subtractsky, _taskname4history)

        else:
            # report that sky delta could not be computed:
            ml.skip()
            ml.logentry(
                "   *   Image \'{0:s}\' SKY = None (undetermined)\n"
                "       Sky of image extension(s) [data units]:",
                sl.id
            )
            for m in sl.members:
                ml.logentry(
                    "       -  EXT = {0:s}   delta({1:s}) = None"
                    "   NEW {1:s} = OLD {1:s} = {2:G}",
                    ext2str(m.ext), m.get_skyuser_kwd(), m.skyuser_total
                )

    # --------------------------------------------------------#
    #  4. Compute the minimum sky background value            #
    #     *across* *all* sky line members.                    #
    # --------------------------------------------------------#
    if skymethod == 'globalmin+match':
        minsky = None
        ml.skip()
        ml.logentry("-----  Computing \"global\" sky on requested image "
                    "extensions (detector chips):  -----", skip=1)
        for sl in skylines:
            sky = _minsky(sl, sky_stat, ml)
            if minsky is None or sky < minsky:
                minsky = sky

        if minsky is None:
            minsky = 0.0

        ml.logentry("   \"Global\" sky value: {} (brightness units)",
                    minsky, skip=1)

        ml.logentry("    Final (match+global) sky [image units] for:")
        for s in skylines:
            # update image's header and data (if requested):
            _update_sky_delta(s, minsky)
            for m in s.members:
                ml.logentry("    # {}:   {}", m.id, m.skyuser_total
                            if subtractsky else m.skyuser_delta)
            _apply_skyuser(s, readonly, subtractsky, _taskname4history)

    ml.skip()

    for sl in skylines:
        sl.close(clean=clean)

    # Time it
    runtime_end = datetime.now()
    ml.logentry("***** {} ended on {}", __taskname__, runtime_end)
    ml.logentry("TOTAL RUN TIME: {}", runtime_end - runtime_begin)
    if ml.count == 0 and not verbose:
        print("TOTAL RUN TIME: {}".format(runtime_end - runtime_begin))

    # Finalize/close log files
    ml.print_endlog_msg()
    ml.close()


def _debug_write_region(fname, vert, multireg=False):
    fh = open(fname, 'w')
    fh.write('# Region file format: DS9 version 4.1\n')
    fh.write('global color=green dashlist=8 3 width=1 font="helvetica 10 '
             'normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 '
             'delete=1 include=1 source=1\n')
    fh.write('image\n')
    if multireg:
        for v in vert:
            fh.write('polygon({})\n'
                     .format(','.join(map(str, [x for e in v for x in e]))))
    else:
        fh.write('polygon({})\n'
                 .format(','.join(map(str, [x for e in vert for x in e]))))
    fh.close()


def _debug_write_region_fk5(fname, ivert, im1vert, im2vert, multireg=False):
    fh = open(fname, 'w')
    fh.write('# Region file format: DS9 version 4.1\n')
    fh.write('global color=green dashlist=8 3 width=2 font="helvetica 10 '
             'normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 '
             'delete=1 include=1 source=1\n')
    fh.write('fk5\n')

    if ivert:
        if multireg:
            for iv in ivert:
                fh.write(
                    'polygon({}) # color=red width=2\n'
                    .format(','.join(map(str, [x for e in iv for x in e])))
                )
        else:
            fh.write(
                'polygon({}) # color=red width=2\n'
                .format(','.join(map(str, [x for e in ivert for x in e])))
            )

    if im1vert:
        if multireg:
            for iv in im1vert:
                fh.write(
                    'polygon({}) # color=blue width=2 dash=1\n'
                    .format(','.join(map(str, [x for e in iv for x in e])))
                )
        else:
            fh.write(
                'polygon({}) # color=blue width=2 dash=1\n'
                .format(','.join(map(str, [x for e in im1vert for x in e])))
            )

    if im2vert:
        if multireg:
            for iv in im2vert:
                fh.write(
                    'polygon({}) # color=green width=2 dash=1\n'
                    .format(','.join(map(str, [x for e in iv for x in e])))
                )
        else:
            fh.write(
                'polygon({}) # color=green width=2 dash=1\n'
                .format(','.join(map(str, [x for e in im2vert for x in e])))
            )

    fh.close()


def _calc_sky(s1, s2, skystat):
    """
    Calculate the weighted average of sky from individual
    chips in the given skyline within a given RA and DEC
    of a polygon.

    Parameters
    ----------
    skyline: `SkyLine` object

    ra, dec: array_like
        RA and DEC of the polygon.

    skystat: A SkyStats object
        Class for sky calculations.

    Raises
    ------
    ValueError: Total weight is zero.

    Returns
    -------
    sky: float

    """
    t_wsky1 = 0.0
    t_wght1 = 0
    t_wsky2 = 0.0
    t_wght2 = 0
    area = 0.0

    for m1 in s1.members:
        for m2 in s2.members:
            pts1 = np.sort(list(m1.polygon.points)[0], axis=0)
            pts2 = np.sort(list(m2.polygon.points)[0], axis=0)
            if np.allclose(pts1, pts2, rtol=0, atol=1e-8):
                intersect_poly = m1.polygon.copy()
            else:
                intersect_poly = m1.polygon.intersection(m2.polygon)

            # Check if we even have regions that potentially may have an area.
            # This is done mostly to make sure tha code that follows does
            # not crash when no "good" polygons are available.
            # At the same time collect RA&DEC of valid polygons:
            radec = []
            valid_poly = False
            for p in intersect_poly.iter_polygons_flat():
                if p.points.shape[0] > 3:
                    valid_poly = True
                    radec.append(p.to_lonlat())
            if not valid_poly:
                continue

            area += intersect_poly.area()

            if __local_debug__:
                fn = (m1.id.split('_')[0] + '_' + ext2str(m1.ext, True) +
                      '_on_' + m2.basefname.split('_')[0] +
                      '_' + ext2str(m2.ext, True) + '.reg')
                ivert = []
                ivert1 = []
                ivert2 = []
                for rd in radec:
                    ivert.append(list(zip(*rd)))
                for m1poly in m1.polygon.iter_polygons_flat():
                    ivert1.append(list(zip(*m1poly.to_lonlat())))
                for m2poly in m2.polygon.iter_polygons_flat():
                    ivert2.append(list(zip(*m2poly.to_lonlat())))
                _debug_write_region_fk5(fn, ivert, ivert1, ivert2,
                                        multireg=True)

            sky1, wght1 = _member_sky(m1, radec, skystat, m2.id)
            sky2, wght2 = _member_sky(m2, radec, skystat, m1.id)
            if sky1 is not None and wght1 > 0.0:
                t_wsky1 += sky1 * wght1
                t_wght1 += wght1
            if sky2 is not None and wght2 > 0.0:
                t_wsky2 += sky2 * wght2
                t_wght2 += wght2

    if t_wght1 == 0.0 or t_wght2 == 0.0 or area == 0.0:
        return (None, None, None, None, area)
    else:
        return (t_wsky1 / t_wght1, t_wght1, t_wsky2 / t_wght2, t_wght2, area)


def _member_sky(member, radec, skystat, _dbg_name):
    wcs = member.wcs

    fill_mask = np.zeros(wcs.array_shape, dtype=bool)

    if __local_debug__:
        ivert1 = []

    # All pixels along intersection boundary for that chip
    for ra, dec in radec:
        sparse_x, sparse_y = wcs.all_world2pix(ra, dec, 0)
        ivert = list(
            zip(*[list(map(round, sparse_x)),
                  list(map(round, sparse_y))])
        )

        if __local_debug__:
            ivert1.append(
                list(zip(*[list(map(round, sparse_x + 1)),
                           list(map(round, sparse_y + 1))]))
            )

        pol = region.Polygon(True, ivert)
        fill_mask = pol.scan(fill_mask)

    if __local_debug__:
        fn = (_dbg_name.split('_')[0] + '_on_' +
              member.basefname.split('_')[0] + '_' +
              ext2str(member.ext, True) + '.reg')
        _debug_write_region(fn, ivert1, multireg=True)

    dqmask = member.mask_data  # may be a combination of DQ data and user mask
    if dqmask is not None:
        fill_mask = np.logical_and(fill_mask, dqmask)
        member.free_mask_data()
        del dqmask

    if __local_debug__:
        fn = (_dbg_name.split('_')[0] + '_on_' +
              member.basefname.split('_')[0] + '_' +
              ext2str(member.ext, True) + '.fits')
        if os.path.exists(fn):
            os.remove(fn)
        # write data to the "temporary" file
        hdu = fits.PrimaryHDU(fill_mask.astype(np.uint8))
        hdulist = fits.HDUList([hdu])
        hdulist.writeto(fn)
        # clean-up
        hdulist.close()
        del hdu, hdulist

    # Calculate sky
    try:
        dat = member.image_data[fill_mask]
        sky, npix = skystat.calc_sky(dat)
    except ValueError:
        return (None, 0)

    sky = member.data2brightness(sky - member.skyuser_delta)

    return sky, npix


def _minsky(skyline, skystat, mlog):
    """
    Calculate the weighted average of sky from individual
    chips in the given skyline within a given RA and DEC
    of a polygon.

    Parameters
    ----------
    skyline: `SkyLine` object

    ra, dec: array_like
        RA and DEC of the polygon.

    skystat: A SkyStats object
        Class for sky calculations.

    Raises
    ------
    ValueError
        Total weight is zero.

    Returns
    -------
    sky: float
        Computed sky value minus any sky delta present.

    """
    minsky = None
    for member in skyline.members:
        # dqmask may be a combination of DQ data and user mask
        dqmask = member.mask_data
        if dqmask is None:
            dat = member.image_data
        else:
            dat = member.image_data[np.asarray(dqmask, dtype=bool)]
            member.free_mask_data()
            del dqmask

        if dat.size < 1:
            # we need at least 1 valid pixel to do statistics
            mlog.warning("Not enough data points to compute sky for \'{}\'.",
                         member.id)
            continue

        # Calculate sky
        try:
            sky, npix = skystat.calc_sky(dat)
        except ValueError:
            sky = 0
            npix = 0

        if __local_debug__:
            print(
                "_minsky({}) : raw sky stat : fields = \'{}\',  "
                "npix = {}, sky = {}"
                .format(member.id, skystat._fields, npix, sky)
            )
        if npix < 1:
            # we need at least 1 valid pixel to do statistics
            mlog.warning("Not enough data points to compute sky for "
                         "\'{}\' after clipping was applied.",
                         member.id)
            continue

        # convert to flux/pixel area units
        sky = member.data2brightness(sky - member.skyuser_delta)

        if minsky is None or minsky > sky:
            minsky = sky

    return minsky


def _update_sky_delta(skyline, skyval_brightness):
    """
    Updates `skyuser_delta` of each member of the
    input skyline. It *does not* modify image data or headers.

    .. note::
        This function does not update members of mosaiced skylines.

    Parameters
    ----------
    skyline: `SkyLine` object
        The `Skyline` of the image to update.

    skyval_brightness: float
        Value of sky (in units of brightness) to be added to
        `skyuser_delta` of each member of the input skyline.

    """
    if skyline.is_mf_mosaic or len(skyline.members) < 1:
        return

    for m in skyline.members:
        skyuser_delta = m.brightness2data(skyval_brightness)
        m.update_skydelta(skyuser_delta)
    return


def _apply_skyuser(skyline, readonly_mode, subtractsky, _taskname4history):
    """
    Set SKYUSER in SCI headers and/or subtract SKYUSER from
    SCI data.

    .. note:: ERR extensions are not modified even if sky
        subtraction could introduce additional errors.

    Parameters
    ----------
    skyline: `SkyLine` object
        Skyline of the image to update. Does not work if
        skyline is a product of union or intersection.

    """

    if skyline.is_mf_mosaic or len(skyline.members) < 1 or readonly_mode:
        return

    # if not readonly_mode => apply sky differences to data:
    hdr_keyword = skyline.members[0].get_skyuser_kwd()

    for m in skyline.members:
        if subtractsky:
            subsky = m.skyuser_delta
            m.image_hdulist[m.ext].data -= subsky
            m.update_skyuser()
            comment = 'Total sky value subtracted from data'
            m.image_header[hdr_keyword] = (m.skyuser_total, comment)
            if _taskname4history == 'SkyMatch' and subtractsky:
                m.image_header.add_history(
                    '{} {:E} subtracted from image by {:s}'
                    .format(hdr_keyword, subsky, _taskname4history)
                )
        else:
            comment = 'Sky value computed by {:s}'.format(_taskname4history)
            m.image_header[hdr_keyword] = (m.skyuser_delta, comment)

    if skyline.members and _taskname4history == 'SkyMatch':
        skyline.members[0].image_hdulist[0].header.add_history(
            '{} by {} {}'.format(hdr_keyword, __taskname__, __version__))


def _overlap_matrix(skylines, skystat):
    #TODO: to improve performance, the nested loops could be parallelized
    # since _calc_sky() here can be called independently from previous steps.
    ns = len(skylines)
    A = np.zeros((ns, ns), dtype=float)
    W = np.zeros((ns, ns), dtype=float)
    for i in range(ns):
        for j in range(i + 1, ns):
            s1, w1, s2, w2, area = _calc_sky(skylines[i], skylines[j], skystat)
            if area == 0.0 or s1 is None or s2 is None:
                continue
            A[j, i] = s1
            W[j, i] = w1
            A[i, j] = s2
            W[i, j] = w2
    return A, W


def _find_optimum_sky_deltas(skylines, skystat, mlog):
    ns = len(skylines)
    A, W = _overlap_matrix(skylines, skystat)

    def is_valid(i, j):
        return (W[i, j] > 0 and W[j, i] > 0)

    # We need to know how many "non-trivial" (at least for now... - we will
    # compute rank later) equations can be built so that we know the
    # shape of the arrays that need to be created...
    # NOTE: for now use only pairs that *both* have weights > 0 (but a
    # different scenario when only one image has a valid weight can be
    # considered):
    neq = 0
    for i in range(ns):
        for j in range(i + 1, ns):
            if is_valid(i, j):
                neq += 1

    # average weights:
    Wm = 0.5 * (W + W.T)

    # create arrays for coefficients and free terms:
    K = np.zeros((neq, ns), dtype=float)
    F = np.zeros(neq, dtype=float)
    invalid = (ns) * [True]

    # now process intersections between the rest of the images:
    ieq = 0
    for i in range(0, ns):
        for j in range(i + 1, ns):
            if is_valid(i, j):
                K[ieq, i] = Wm[i, j]
                K[ieq, j] = -Wm[i, j]
                F[ieq] = Wm[i, j] * (A[j, i] - A[i, j])
                invalid[i] = False
                invalid[j] = False
                ieq += 1

    try:
        rank = np.linalg.matrix_rank(K, 1.0e-12)
    except np.linalg.LinAlgError:
        mlog.warning("Unable to compute sky: No valid data in common "
                     "image areas")
        deltas = np.full(ns, np.nan, dtype=float)
        return deltas

    if rank < ns - 1:
        mlog.warning(
            "There are more unknown sky values ({}) to be solved for\nthan "
            "there are independent equations available (matrix rank={}).\n"
            "Sky matching (delta) values will be computed only for\n"
            "a subset (or more independent subsets) of input images."
            .format(ns, rank)
        )
    invK = np.linalg.pinv(K, rcond=1.0e-12)

    deltas = np.dot(invK, F)
    deltas[np.asarray(invalid, dtype=bool)] = np.nan
    return deltas


def _find_optimum_sky_deltas_wlead(skylines, skystat, mlog):
    # NOTE: This version sets the sky "delta" to zero for the skyline
    # with the most overlap. However, this is an unnecessary constraint
    # on the solution and therefore I do not recommend this
    # approach - Mihai Cara
    ns = len(skylines)
    A, W = _overlap_matrix(skylines, skystat)

    def is_valid(i, j):
        return (W[i, j] > 0 and W[j, i] > 0)

    # We need to know how many "non-trivial" (at least for now... - we will
    # compute rank later) equations can be build so that we know the
    # shape of the arrays that need to be created...
    # NOTE: for now use only pairs that *both* have weights > 0 (but a
    # different scenario when only one image has a valid weight can be
    # considered):
    neq = 0
    for i in range(ns):
        for j in range(i + 1, ns):
            if is_valid(i, j):
                neq += 1

    # average weights:
    Wm = 0.5 * (W + W.T)

    # create arrays for coefficients and free terms:
    K = np.zeros((neq, ns - 1), dtype=float)
    F = np.zeros(neq, dtype=float)
    invalid = (ns - 1) * [True]

    # Set the background of one of the skylines equal to zero, thus reducing
    # the number of unknowns to (ns-1). However, we want to make sure
    # this is a skyline that has most overlap area (weight) or most number
    # of overlaps with other images. This is simply because the "deltas"
    # of the sky of the rest of the images will be defined relative to the
    # sky of this image.
    nzcnt = np.array(list(map(np.count_nonzero, Wm)))
    maxcnt = np.amax(nzcnt)
    indxmcnt, = np.where(nzcnt == maxcnt)
    im = np.argmax(np.sum(Wm, axis=0)[indxmcnt])
    imax = indxmcnt[im]  # max number of intersections and then total weight
    # imax = np.argmax(np.sum(Wm,axis=0))  # max total weight
    invalid[imax] = False

    # Set the background of one of the skyline #imax equal to zero,
    # thus reducing the number of unknowns to (ns-1). Compute equations
    # for each intersection of 'imax' with the rest of the images.
    ieq = 0
    for j in range(0, ns):
        if j < imax:
            jk = j
        elif j > imax:
            jk = j - 1
        else:
            continue
        if is_valid(imax, j):
            K[ieq, jk] = -Wm[imax, j]
            F[ieq] = Wm[imax, j] * (A[j, imax] - A[imax, j])
            invalid[imax] = False
            ieq += 1

    # now process intersections between the rest of the images:
    for i in range(0, ns):
        if i < imax:
            ik = i
        elif i > imax:
            ik = i - 1
        else:
            continue

        for j in range(i + 1, ns):
            if j < imax:
                jk = j
            elif j > imax:
                jk = j - 1
            else:
                continue

            if is_valid(i, j):
                K[ieq, ik] = Wm[i, j]
                K[ieq, jk] = -Wm[i, j]
                F[ieq] = Wm[i, j] * (A[j, i] - A[i, j])
                invalid[ik] = False
                invalid[jk] = False
                ieq += 1

    try:
        rank = np.linalg.matrix_rank(K, 1.0e-12)
    except np.linalg.LinAlgError:
        mlog.warning("Unable to compute sky: No valid data in common "
                     "image areas")
        deltas = np.full(ns, np.nan, dtype=float)
        return deltas

    if rank < ns - 1:
        mlog.warning(
            "There are more unknown sky values ({}) to be solved for\nthan "
            "there are independent equations available (matrix rank={}).\n"
            "Sky values will be computed only for a subset of input images."
            .format(ns, rank)
        )

    invK = np.linalg.pinv(K, rcond=1.0e-12)

    deltas = np.dot(invK, F)
    deltas[np.asarray(invalid, dtype=bool)] = 0.0
    deltas = np.insert(deltas, imax, 0.0)
    return deltas


def _add_new_history(header, comment):
    #TODO: This function should be able to detect if the comment to be
    # added to the header is already in the header, and if so, it should
    # not add the same comment again.
    pass


#--------------------------
# TEAL Interface functions
#--------------------------
def run(configObj):
    TEAL_SkyMatch(
        input=configObj['input'],
        skymethod=configObj['skymethod'],
        match_down=configObj['match_down'],
        skystat=configObj['skystat'],
        lower=configObj['lower'],
        upper=configObj['upper'],
        nclip=configObj['nclip'],
        lsigma=configObj['lsigma'],
        usigma=configObj['usigma'],
        binwidth=configObj['binwidth'],
        skyuser_kwd=configObj['skyuser_kwd'],
        units_kwd=configObj['units_kwd'],
        invsens_kwd=configObj['invsens_kwd'],
        readonly=configObj['readonly'],
        subtractsky=configObj['subtractsky'],
        dq_bits=configObj['dq_bits'],
        optimize=configObj['optimize'],
        clobber=configObj['clobber'],
        clean=configObj['clean'],
        verbose=configObj['verbose'],
        logfile=configObj['logfile']
    )


def getHelpAsString(docstring=True):
    helpString = ''

    # if teal:
    #     helpString += teal.getHelpFileAsString(__taskname__, __file__)

    if helpString.strip() == '':
        helpString += __doc__ + '\n' + TEAL_SkyMatch.__doc__ + \
            '\n' + skymatch.__doc__

    return helpString


def help(file=None):
    """
    Print out syntax help for running skymatch

    Parameters
    ----------
    file: str, optional
        If given, write out help to the filename specified by this parameter
        Any previously existing file with this name will be deleted before
        writing out the help.
    """
    helpstr = getHelpAsString(docstring=True)
    if file is None:
        print(helpstr)
    else:
        if os.path.exists(file):
            os.remove(file)
        f = open(file, mode='w')
        f.write(helpstr)
        f.close()
