#!/usr/bin/env python
#

"""
This file is part of ESA's XMM-Newton Scientific Analysis System (SAS).
SAS 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.

SAS 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 SAS.  If not, see .


Name    = 'slewFootCalc'
VERSION = '1.1';
author  = 'Ole Koenig, Jose Vicente Perea';
date    = '2018-07-30';


ChangeLog
=========

Version 1.1 - 2018-07-30 (JVP)
------------

+ Adaptation for Pipeline

Version 1.0 - 2018-06-26 (Ole Koenig)
------------

+ Initial version
  Task for Slew footprints
  This code aims on determining the footprints of the XMM exposure maps
  Author: Ole Koenig ole.koenig@fau.de 26/06/2018

"""

__version__   = "1.1"
__author__    = "Jose V. Perea"

__copyright__ = "ESA (C) 2000-2018"
__license__   = "GPL"
__progname__  = "slewFootCalc"

#------------------------------------------------------------------------------------------------

def usage():

#------------------------------------------------------------------------------------------------
   '''
      Help.
   '''
   print("""
usage:  slewFootCalc.py  --expMap=<...>  --ds9footprint=<...> 

This code aims on determining the footprints of the XMM exposure maps
Author: Ole Koenig ole.koenig@fau.de 26/06/2018

Basic idea of the code:
Axioms:
The XMM footprint consists of a minimum of 4 and a maximum of 6 points.
It has a clear rhombohedral shape with straight edges.
The four cornerpoints of the image frame are (clock-wise, in pixels):
      C1 = bottom left (1,1)
      C2 = upper left (1,max_y)
      C3 = upper right (max_x,max_y)
      C4 = bottom right (max_x,1)

1) Load the fits file
2) Detect NON-illuminated corner points (there are in total 4 corner points).
   Reason: By walking the edges from the NON-illuminated corners the code avoids
   landing in "gaps" or "stripes of 0 value" in the exposure map.
3) From the corner point walk in x and y direction through the arrays until 
   finding a value >0. This is the footprint point.
   Repeat this in a clock-wise rotation starting in the lower left corner.
4) In total there should be a maximum of 6 footprint points (and minimum of 4).
   If the code gives out less then 4 or more then 6 corner points it appends the
   filename to "list_of_weird_images.txt"
5) Tranform pixels to WCS world coordinates and write to file
   "xmm_expmaps_footprints.txt" in format filename | Polygon(ra1,dec1,ra2,dec2...)

                            ###### DS9 output ######
Additional one can set the option --ds9 which creates a .reg file for ds9
displaying the footprint points as polygon.
                            ########################

The code can be either used with an input file (list of filenames, delimiter=\n)
                python3 slewFootCalc.py --file filelist.txt [--ds9]

or with
                python3 slewFootCalc.py --file name.fits [--ds9] 

which promts the output to stdout.


What the code is NOT capable of:

- Find footprint points in the center of the image (it only walks on the edges)
- Calculate if the shape is a correct parallelogonal hexagon
- handle wavy lines within the image: it assumes a strict rhombohedral shape
- handle gaps (only walks on the edges of the image!)

Note that the x,y pixel definition in fits files runs from 1..n (actually, the
only really reasonable choice!) and not from 0..n-1 as it is in python.

Modified by:
- if either of the input image dimensions is 1 (or 0), exits without producing footprint.
- if the code cannot find the footprints by stepping just use the four corners of the 
input image

===================================================================================================
   """)

#------------------------------------------------------------------------------------------------
def error(string):
   '''
      prints the string in red.
   '''
   sys.stdout.write( string )					# No color string


def version(version_num):
   '''
      prints version number.
   '''
   versionString = '{0} ({0}-{1}) [-]'.format(sys.argv[0].split('/')[-1], version_num)
   print(versionString)



# ----------------------------------------------------------------------------------------
def getFootprint(filename, outFootprint):
# {{{
    ######## PARAMETERS ########
    # {{{
    # Pixel steps with which the algorithm walks. A larger threshold will
    # result in slightly faster computation, however more unprecise footprints
    thr = 1
    
    
    # If image_stopped = True the code will give error messages in case that one
    # edge of the image is illuminated and the opposite one NOT (in both x and y
    # direction. This was explicetly implemented to filter the xmm newton exposure
    # maps which ended on *EXPMAP8000.FTZ because they tend to stop in the middle of
    # the image. However, there are also many cases (see e.g. P9031400004PNS003EXPMAP8014.FTZ)
    # which fulfill this criterion (so I split the *8000 images from the rest and
    # calculated it separately.
    
    # RDS - not used in production code
    #image_stopped = False
    ############################
    # }}}

# main() returns two arguments:
# output  = fk5 ra,dec values of footprints
# damaged = if True the image is added to a "list_of_weird_images.txt"
    color_red      = '\x1b[6;15;41m'
    color_normal   = '\x1b[0m'

    from astropy.io import fits
    from astropy.wcs import WCS

    output, damaged = "",False
    
    # Load the exposure map (this is likely the most time consuming operation)
    with fits.open (filename) as hdulist:

      # Get x,y axis length (in pixel)
      NAXIS1 = hdulist[0].header['NAXIS1']
      NAXIS2 = hdulist[0].header['NAXIS2']

      # Load WCS information
      wcs =  WCS(hdulist[0].header)

      # Image[y][x] rows give x-coord, columns give y-coord
      image = hdulist[0].data
    #hdulist.close()

    # RDS - handle the case where one of the axes has only one pixel
    if NAXIS1<2 or NAXIS2<2:
       print("File " + filename + "has too few pixels in at least one dimension to "
                     + "calculate a footprint")
       damaged = True
       return " ",damaged

    # Define the four corner points
    C1, C2, C3, C4 = (0,0), (0,NAXIS2-1), (NAXIS1-1,NAXIS2-1), (NAXIS1-1,0)

    # At beginning all corner points are footprint points and then remove
    # illuminated pixels in get_start()
    startp = [C1,C2,C3,C4] # foot print points
    
    x_pixel_max = len (image[0]) -1
    y_pixel_max = len (image) -1

    # Get NON-illuminated pixels as starting points for walking section
    # ill are the remaining footprintpoints (illuminated corners) after deleting
    # the non-illuminated starting points (non_ill)
    non_ill, ill = get_start(image, x_pixel_max, y_pixel_max, C1, C2, C3, C4)

    # This code not run unless boolean changed at to pof code
    #if image_stopped == True:
        # First error handling: The *EXPMAP8000.FTZ tend to stop in the middle of
        # image. Since the code can't handle footprints within the image one has to
        # find a way to filter them and append to list_of_weird_images
        
        # If the the upper edge is illuminated somewhere, the bottom edge must be
        # illuminated, too! If this is not the case --> weird!
        #if max(image[0])==0.0 and max(image[-1])!=0.0 or max(image[0])!=0.0 and max(image[-1])==0.0:
         #   print("     " + color_red + "WARNING" + color_normal + ": One x-edge is illuminated, the other one not, appending to list of weird images.")
          #  damaged = True
    
        # Left edge illuminated and right edge not AND only one corner illuminated
        #elif image.max(axis=0)[0]==0.0 and image.max(axis=0)[-1]!=0.0 and len(ill)==1 or image.max(axis=0)[0]!=0.0 and image.max(axis=0)[-1]==0.0 and len(ill)==1:
         #   print("     " + color_red + "WARNING" + color_normal + ": One y-edge is illuminated, the other one not, appending to list of weird images.")
          #  damaged = True

    # Actual footprintfinding
    fpp=[]
    for ii, p in enumerate(non_ill,start=1):
        ####### WALKING FROM C1 #######
        if p=="C1":
        # Start at corner 1 (pixel 1/1) and walk first (y=0) row to the right
            for i in range(0, x_pixel_max, thr):
                if image[0][i] > 0.0:
                    fpp.append( (i,0) )
                    break

        if C1 in ill and C1 not in fpp and ii==1:
            fpp.append( C1 )
            
        if p=="C1":
        # Walk from C1 upwards (x=0)
            for i in range(0, y_pixel_max, thr):
                if image[i][0] > 0.0:
                    fpp.append( (0,i) )
                    break

                
        ####### WALKING FROM C2 #######
        if p=="C2":
        # Walk down on left edge (x=0)
            for i in reversed( range(0, y_pixel_max, thr)):
                if image[i][0] > 0.0:
                    fpp.append( (0,i) )
                    break
                
                
        if C2 in ill and C2 not in fpp and ii==2:
            fpp.append( C2 )
            
        if p=="C2":
        # Walk to the right on upper edge (y=1008)
            for i in range(0, x_pixel_max, thr):
                if image[-1][i] > 0.0:
                    fpp.append( (i,NAXIS2-1) )
                    break

                
        ####### WALKING FROM C3 #######
        if p=="C3":
        # Walk left on upper edge (y=1008)
            for i in reversed( range(0, x_pixel_max, thr)):
                if image[-1][i] > 0.0:
                    fpp.append( (i,NAXIS2-1) )
                    break

        if C3 in ill and C3 not in fpp and ii==3:
            fpp.append( C3 )

        if p=="C3":
        # Walk down on right edge (x=637)
            for i in reversed( range(0, y_pixel_max, thr)):
                if image[i][-1] > 0.0:
                    fpp.append( (NAXIS1-1,i) )
                    break


        ####### WALKING FROM C4 #######
        if p=="C4":
        # Walk up on right edge (x=637)
            for i in range(0, y_pixel_max, thr):
                if image[i][-1] > 0.0:
                    fpp.append( (NAXIS1-1,i) )
                    break

        if C4 in ill and C4 not in fpp and ii==4:
            fpp.append( C4 )
                
        if p=="C4":
        # Walk left on bottom edge (y=0)
            for i in reversed( range(0, x_pixel_max, thr)):
                if image[0][i] > 0.0:
                    fpp.append( (i,0) )
                    break

    #print ("Footprint points in pixels are: ", fpp)
                
    # Transform pixels to world coordinates if suitable number (at least 4) of footprints found
    if len(fpp) >= 4:# and len(fpp) <= 6:
        world = wcs.all_pix2world(fpp,0)

    else:
        # Error handling
        #damaged = True   - RDS can we do this ?
        if len(fpp)==0:
            # This means that no illuminated corners where found
            # --> The footprint is the whole image (all four corners)
            print("     " + color_red + "WARNING" + color_normal + ": No corner is illuminated! using image corners")
            world = wcs.all_pix2world(startp,0)
        elif 06
            # Less then 4 footprintpoints found. This is a weird shape!
            print("     " + color_red + "WARNING" + color_normal + ": Number of footprints smaller then 4, using image corners.")
            ## RDS - world = wcs.all_pix2world(fpp,0) # Default to four corners
            world = wcs.all_pix2world(startp,0)
            
#    # Print output in format filename | Polygon(ra1,dec1,ra2,dec2,...) |{(ra_rad1, dec_rad1),(ra_rad2, dec_rad2),...}
#    output += os.path.basename(filename) + " | Polygon("
#    for coord in world:
#        output += str(coord[0]) + "," + str(coord[1]) + ","
#    output = output[:len(output)-1] + ") |{" # remove the last ,
#    for coord in world:
#        output += "(" + str(coord[0]*pi/180) + "," + str(coord[1]*pi/180) + "),"
#    output=output[:len(output)-1] + "}" # Change the last , for a }
#
#    if make_ds9_file == True:
#        try:
#            os.mkdir ('ds9')
#        except OSError:
#            pass
#        #print ("     Creating region file.")

    writeds9output( outFootprint, world )

    
    return output, damaged
# }}}

# ----------------------------------------------------------------------------------------
# Function to determine the non-illuminated starting corners    
def get_start (image, x_pixel_max, y_pixel_max, C1, C2, C3, C4):
# {{{
    footprintpoints = [C1,C2,C3,C4]
    non_ill = ["","","",""]
    if image[0][0]==0.0:
        footprintpoints.remove( C1 )
        non_ill[0] = "C1"
    if image[y_pixel_max][0]==0.0:
        footprintpoints.remove( C2 )
        non_ill[1] = "C2"
    if image[y_pixel_max][x_pixel_max]==0.0:
        footprintpoints.remove( C3 )
        non_ill[2] = "C3"
    if image[0][x_pixel_max]==0.0:
        footprintpoints.remove( C4 )
        non_ill[3] = "C4"

    return non_ill, footprintpoints
# }}}

# ----------------------------------------------------------------------------------------
# Function to write a ds9 region file where footprint points are displayed
def writeds9output (outFootprint, fpp):
# {{{
    import os

    ### fp = open ('ds9/' + os.path.splitext (os.path.basename(filename))[0] + '_ds9.reg', 'w')
    print('{0}: Slew footprint ds9 file : {1}'.format(__progname__, outFootprint))
    fp = open (outFootprint, 'w')

    fp.write("# Region file format: DS9 version 4.0\n# Filename: %s\n" % outFootprint)
    fp.write("# global color=green font='helvetica 12 bold' select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\n")
    fp.write("linear\npolygon(")
    for i in range (len(fpp)):
        fp.write ("%f, %f," % (fpp[i][0], fpp[i][1])) # ra, dec
    #remove the comma at end of spoly
    cur_pos = fp.tell()
    fp.seek (cur_pos-1, os.SEEK_SET)
    fp.truncate()
    fp.write(") # color=green width=2")
    """
    # Circles
    fp.write("FK5\n")
    for i in range( len(fpp)):
        fp.write ("Circle(%f, %f, 0.005)\n" % (fpp[i][0], fpp[i][1])) # x,y (pixels)
    """
    fp.close()
# }}}


# ----------------------------------------------------------------------------------------
def main(arguments):
   '''
      This is the main function doing the full loop.

   '''
   import getopt, os

   import ModuleUtil
   
   version_num    = __version__		# Version number to be updated
   verbose        = False

   expMap         = ''
   outFootprint   = 'expMapfootprint.ds9'

   arguments = arguments.split()
   progName = arguments[0]
   argsList = arguments[1:]

   try:
      opts, args = getopt.getopt(argsList, "vo:Vh", ["expMap=", "outFootprint=", \
                                                   "version", "verbose=", "help"])
   except getopt.GetoptError as err:
      usage()
      error("ERROR: Could not understand the options!\n\n")
      print(err)  # will print something like "option -a not recognized"
      print("================================================================================================")
      return

   for o, a in opts:
      if o in ("-v", "--version"):
         version(version_num)
         return 0

      if o in ("-V", "--verbose"):
         verbose = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()

      if o in ("-h", "--help"):
         usage()
         return 0

      if o in ("--expMap"):
         expMap = a
      
      if o in ("-o", "--outFootprint"):
         outFootprint = a


   # Get the Slew footprint
   if os.path.isfile( expMap ) and os.access(expMap, os.R_OK):
      print('{0}: Creating footPrint for Slew exp map : {1}'.format(__progname__, expMap))
      output, damaged = getFootprint( expMap, outFootprint )
   else:
      usage()
      error("ERROR: Input exposure map does not exist\n\n")
      return

# ----------------------------------------------------------------------------------------




# ------------------------------------------------------------------------------------------------------------------------
if __name__=='__main__':
   
   import sys
   # start = time.time()
   
   main(" ".join(sys.argv))
   
   # print ("\nTotal runtime = ", time.time() - start, "s")
# ------------------------------------------------------------------------------------------------------------------------