/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA       *
 *                                                                              */

/*
 * $Author: amodigli $
 * $Date: 2010-09-24 09:32:03 $
 * $Revision: 1.10 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.8  2010/02/13 12:22:31  amodigli
 * removed inlines (let's do work to compiler)
 *
 * Revision 1.7  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.6  2007/05/02 13:17:23  jmlarsen
 * Allow specifying offset in optimal extraction
 *
 * Revision 1.5  2007/02/21 12:44:58  jmlarsen
 * uves_iterate_increment: Avoid recursion, fixed bug at beginning of new order (y limits were not recalculated)
 *
 * Revision 1.4  2006/11/20 08:01:46  jmlarsen
 * Changed format of pointer in uves_iterate_dump
 *
 * Revision 1.3  2006/11/16 09:48:30  jmlarsen
 * Renamed data type position -> uves_iterate_position, for namespace reasons
 *
 * Revision 1.2  2006/09/11 14:19:28  jmlarsen
 * Updated documentation
 *
 * Revision 1.1  2006/09/08 14:04:00  jmlarsen
 * Simplified code by using iterators, sky subtraction much optimized
 *
 *
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_extract_iterate   Iterate over an image
 *
 * This class exists in order to be able to write
 *
 *  @code
 *      for (uves_iterate_set_first(pos, ...);
 *           !uves_iterate_finished(pos);
 *           uves_iterate_increment(pos))
 *       {
 *             do something
 *       }
 *  @endcode
 *
 *  instead of repeating the same structure of nested loops + selections
 *  again and again in the extraction module. I.e. to avoid repetitions of
 *  (in pseudocode)
 *   
 *  @code
 *
 *  for each order ...
 *     for each bin ...
 *        calculate limits for this bin
 *        for each position on the slit ...
 *           if pixel is good ...
 *
 *              do something
 *
 *  @endcode
 *
 *  where only the "do something" part is different
 *
 */
/*----------------------------------------------------------------------------*/
/**@{*/


/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/

#include <uves_extract_iterate.h>

#include <uves_utils.h>

#include <cpl.h>

/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/
static
bool illegal_position(const uves_iterate_position *p);

/*-----------------------------------------------------------------------------
                            Implementation
 -----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------*/
/**
   @brief    Allocate iterator
   @param    image           the image to iterate over
   @param    order_locations order polynomial
   @param    minorder        first order (inclusive) in image
   @param    maxorder        last order (inclusive) in image
   @param    slit_geometry   extraction slit length and offset
   @return   newly allocated iterator which must be deallocated using
             uves_iterate_delete

   The iterator, pos, should be used in constructs like
   @code
    for (uves_iterate_set_first(pos,
                    xmin,
                    xmax,
                    ordermin, ordermax,
                    current_image_bpm, true);
         !uves_iterate_finished(pos);
         uves_iterate_increment(pos))
         {
                flux += image[p->x, p->y];
         }
   @endcode

   The struct members

    int order;
    int x, y;
    double ycenter;
    int ylow;
    int yhigh;

   can be accessed inside the loop.

*/
/*----------------------------------------------------------------------------*/
uves_iterate_position *
uves_iterate_new(int nx, int ny,
         const polynomial *order_locations,
         int minorder, int maxorder,
         slit_geometry sg)
{
    uves_iterate_position *p = cpl_calloc(1, sizeof(uves_iterate_position));

    p->nx = nx;
    p->ny = ny;
    p->order_locations = order_locations;
    p->minorder = minorder;
    p->maxorder = maxorder;
    p->sg = sg;

    return p;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Deallocate iterator and set pointer to NULL
   @param    p             iterator to delete
*/
/*----------------------------------------------------------------------------*/
void
uves_iterate_delete(uves_iterate_position **p)
{
    if (p != NULL)
    {
        cpl_free(*p);
        *p = NULL;
    }
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Initialize iteration
   @param    p             iterator to initialize
   @param    xmin          minimum x bin
   @param    xmax          maximum x bin
   @param    ordermin      first order to iterate over
   @param    ordermax      last order to iterate over
   @param    bpm           image bad pixel map to use for this iteration. If
                           NULL all pixels are considered good
   @param    loop_y        if true, the iteration is over (order, x, y).
                           If false, the iteration is over (order, x) only,
               and the bpm is not used
*/
/*----------------------------------------------------------------------------*/
void
uves_iterate_set_first(uves_iterate_position *p,
               int xmin, int xmax, 
               int ordermin, int ordermax,
               const cpl_binary *bpm,
               bool loop_y)
{
    /* Limits for this iteration */
    p->xmin = xmin;
    p->xmax = xmax;
    p->ordermax = ordermax;
    p->bpm = bpm;
    p->loop_y = loop_y;
    p->end = false;

    /* Set first postion */
    p->x = xmin;
    p->order = ordermin;

    p->ycenter = uves_polynomial_evaluate_2d(p->order_locations, p->x, p->order)
        + p->sg.offset;
    p->yhigh = uves_round_double(p->ycenter + p->sg.length/2);
    p->ylow  = uves_round_double(p->ycenter - p->sg.length/2);
    if (loop_y) 
    {
        p->y = p->ylow;
    }

    /* Go to first good pixel */
    while (illegal_position(p) && !uves_iterate_finished(p))
    {
        uves_iterate_increment(p);
    }
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Get next position
   @param    p             iterator to increment

   p is changed to contain the next position
*/
/*----------------------------------------------------------------------------*/
void
uves_iterate_increment(uves_iterate_position *p)
{
    do {
    if (p->loop_y && p->y < p->yhigh)
        {
        (p->y)++;
        }
    else if (p->x < p->xmax)
        {
        (p->x)++;
        
        p->ycenter = 
                    uves_polynomial_evaluate_2d(p->order_locations, 
                                                p->x, p->order)
                    + p->sg.offset;

        p->yhigh = uves_round_double(p->ycenter + p->sg.length/2);
        p->ylow  = uves_round_double(p->ycenter - p->sg.length/2);
        if (p->loop_y) p->y = p->ylow;
        }
    else if (p->order < p->ordermax)
        {
        (p->order)++;
        p->x = p->xmin;
        
        p->ycenter = 
                    uves_polynomial_evaluate_2d(p->order_locations,
                                                p->x, p->order)
                    + p->sg.offset;

        p->yhigh = uves_round_double(p->ycenter + p->sg.length/2);
        p->ylow  = uves_round_double(p->ycenter - p->sg.length/2);
        if (p->loop_y) p->y = p->ylow;
        }
    else
        {
        p->end = true;
        }
    } while (illegal_position(p) && !uves_iterate_finished(p));

    return;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Finished iterating?
   @param    p             iterator to check
   @return false if p contains a valid position, false if the iteration
   is over

   This function should always be checked after calling uves_iterate_increment()
   and uves_iterate_set_first()
*/
/*----------------------------------------------------------------------------*/
bool
uves_iterate_finished(const uves_iterate_position *p)
{
    return p->end;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    dump iterator (for debugging)
   @param    p             iterator to check
   @param    stream        stream to write to
*/
/*----------------------------------------------------------------------------*/
void
uves_iterate_dump(const uves_iterate_position *p, FILE *stream)
{
    fprintf(stream, "Position:\n");
    fprintf(stream, "order       = %d\n", p->order);
    fprintf(stream, "x           = %d\n", p->x);
    fprintf(stream, "y           = %d\n", p->y);
    fprintf(stream, "ycenter     = %f\n", p->ycenter);
    fprintf(stream, "ylow, yhigh = %d, %d\n", p->ylow, p->yhigh);
    fprintf(stream, "Limits:\n");
    fprintf(stream, "xmin, xmax = %d, %d\n", p->xmin, p->xmax);
    fprintf(stream, "ordermax   = %d\n", p->ordermax);
    fprintf(stream, "bpm        = %d\n", p->bpm != NULL ? 1 : 0);
    fprintf(stream, "loop_y     = %s\n", p->loop_y ? "true" : "false");
    fprintf(stream, "end        = %s\n", p->end    ? "true" : "false");
    fprintf(stream, "Geometry:\n");
    fprintf(stream, "nx, ny             = %d, %d\n", p->nx, p->ny);
    fprintf(stream, "minorder, maxorder = %d, %d\n", p->minorder, p->maxorder);
    fprintf(stream, "order_locations    = %d\n", p->order_locations != NULL ? 1 : 0);
    fprintf(stream, "slit length        = %f\n", p->sg.length);
    fprintf(stream, "slit offset        = %f\n", p->sg.offset);

    return;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Find out if the current position is valid
   @param    p             iterator to check
   @return true iff the current position is a good pixel and the entire
   slit is inside the image
*/
/*----------------------------------------------------------------------------*/
static
bool illegal_position(const uves_iterate_position *p)
{
    return p->ylow < 1 || p->yhigh > p->ny ||
    (p->loop_y && p->bpm != NULL &&
     p->bpm[(p->x-1) + (p->y-1)*p->nx] != CPL_BINARY_0);
}

/**@}*/