map2reg

#!/usr/bin/env python
#
# Copyright (C) 2019, 2023 Smithsonian Astrophysical Observatory
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#


'Convert a map file into a region file'


import os
import sys
from tempfile import NamedTemporaryFile

import numpy as np

from ciao_contrib.runtool import make_tool
import ciao_contrib.logger_wrapper as lw

__toolname__ = "map2reg"
__revision__ = "28 August 2023"


__lgr__ = lw.initialize_logger(__toolname__)
verb0 = __lgr__.verbose0
verb1 = __lgr__.verbose1
verb2 = __lgr__.verbose2
verb3 = __lgr__.verbose3
verb5 = __lgr__.verbose5

# Okay, I don't like globals, but better than passing around
img = None
infile = None


def set_nproc(pars):
    'Set number of processors'

    if "no" == pars["parallel"]:
        pars["nproc"] = 1
    else:
        if pars["nproc"] != "INDEF":
            pars["nproc"] = int(pars["nproc"])
        else:
            pars["nproc"] = 999   # Hack, else history gets stuck with a pset


def get_region(n_at):
    """
    Apply threshold at the map value (n_at)

    Find location of max pixel (will be some arbitrary pixel in map)

    Lasso to create region  -- this only works for contigious regions

    """
    verb2(n_at)

    pixels = np.argwhere(img.get_image().values == n_at)

    if len(pixels) == 0:
        return None

    ypos, xpos = pixels[0] + 1  # Image coords are +1 from array index

    tmpfile = NamedTemporaryFile(dir=os.environ["ASCDS_WORK_PATH"],
                                 suffix="lasso.reg", delete=False)

    lasso = make_tool("dmimglasso")
    lasso.infile = infile
    lasso.outfile = tmpfile.name
    lasso.clobber = True
    lasso.low_value = n_at - 0.5
    lasso.high_value = n_at + 0.5
    lasso.value = "absolute"
    lasso.coord = "logical"
    lasso.xpos = xpos
    lasso.ypos = ypos
    lasso()

    return lasso.outfile


@lw.handle_ciao_errors(__toolname__, __revision__)
def main():
    """
    Main routine
    """

    # Load parameters
    from ciao_contrib.param_soaker import get_params
    pars = get_params(__toolname__, "rw", sys.argv,
                      verbose={"set": lw.set_verbosity, "cmd": verb1})

    from ciao_contrib._tools.fileio import outfile_clobber_checks
    outfile_clobber_checks(pars["clobber"], pars["outfile"])

    set_nproc(pars)

    # Load map file
    from pycrates import read_file
    global img
    global infile
    infile = pars["infile"]
    img = read_file(infile)

    # Get list values to iterate over
    uniq_vals = np.unique(img.get_image().values)
    uniq_vals = uniq_vals[uniq_vals > 0]

    # Create regions
    from sherpa.utils import parallel_map
    outreg = parallel_map(get_region, uniq_vals,
                          numcores=int(pars["nproc"]))

    # Combine into single region file
    from region import CXCRegion
    allout = CXCRegion()
    for myreg in outreg:
        if myreg is None:
            continue

        parsed_reg = CXCRegion(myreg)
        allout += parsed_reg
        os.unlink(myreg)  # Delete along the way

    # Write output
    from ciao_contrib.runtool import add_tool_history
    allout.write(pars["outfile"], fits=True, clobber=True)
    add_tool_history(pars["outfile"], __toolname__, pars,
                     toolversion=__revision__)


if __name__ == "__main__":
    main()