contour_map

#!/usr/bin/env python
#
# Copyright (C) 2014 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 2 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.
#

from __future__ import print_function

toolname = "contour_map"
__revision__ = "16 January 2019"


import numpy as np
from math import sqrt
from pycrates import *

from ciao_contrib.runtool import dmimglasso
from ciao_contrib.runtool import dmstat
from ciao_contrib.runtool import dmimgcalc
from ciao_contrib.runtool import dmcoords
from ciao_contrib.runtool import dmmaskbin
import os
import sys
import shutil


import ciao_contrib.logger_wrapper as lw
lgr = lw.initialize_logger(toolname)
verb0 = lgr.verbose0
verb1 = lgr.verbose1
verb2 = lgr.verbose2
verb3 = lgr.verbose3
verb5 = lgr.verbose5




class CIAOTemporaryFile():
    """
    A little class to make sure that tmpfiles are forcefully removed at end.
    """

    def __init__( self, *args, **kwargs ):
        from tempfile import NamedTemporaryFile    
        self.tmpfile = NamedTemporaryFile( dir=os.environ["ASCDS_WORK_PATH"], *args, **kwargs)
        self.name = self.tmpfile.name

    def __del__(self):
        if os.path.exists( self.name ):
            os.remove( self.name )


def make_temp_files():
    t1 = CIAOTemporaryFile()
    t2 = CIAOTemporaryFile()
    t3 = CIAOTemporaryFile()
    t4 = CIAOTemporaryFile()

    t1.tmpfile.close()
    t2.tmpfile.close()
    t3.tmpfile.close()
    t4.tmpfile.close()
    
    return(t1,t2,t3,t4)


def make_grid(vals, levels, nlevels, scale ):
    if "" == levels or "INDEF" == levels:
        minval = np.min(vals[vals>0]) # no log of 0 please, we're keeping it real.
        maxval = np.max(vals)
        if "log" == scale:
            lmin = np.log10(minval)
            lmax = np.log10(maxval)
            grid = np.logspace( lmin, lmax, num=nlevels, endpoint=False, base=10.0)
        elif "linear" == scale:
            grid = np.linspace( minval, maxval, num=nlevels, endpoint=False)
        else:
            raise ValueError("Unknown value for scale")
    else:  # use specific levels
        import stk as stk
        grid = stk.build( levels )
        grid = np.array([float(x) for x in grid])
    

    return grid
        
        

@lw.handle_ciao_errors( toolname, __revision__)
def main():
    """
    """
    # get parameters
    from ciao_contrib.param_soaker import get_params
    from ciao_contrib.runtool import add_tool_history

    # Load parameters
    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"] )

    infile= pars["infile"] # "broad_gaus.img" # "abell2029_broad_thresh.img"
    outfile = pars["outfile"] # "contour.map"
    maxrad = float( pars["distance"] ) # 75
    maxshape = pars["shape"]  # "circle"
    levels = pars["levels"]
    nlevels = int( pars["nlevels"] )
    scale = pars["scale"]
    maxcntrs = int(pars["maxcontours"] ) # 1000


    # Load input image
    img2 = read_file(infile)
    vals = img2.get_image().values*1.0

    # Make list of contour levels
    grid = make_grid( vals, levels, nlevels, scale )

    # Create output file        
    out_mask = np.zeros_like( vals )
    img2.get_image().values = out_mask
    write_file( img2, outfile, clobber=True)

    t1,t2,t3,t4 = make_temp_files()
    use_img = infile
    maxval = max(grid)
    num_contours = 0
    stopit = (np.nan, np.nan)
    while num_contours < maxcntrs and maxval > min(grid):
        # Trying to do some of this with crates left too many files open

        dmstat( use_img, centroid=False, sigma=False, median=False)
        maxval = float( dmstat.out_max )
        if maxval <= min(grid):
            break
        
        # location of pixel == maxval
        xc = dmstat.out_max_loc.split(",")[0]
        yc = dmstat.out_max_loc.split(",")[1]
    
        # convert to logical coords
        dmcoords( infile, asol="", op="sky", x=xc, y=yc )
        lx = int(float(dmcoords.logicalx)+0.5)
        ly = int(float(dmcoords.logicaly)+0.5)

        if stopit == (lx,ly):
            # Hey,I've been here before.  Didn't work the last time
            # so just bail out.
            break
        else:
            stopit = (lx,ly)

        if "circle" == maxshape:
            ff = "[sky=circle({},{},{})][opt full]".format( xc,yc,maxrad)
        elif "box" == maxshape:
            ff = "[sky=box({0},{1},{2},{2})][opt full]".format( xc,yc,maxrad*2.0)
        else:
            raise ValueError("Unknown shape")

        low = grid[grid<maxval]
        if len(low) == 0:
            # pixels in image below lowest contour level, can't reach 'em
            break
        low = low[-1]

        num_contours += 1

        verb1("Contour lower limit: {} ({},{})".format(low,lx,ly))

        dmimglasso( use_img+ff, t1.name, xpos=lx, ypos=ly, coord="logical", 
            low_val=low, value="absolute", clobber=True)

        dmimgcalc( use_img+","+t1.name, "none", outfile=t3.name,
            op="imgout=(img1-(img1*img2))", clobber=True)

        shutil.move( t3.name, t2.name)
        use_img = t2.name
        
        dmimgcalc( outfile+","+t1.name, "none", outfile=t4.name+"[CONTMAP]", 
            op="imgout=(img1+(img2*{}))".format( num_contours), clobber=True)
        shutil.move( t4.name, outfile)

        
    # End while

    add_tool_history( outfile, toolname, pars, toolversion=__revision__)
    if len(pars["binimg"])>0 and "none" != pars["binimg"].lower():
        dmmaskbin( pars["infile"], outfile+"[opt type=i4]", pars["binimg"], clobber=True)
        add_tool_history( pars["binimg"], toolname, pars, toolversion=__revision__)


        

if __name__ == "__main__":
    try:
        main()
    except Exception as E:
        print("\n# "+toolname+" ("+__revision__+"): ERROR "+str(E)+"\n", file=sys.stderr)
        sys.exit(1)
    sys.exit(0)