/* $Id: cpl_detector.c,v 1.10 2008/02/26 09:53:27 yjung Exp $
 *
 * This file is part of the ESO Common Pipeline Library
 * Copyright (C) 2001-2004 European Southern 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: yjung $
 * $Date: 2008/02/26 09:53:27 $
 * $Revision: 1.10 $
 * $Name:  $
 */

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

/*-----------------------------------------------------------------------------
                                   Define
 -----------------------------------------------------------------------------*/

#define CPL_CLASS_NONE      0
#define CPL_CLASS_DOUBLE    1
#define CPL_CLASS_FLOAT     2
#define CPL_CLASS_INT       3

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <math.h>

#include <cpl_msg.h>
#include <cpl_stats.h>
#include <cpl_mask.h>
#include <cpl_matrix.h>
#include <cpl_tools.h>
#include <cpl_math_const.h>

#include "cpl_detector.h"

/*-----------------------------------------------------------------------------
                                   Define
 -----------------------------------------------------------------------------*/

#define CPL_DET_CLEAN_BAD_PIX   0

/* Computes the square of an euclidean distance bet. 2 points */
#define pdist(x1,y1,x2,y2) (((x1-x2)*(x1-x2))+((y1-y2)*(y1-y2)))
/** Computes the square of an euclidean distance bet. 2 points (polar coord) */
#define qdist(r1,t1,r2,t2) \
  ((r1*r1)+(r2*r2)-2*r1*r2*cos((t1-t2)*CPL_MATH_PI_2/90))

/*-----------------------------------------------------------------------------
                            Function codes
 -----------------------------------------------------------------------------*/

static cpl_bivector * cpl_bivector_gen_rect_poisson(const int *, const int,
        const int) ;
static cpl_bivector * cpl_bivector_gen_ring_poisson(const double *, const int,
        const int) ;

/*----------------------------------------------------------------------------*/
/**
 * @defgroup cpl_detector   High-level functions to compute detector features
 *
 * @par Synopsis:
 * @code
 *   #include "cpl_detector.h"
 * @endcode
 */
/*----------------------------------------------------------------------------*/
/**@{*/

/*-----------------------------------------------------------------------------
                            Function codes
 -----------------------------------------------------------------------------*/

#define RECT_RON_HS         4
#define RECT_RON_SAMPLES    100
/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the readout noise in a rectangle.
  @param    diff        Input image, usually a difference frame.
  @param    zone_def    Zone where the readout noise is to be computed.
  @param    ron_hsize   to specify half size of squares (<0 to use default)
  @param    ron_nsamp   to specify the nb of samples (<0 to use default)
  @param    noise       Output parameter: noise in the frame.
  @param    error       Output parameter: error on the noise.
  @return   the #_cpl_error_code_ or CPL_ERROR_NONE
  @note     This function uses the pseudo-random number generator.

  This function is meant to compute the readout noise in a frame by means of a
  MonteCarlo approach. The input is a frame, usually a difference between two
  frames taken with the same settings for the acquisition system, although no
  check is done on that, it is up to the caller to feed in the right kind of
  frame.

  The provided zone is an array of four integers specifying the zone to take
  into account for the computation. The integers specify ranges as xmin, xmax,
  ymin, ymax, where these coordinates are given in the FITS notation (x from 1
  to lx, y from 1 to ly and bottom to top). Specify NULL instead of an array of
  four values to use the whole frame in the computation.

  The algorithm will create typically 100 9x9 windows on the frame, scattered
  optimally using a Poisson law. In each window, the standard deviation of all
  pixels in the window is computed and this value is stored. The readout noise
  is the median of all computed standard deviations, and the error is the
  standard deviation of the standard deviations.

  Both noise and error are returned by modifying a passed double. If you do
  not care about the error, pass NULL.

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT if diff or noise is NULL
  - CPL_ERROR_ILLEGAL_INPUT if the specified window (zone_def) is invalid
 */
/*----------------------------------------------------------------------------*/
cpl_error_code cpl_flux_get_noise_window(
        const cpl_image *   diff,
        const int       *   zone_def,
        int                 ron_hsize,
        int                 ron_nsamp,
        double          *   noise,
        double          *   error)
{
    const int hsize    = ron_hsize < 0 ? RECT_RON_HS      : ron_hsize;
    const int nsamples = ron_nsamp < 0 ? RECT_RON_SAMPLES : ron_nsamp;
    cpl_bivector    *   sample_reg ;
    cpl_vector      *   rms_list ;
    int                 rect[4] ;
    int                 zone[4];
    double          *   px;
    double          *   py;
    double          *   pr;
    int                 i ;

    /* Test entries */
    cpl_ensure_code(diff && noise, CPL_ERROR_NULL_INPUT) ;

    /* Generate nsamples window centers in the image */
    if (zone_def != NULL) {
        rect[0] = zone_def[0] + hsize + 1 ; /* xmin */
        rect[1] = zone_def[1] - hsize - 1 ; /* xmax */
        rect[2] = zone_def[2] + hsize + 1 ; /* ymin */
        rect[3] = zone_def[3] - hsize - 1 ; /* ymax */
    } else {
        rect[0] = hsize + 1 ;                /* xmin */
        rect[1] = cpl_image_get_size_x(diff) - hsize - 1 ;    /* xmax */
        rect[2] = hsize + 1 ;                /* ymin */
        rect[3] = cpl_image_get_size_y(diff) - hsize - 1 ;    /* ymax */
    }

    cpl_ensure_code( rect[0] < rect[1] && rect[2] < rect[3],
                CPL_ERROR_ILLEGAL_INPUT) ;

    /* Generate n+1 regions, because the first region is always at (0,0) */
    /* and it would bias the measurement. */
    sample_reg = cpl_bivector_gen_rect_poisson(rect, nsamples+1, nsamples+1);
    cpl_ensure( sample_reg != NULL, CPL_ERROR_ILLEGAL_INPUT,
                CPL_ERROR_ILLEGAL_INPUT);

    px = cpl_bivector_get_x_data(sample_reg);
    py = cpl_bivector_get_y_data(sample_reg);

    /* Now, for each window center, extract a vignette and compute the */
    /* signal RMS in it. Store this rms into a table. */
    rms_list = cpl_vector_new(nsamples);
    cpl_ensure( rms_list != NULL, CPL_ERROR_NULL_INPUT,
                CPL_ERROR_NULL_INPUT) ;
    pr = cpl_vector_get_data(rms_list);

    for (i=0 ; i<nsamples ; i++) {
        zone[0] = (int)px[i+1] - hsize ;
        zone[1] = (int)px[i+1] + hsize ;
        zone[2] = (int)py[i+1] - hsize ;
        zone[3] = (int)py[i+1] + hsize ;
        pr[i] = cpl_image_get_stdev_window(diff, 
                zone[0], zone[2], zone[1], zone[3]) ;
    }
    cpl_bivector_delete(sample_reg);

    /* The error is the rms of the rms */
    if (error != NULL) *error = cpl_vector_get_stdev(rms_list);

    /* The final computed RMS is the median of all values.  */
    /* This call will modify the rms_list */
    *noise = cpl_vector_get_median(rms_list) ;

    cpl_vector_delete(rms_list);

    return CPL_ERROR_NONE;
}
#undef RECT_RON_HS
#undef RECT_RON_SAMPLES

#define RING_RON_HS         4    
#define RING_RON_SAMPLES    50
/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the readout noise in a ring.
  @param    diff    Input image, usually a difference frame.
  @param    zone_def    Zone where the readout noise is to be computed.
  @param    ron_hsize   to specify half size of squares (<0 to use default)
  @param    ron_nsamp   to specify the nb of samples (<0 to use default)
  @param    noise       Output parameter: noise in the frame.
  @param    error       Output parameter: error on the noise.
  @return   the #_cpl_error_code_ or CPL_ERROR_NONE
  @see      cpl_flux_get_noise_window 
  @note     This function uses the pseudo-random number generator.

  The provided zone is an array of four integers specifying the zone to take
  into account for the computation. The integers specify a ring as x, y,
  r1, r2 where these coordinates are given in the FITS notation (x from 1
  to nx, y from 1 to ny and bottom to top). 

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT if diff noise or zone_def is NULL
  - CPL_ERROR_ILLEGAL_INPUT if the internal radius is bigger than the
    external one in zone_def
 */
/*----------------------------------------------------------------------------*/
cpl_error_code cpl_flux_get_noise_ring(
        const cpl_image *   diff,
        const double    *   zone_def,
        int                 ron_hsize,
        int                 ron_nsamp,
        double          *   noise,
        double          *   error)
{
    cpl_bivector    *   sample_reg ;
    double          *   px;
    double          *   py;
    int                 zone[4];
    const int           hsize = ron_hsize < 0 ? RING_RON_HS      : ron_hsize;
    const int           nsamp = ron_nsamp < 0 ? RING_RON_SAMPLES : ron_nsamp;
    cpl_vector      *   rms_list;
    int                 nx, ny;
    int                 i, rsize;

    cpl_ensure_code(diff,     CPL_ERROR_NULL_INPUT) ;
    cpl_ensure_code(noise,    CPL_ERROR_NULL_INPUT) ;
    cpl_ensure_code(zone_def, CPL_ERROR_NULL_INPUT) ;
    cpl_ensure_code(zone_def[2] < zone_def[3], CPL_ERROR_ILLEGAL_INPUT);

    nx = cpl_image_get_size_x(diff);
    ny = cpl_image_get_size_y(diff);

    cpl_ensure_code(2 * ron_hsize < nx, CPL_ERROR_ILLEGAL_INPUT);
    cpl_ensure_code(2 * ron_hsize < ny, CPL_ERROR_ILLEGAL_INPUT);

    /* Generate n+1 regions, because the first region is always at a given */
    /* position and it would bias the measurement. */
    sample_reg = cpl_bivector_gen_ring_poisson(zone_def, nsamp+1, nsamp+1);
    cpl_ensure_code(sample_reg != NULL, CPL_ERROR_NULL_INPUT);
    px = cpl_bivector_get_x_data(sample_reg);
    py = cpl_bivector_get_y_data(sample_reg);

    /* Now, for each valid sample, extract a vignette and compute the RMS */
    rms_list = cpl_vector_new(nsamp + 1);
    cpl_ensure_code(rms_list, CPL_ERROR_NULL_INPUT);
    rsize = 0 ;
    for (i=0 ; i < nsamp + 1 ; i++) {
        /* Add 0.5 to round off to nearest (positive) integer */
        const double dx = 0.5 + zone_def[0] + px[i] * cos(py[i]);
        const double dy = 0.5 + zone_def[1] + px[i] * sin(py[i]);

        zone[0] = (int)(dx - hsize);
        if (zone[0] < 1) continue;

        zone[1] = (int)(dx + hsize); 
        if (zone[1] > nx) continue;

        zone[2] = (int)(dy - hsize);
        if (zone[2] < 1) continue;

        zone[3] = (int)(dy + hsize);
        if (zone[3] > ny) continue;

        cpl_vector_set(rms_list, rsize, cpl_image_get_stdev_window(diff,
                                                                   zone[0],
                                                                   zone[2],
                                                                   zone[1],
                                                                   zone[3]));
        rsize++;
        /* At this point rms_list has rsize elements */
    }

    cpl_bivector_delete(sample_reg);

    if (4 * rsize < nsamp) {
        cpl_vector_delete(rms_list);
        return cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
                                           __FILE__, __LINE__,
                                           "Need at least %d/4 (not %d) samples"
                                           " to compute noise", nsamp, rsize);
    }

    /* Ignore all but the first rsize elements */
    rms_list = cpl_vector_wrap(rsize, (double*)cpl_vector_unwrap(rms_list));

    /* The error is the rms of the rms */
    if (error != NULL) *error = cpl_vector_get_stdev(rms_list);
    
    /* The final computed RMS is the median of all values.  */
    /* This call will modify the rms_list */
    *noise = cpl_vector_get_median(rms_list) ;
    
    cpl_vector_delete(rms_list);

    return CPL_ERROR_NONE;
}
#undef RING_RON_HLX
#undef RING_RON_HLY
#undef RING_RON_SAMPLES

#define CPL_OPERATION CPL_DET_CLEAN_BAD_PIX
/*----------------------------------------------------------------------------*/
/**
  @brief    Clean the bad pixels in an image
  @param    in      The image to clean
  @return   the #_cpl_error_code_ or CPL_ERROR_NONE

  Images can be CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE.
  
  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL
  - CPL_ERROR_TYPE_MISMATCH if the passed image list type is not supported
 */
/*----------------------------------------------------------------------------*/
cpl_error_code cpl_detector_interpolate_rejected(cpl_image * in) 
{
    int                 nx, ny ;
    int                 nbad ;
    cpl_mask        *   bpm ;
    cpl_binary      *   pbpm ;
    const cpl_mask  *   bpm_cur ;
    const cpl_binary*   pbpm_cur ;
    int                 pos ;
    double              sum_pix ;
    int                 nb_pix ;
    int                 i, j ;
        
    /* Initialise */
    nx = cpl_image_get_size_x(in) ;
    ny = cpl_image_get_size_y(in) ;
    nbad = cpl_image_count_rejected(in) ;

    /* Check entries */
    cpl_ensure_code(in, CPL_ERROR_NULL_INPUT) ;
    if (nbad==0) return CPL_ERROR_NONE ;

    /* Switch on image type */
    switch (cpl_image_get_type(in)) {
#define CPL_CLASS CPL_CLASS_DOUBLE
#include "cpl_detector_body.h"
#undef CPL_CLASS
#define CPL_CLASS CPL_CLASS_FLOAT
#include "cpl_detector_body.h"
#undef CPL_CLASS
        default:
          cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH) ;
    }
    return CPL_ERROR_NONE ;
}
#undef CPL_OPERATION

/**@}*/

/*----------------------------------------------------------------------------*/
/*
  @brief    Generate points with a Poisson scattering property in a rectangle.
  @param    r       Array of 4 integers as [xmin,xmax,ymin,ymax]
  @param    np      Number of points to generate.
  @param    homog   Homogeneity factor.
  @return   Newly allocated cpl_bivector object or NULL in error case

  The returned object must be deallocated using cpl_bivector_delete().

  <b>POISSON POINT GENERATION</b>

  Without homogeneity factor, the idea is to generate a set of np points within
  a given rectangle defined by (xmin xmax ymin ymax).
  All these points obey a Poisson law, i.e. no couple of points is closer to
  each other than a minimal distance. This minimal distance is defined as a
  function of the input requested rectangle and the requested number of points
  to generate. We apply the following formula:

  @f$ d_{min} = \sqrt{\frac{W \times H}{\sqrt{2}}} @f$

  Where W and H stand for the rectangle width and height. Notice that the
  system in which the rectangle vertices are given is completely left
  unspecified, generated points will have coordinates in the specified x and y
  ranges.

  With a specified homogeneity factor h (0 < h <= np), the generation algorithm
  is different.

  the Poisson law applies for any h consecutive points in the final output, but
  not for the whole point set. This enables us to generate groups of points
  which statisfy the Poisson law, without constraining the whole set. This
  actually is equivalent to dividing the rectangle in h regions of equal
  surface, and generate points randomly in each of these regions, changing
  region at each point.

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL
  - CPL_ERROR_ILLEGAL_INPUT if np is not strictly positive
 */
/*----------------------------------------------------------------------------*/
static cpl_bivector * cpl_bivector_gen_rect_poisson(
        const int   *   r,
        const int       np,
        const int       homog)
{
    double              min_dist ;
    int                 i ;
    int                 gnp ;
    cpl_bivector    *   list ;
    double              cand_x, cand_y ;
    int                 ok ;
    int                 start_ndx ;
    int                 xmin, xmax, ymin, ymax ;
    /* Corrected Homogeneity factor */
    const int           homogc = 0 < homog && homog < np ? homog : np;
    double          *   px ;
    double          *   py ;

    /* error handling: test arguments are correct */
    cpl_ensure(r, CPL_ERROR_NULL_INPUT, NULL) ;
    cpl_ensure( np > 0, CPL_ERROR_ILLEGAL_INPUT, NULL) ;

    list = cpl_bivector_new(np);
    cpl_ensure(list, CPL_ERROR_NULL_INPUT,NULL) ;
    px = cpl_bivector_get_x_data(list) ;
    py = cpl_bivector_get_y_data(list) ;

    xmin = r[0] ;
    xmax = r[1] ;
    ymin = r[2] ;
    ymax = r[3] ;

    min_dist   = CPL_MATH_SQRT1_2*((xmax-xmin)*(ymax-ymin) /(double)(homogc+1));
    gnp        = 1 ;
    px[0]      = 0 ;
    py[0]      = 0 ;

    /* First: generate <homog> points */
    while (gnp < homogc) {
        /* Pick a random point within requested range */
        cand_x = cpl_drand() * (xmax - xmin) + xmin ;
        cand_y = cpl_drand() * (ymax - ymin) + ymin ;

        /* Check the candidate obeys the minimal Poisson distance */
        ok = 1 ;
        for (i=0 ; i<gnp ; i++) {
            if (pdist(cand_x, cand_y, px[i], py[i])<min_dist) {
                /* does not check Poisson law: reject point */
                ok = 0 ;
                break ;
            }
        }
        if (ok) {
            /* obeys Poisson law: register the point as valid */
            px[gnp] = cand_x ;
            py[gnp] = cand_y ;
            gnp ++ ;
        }
    }

    /* Iterative process: */
    /* Pick points out of Poisson distance of the last <homogc-1> points. */
    start_ndx=0 ;
    while (gnp < np) {
        /* Pick a random point within requested range */
        cand_x = cpl_drand() * (xmax - xmin) + xmin ;
        cand_y = cpl_drand() * (ymax - ymin) + ymin ;

        /* Check the candidate obeys the minimal Poisson distance */
        ok = 1 ;
        for (i=0 ; i<homogc ; i++) {
            if (pdist(cand_x,
                      cand_y,
                      px[start_ndx+i],
                      py[start_ndx+i]) < min_dist) {
                /* does not check Poisson law: reject point */
                ok = 0 ;
                break ;
            }
        }
        if (ok) {
            /* obeys Poisson law: register the point as valid */
            px[gnp] = cand_x ;
            py[gnp] = cand_y ;
            gnp ++ ;
        }
    }

    /* Iterative process: */
    /* Pick points out of Poisson distance of the last <homogc-1> points. */
    start_ndx=0 ;
    while (gnp < np) {
        /* Pick a random point within requested range */
        cand_x = cpl_drand() * (xmax - xmin) + xmin ;
        cand_y = cpl_drand() * (ymax - ymin) + ymin ;

        /* Check the candidate obeys the minimal Poisson distance */
        ok = 1 ;
        for (i=0 ; i<homogc ; i++) {
            if (pdist(cand_x,
                      cand_y,
                      px[start_ndx+i],
                      py[start_ndx+i]) < min_dist) {
                /* does not check Poisson law: reject point */
                ok = 0 ;
                break ;
            }
        }
        if (ok) {
            /* obeys Poisson law: register the point as valid */
            px[gnp] = cand_x ;
            py[gnp] = cand_y ;
            gnp ++ ;
            start_ndx ++ ;
        }
    }
    return list ;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Generate points with a Poisson scattering property in a ring.
  @param    r       Array of 4 doubles as [x,y,r1,r2]
  @param    np      Number of points to generate.
  @param    homog   Homogeneity factor.
  @return   Newly allocated cpl_bivector object or NULL in error case
  @see      cpl_bivector_gen_rect_poisson

  The returned object must be deallocated using cpl_bivector_delete().

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL
  - CPL_ERROR_ILLEGAL_INPUT if np is not strictly positive
 */
/*----------------------------------------------------------------------------*/
static cpl_bivector * cpl_bivector_gen_ring_poisson(
        const double    *   r,
        const int           np,
        const int           homog)
{
    double              min_dist ;
    int                 gnp ;
    cpl_bivector    *   list ;
    double              cand_r, cand_t ;
    int                 ok ;
    int                 start_ndx ;
    int                 r1, r2 ;
    int                 i ;
    /* Corrected Homogeneity factor */
    const int           homogc = 0 < homog && homog < np ? homog : np;
    double          *   px ;
    double          *   py ;

    /* Check entries */
    cpl_ensure(r, CPL_ERROR_NULL_INPUT, NULL) ;
    cpl_ensure(np > 0, CPL_ERROR_ILLEGAL_INPUT, NULL) ;

    list = cpl_bivector_new(np);
    cpl_ensure(list, CPL_ERROR_NULL_INPUT,NULL) ;
    px = cpl_bivector_get_x_data(list) ;
    py = cpl_bivector_get_y_data(list) ;

    r1 = (int)(r[2]) ;
    r2 = (int)(r[3]) ;

    min_dist    = (CPL_MATH_PI_2/CPL_MATH_SQRT1_2)*((r2*r2)-(r1*r1)) /
        (double)(homogc+1);
    gnp         = 1 ;
    px[0]       = r1 ;
    py[0]       = 0 ;

    /* First: generate <homog> points */
    while (gnp < homogc) {
        /* Pick a random point within requested range */
        cand_r = cpl_drand() * (r2 - r1) + r1 ;
        cand_t = cpl_drand() * 360 ;

        /* Check the candidate obeys the minimal Poisson distance */
        ok = 1 ;
        for (i=0 ; i<gnp ; i++) {
            if (qdist(cand_r, cand_t, px[i], py[i])<min_dist) {
                /* does not check Poisson law: reject point */
                ok = 0 ;
                break ;
            }
        }
        if (ok) {
            /* obeys Poisson law: register the point as valid */
            px[gnp] = cand_r ;
            py[gnp] = cand_t ;
            gnp ++ ;
        }
    }

    /* Iterative process: */
    /* Pick points out of Poisson distance of the last <homogc-1> points. */
    start_ndx=0 ;
    while (gnp < np) {
        /* Pick a random point within requested range */
        cand_r = cpl_drand() * (r2 - r1) + r1 ;
        cand_t = cpl_drand() * 360 ;

        /* Check the candidate obeys the minimal Poisson distance */
        ok = 1 ;
        for (i=0 ; i<homogc ; i++) {
            if (qdist(cand_r,
                      cand_t,
                      px[start_ndx+i],
                      py[start_ndx+i]) < min_dist) {
                /* does not check Poisson law: reject point */
                ok = 0 ;
                break ;
            }
        }
        if (ok) {
            /* obeys Poisson law: register the point as valid */
            px[gnp] = cand_r ;
            py[gnp] = cand_t ;
            gnp ++ ;
            start_ndx ++ ;
        }
    }
    return list ;
}