/* $Id: fors_extract.c,v 1.45 2011-10-13 14:29:24 cgarcia Exp $
 *
 * This file is part of the FORS Library
 * Copyright (C) 2002-2010 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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
 */

/*
 * $Author: cgarcia $
 * $Date: 2011-10-13 14:29:24 $
 * $Revision: 1.45 $
 * $Name: not supported by cvs2svn $
 */

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

#include <fors_extract.h>
#include <fors_star.h>

#include <fors_tools.h>
#include <fors_dfs.h>
#include <fors_pfits.h>
#include <fors_utils.h>

#include <cpl.h>

#include <string.h>
#include <stdbool.h>
#include <math.h>

/**
 * @defgroup fors_extract  Image source extraction
 */

/**@{*/

/*-----------------------------------------------------------------------------
    (Proto)types
 -----------------------------------------------------------------------------*/

struct _extract_method
{
    enum {SEX, TEST} method;
    const char *sex_exe;
    const char *sex_config;
    const char *sex_mag;
    const char *sex_magerr;
    int sex_radius;
};

static fors_star_list *
extract_sex(                                const fors_image *image,
                                            const fors_setting *setting,
                                            const char *sex_exe,
                                            const char *sex_config,
                                            const char *sex_mag,
                                            const char *sex_magerr,
                                            int radius,
                                            double magsyserr,
                                            fors_extract_sky_stats *sky_stats,
                                            cpl_image **background,
                                            cpl_table **extracted_sources);

static fors_star_list *
extract_test(                               fors_extract_sky_stats *sky_stats,
                                            cpl_image **background,
                                            cpl_table **extracted_sources);

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

/*----------------------------------------------------------------------------*/
/**
 * @brief   Check SExtractor object feature: FLAG
 * @param   SExtractor object's binary contamination flags
 * @return  1 if successful, 0 on failure
 * 
 * All FLAGS 1, 2, 4, ..., 128 are severe enough that
 * we do not want to use the source for photometry
 * see SExtractor doc. for the meaning of each flag.
 */
/*----------------------------------------------------------------------------*/
bool
fors_extract_check_sex_flag(                unsigned int    sex_flag)
{
    return (sex_flag == 0x0);
}

/*----------------------------------------------------------------------------*/
/**
 * @brief   Check a sextracted star for validity
 * @param   star    (S)Extracted star
 * @param   ref_img (Optional) reference image (to check for image range)
 * @return  1 if successful, 0 on failure
 * 
 * The following SExtractor failures are caught:
 * - Check for allowed general parameters using fors_star_check_values()
 * - Invalid SExtractor magnitude (can be 99.0)
 * - (x, y) out of range (if @a ref_img != NULL)
 * 
 * @todo
 * - FIXME: FAP: verify the criterion for magnitude rejection
 */
/*----------------------------------------------------------------------------*/
bool
fors_extract_check_sex_star(                const fors_star *star,
                                            const cpl_image *ref_img)
{
    bool    success = 1;
    
    if (star == NULL)
        return 0;
    
    success &= fors_star_check_values(star);
    
    success &= (star->magnitude < 98);
    
    if (ref_img != NULL)
    {
        success &= star->pixel->x >= 1;
        success &= star->pixel->x <= cpl_image_get_size_x(ref_img);
        success &= star->pixel->y >= 1;
        success &= star->pixel->y <= cpl_image_get_size_y(ref_img);
    }
    
    return success;
}

/*----------------------------------------------------------------------------*/
/**
 * @brief    Define recipe parameters
 * @param    parameters     parameter list to fill
 * @param    context        parameters context
 */
/*----------------------------------------------------------------------------*/
void 
fors_extract_define_parameters(             cpl_parameterlist *parameters, 
                                            const char *context)
{
    cpl_parameter *p;
    const char *full_name = NULL;
    const char *name;
    
    /*
    name = "extract_method";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_enum(full_name,
                               CPL_TYPE_STRING,
                               "Source extraction method",
                               context,
                               "sex", 2,
                               "sex", "test");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name);
    */

    name = "sex_exe";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_value(full_name,
                                CPL_TYPE_STRING,
                                "SExtractor executable",
                                context,
                                FORS_SEXTRACTOR_PATH "/sex");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name);
    
    name = "sex_config";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_value(full_name,
                                CPL_TYPE_STRING,
                                "SExtractor configuration file",
                                context,
                                FORS_SEXTRACTOR_CONFIG "/fors.sex");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name);
    
    
    name = "sex_mag";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_value(full_name,
                                CPL_TYPE_STRING,
                                "SExtractor magnitude",
                                context,
                                "MAG_APER");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name);
    
    name = "sex_magerr";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_value(full_name,
                                CPL_TYPE_STRING,
                                "SExtractor magnitude error",
                                context,
                                "MAGERR_APER");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name);

    name = "sex_radius";
    full_name = cpl_sprintf("%s.%s", context, name);
    p = cpl_parameter_new_value(full_name,
                                CPL_TYPE_INT,
                                "Background error map median filter "
                                "radius (unbinned pixels)",
                                context,
                                64);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(parameters, p);
    cpl_free((void *)full_name); full_name = NULL;
    
    return;
}

/*----------------------------------------------------------------------------*/
#undef cleanup
#define cleanup \
do { \
    cpl_free((void *)name); \
    cpl_free((void *)method); \
} while (0)
/**
 * @brief    Get extraction method from parameter list
 * @param    parameters     recipe parameter list
 * @param    context        read extraction method from this context
 * @return   newly allocated extraction method
 *
 * The parameter list should have been previously created using
 * fors_extract_define_parameters()
 */
/*----------------------------------------------------------------------------*/
extract_method *
fors_extract_method_new(                    const cpl_parameterlist *parameters,
                                            const char *context)
{
    extract_method *em = cpl_malloc(sizeof(*em));
    const char *name = NULL;
    const char *method = NULL;

    cpl_msg_info(cpl_func, "Extraction method:");

    cpl_msg_indent_more();
    //"sex" method is the default. The parameter won't appear when calling
    //the recipe from esorex, but it will in the unit tests fors_zeropoint-test,
    //fors_img_science-test.
    name = cpl_sprintf("%s.%s", context, "extract_method");
    if(cpl_parameterlist_find_const(parameters, name) == NULL)
        method = cpl_sprintf("%s", "sex");
    else
        method = cpl_sprintf("%s",dfs_get_parameter_string_const(parameters, name));
    cpl_free((void *)name); name = NULL;
    cpl_msg_indent_less();

    assure( !cpl_error_get_code(), return NULL, NULL );
    assure( method != NULL, return NULL, NULL );

    if (strcmp(method, "sex") == 0) {
        em->method = SEX;
        
        cpl_msg_indent_more();
        name = cpl_sprintf("%s.%s", context, "sex_exe");
        em->sex_exe = dfs_get_parameter_string_const(parameters, 
                                                     name);
        cpl_free((void *)name); name = NULL;
        cpl_msg_indent_less();

        
        cpl_msg_indent_more();
        name = cpl_sprintf("%s.%s", context, "sex_config");
        em->sex_config = dfs_get_parameter_string_const(parameters, 
                                                     name);
        cpl_free((void *)name); name = NULL;
        cpl_msg_indent_less();



        cpl_msg_indent_more();
        name = cpl_sprintf("%s.%s", context, "sex_mag");
        em->sex_mag = dfs_get_parameter_string_const(parameters, 
                                                     name);
        cpl_free((void *)name); name = NULL;
        cpl_msg_indent_less();


        cpl_msg_indent_more();
        name = cpl_sprintf("%s.%s", context, "sex_magerr");
        em->sex_magerr = dfs_get_parameter_string_const(parameters, 
                                                        name);
        cpl_free((void *)name); name = NULL;
        cpl_msg_indent_less();


        cpl_msg_indent_more();
        name = cpl_sprintf("%s.%s", context, "sex_radius");
        em->sex_radius = dfs_get_parameter_int_const(parameters, 
                                           name);
        cpl_free((void *)name); name = NULL;
        cpl_msg_indent_less();
    }
    else if (strcmp(method, "test") == 0) {
        em->method = TEST;
    }
    else {
        assure( false, return NULL, "Unknown extraction method '%s'", method);
    }

    cleanup;
    return em;
}

/*----------------------------------------------------------------------------*/
/**
 * @brief    Deallocate extraction method and set the pointer to NULL
 */
/*----------------------------------------------------------------------------*/
void
fors_extract_method_delete(                 extract_method **em)
{
    if (em && *em) {
        cpl_free(*em); *em = NULL;
    }
    return;
}

/*----------------------------------------------------------------------------*/
#undef cleanup
#define cleanup
/**
 * @brief   Extract sources
 * @param   image                  source image
 * @param   setting                instrument setting (gain)
 * @param   em                     extraction method to use
 * @param   sky_stats              (output) statistics on determined sky
 * @param   background             (output) inferred background image
 * @param   extracted_sources      (output) if non-NULL, table of extracted sources
 * @return  newly allocated list of extracted stars
 */
/*----------------------------------------------------------------------------*/
fors_star_list *
fors_extract(                               const fors_image *image, 
                                            const fors_setting *setting,
                                            const extract_method *em,
                                            double magsyserr,
                                            fors_extract_sky_stats *sky_stats,
                                            cpl_image **background,
                                            cpl_table **extracted_sources)
{
    assure( em != NULL, return NULL, NULL );

    cpl_msg_info(cpl_func, "Extracting sources");
    
    switch (em->method ) {
    case SEX: return extract_sex(image, setting,
                                 em->sex_exe,
                                 em->sex_config,
                                 em->sex_mag,
                                 em->sex_magerr,
                                 em->sex_radius,
                                 magsyserr,
                                 sky_stats,
                                 background,
                                 extracted_sources); break;
    case TEST: return extract_test(sky_stats, background, extracted_sources); break;
    default:
        assure( false, return NULL, "Unknown method %d", em->method );
        break;
    }
}

/*----------------------------------------------------------------------------*/
#undef cleanup
#define cleanup \
do { \
    cpl_table_delete(out); out = NULL; \
    cpl_free((void *)command); \
    cpl_image_delete(work_back); work_back = NULL; \
    cpl_image_delete(bmaxsigma); bmaxsigma = NULL; \
    cpl_image_delete(bsigma); bsigma = NULL; \
    fors_image_delete(&fbsigma); \
} while (0)
/**
 * @brief   Extract sources using SExtractor
 * @param   image           source image
 * @param   setting         instrument setting (gain)
 * @param   sex_exe         SExtractor executable
 * @param   sex_config      SExtractor configuration file
 * @param   sex_mag         SExtractor catalog magnitude
 * @param   sex_magerr      SExtractor catalog magnitude error
 * @param   radius          background error map median filter radius
 * @param   sky_stats       (output) statistics on determined sky
 * @param   background      (output) background image
 * @param   extracted_sources (output) if non-NULL, SExtractor output table
 * @return  newly allocated list of stars
 *
 * Note: The gain given in the setting must describe the
 * image. Therefore, if the provided value of the gain is just
 * the detector gain, the input image must not be stacked or 
 * normalized to e.g. 1s exposure time.
 *
 * A background error map is given to SExtractor. This error map is
 * obtained by applying a median filter to the input image error map 
 * (in order to remove sources).
 */
/*----------------------------------------------------------------------------*/
static fors_star_list *
extract_sex(                                const fors_image *image,
                                            const fors_setting *setting,
                                            const char *sex_exe,
                                            const char *sex_config,
                                            const char *sex_mag,
                                            const char *sex_magerr,
                                            int radius,
                                            double magsyserr,
                                            fors_extract_sky_stats *sky_stats,
                                            cpl_image **background,
                                            cpl_table **extracted_sources)
{
    const char *const filename_data  = "sextract_data.fits";
    const char *const filename_sigma = "sextract_bkg_sigma.fits";
    const char *const filename_cat   = "sextract_cat.fits";
    const char *const filename_bkg   = "sextract_bkg.fits";
    cpl_table *out = NULL;
    const char *command = NULL;
    fors_star_list *stars = NULL;
    cpl_image *work_back = NULL;
    cpl_image *bmaxsigma = NULL;
    cpl_image *bsigma = NULL;
    fors_image *fbsigma = NULL;
    fors_image *cropped = NULL;
    int croplx, croply, cropux, cropuy;
    const char *const filename_data_full  = "sextract_data_full.fits";
    const char *const filename_sigma_full = "sextract_bkg_sigma_full.fits";
    int        vignetted = 1;

    assure( setting != NULL, return NULL, NULL );

    assure( image != NULL, return NULL, NULL );

    assure( sky_stats != NULL, return NULL, NULL );
    assure( background != NULL, return NULL, NULL );


    /*
     * In case of new chips, crop
     */

    if (strcmp(setting->chip_id, "CCID20-14-5-6") == 0) {
        croplx =  380 / setting->binx;
        croply =  626 / setting->biny;
        cropux = 3714 / setting->binx;
        cropuy = 2048 / setting->biny;
    }
    else if (strcmp(setting->chip_id, "CCID20-14-5-3") == 0) {
        croplx =  380 / setting->binx;
        croply =    2 / setting->biny;
        cropux = 3714 / setting->binx;
        cropuy = 1920 / setting->biny;
    }
    else if (strncmp(setting->chip_id, "Marl", 4) == 0) {
        croplx =  380 / setting->binx;
        croply =  694 / setting->biny;
        cropux = 3690 / setting->binx;
        cropuy = 2048 / setting->biny;
    }
    else if (strncmp(setting->chip_id, "Norm", 4) == 0) {
        croplx =  380 / setting->binx;
        croply =    2 / setting->biny;
        cropux = 3690 / setting->binx;
        cropuy = 1894 / setting->biny;
    }
    else {
        vignetted = 0;
    }

    cropped = (fors_image *)image;       /* As a default.... */

    if (vignetted) {
       fors_image_save_sex(image, NULL,
                            filename_data_full,
                            filename_sigma_full,
                            radius);

        assure( !cpl_error_get_code(), return NULL,
                "Could not save image to %s and %s",
                filename_data_full, filename_sigma_full);

        cropped = fors_image_duplicate(image);
        fors_image_crop(cropped, croplx, croply, cropux, cropuy);
    }

    /* provide data and error bars in separate files,
       pass to sextractor */

    fors_image_save_sex(cropped, NULL,
                        filename_data,
                        filename_sigma,
                        radius);
    assure( !cpl_error_get_code(), return NULL,
            "Could not save image to %s and %s",
            filename_data, filename_sigma);

    if (vignetted)
        fors_image_delete(&cropped);

    
    /*
     * A = load filename_sigma
     *
     * M = fors_image_max_filter(A)
     *
     * |A-M|
     */

    command = cpl_sprintf("%s %s,%s "
                          "-c %s "
                          /* Use SExtractor's double mode which is probably
                             more tested and more bugfree than the
                             single image mode */
                          "-GAIN %f "   /* Different terminology here:
                                           SExtractor-gain == ESO-conad,
                                           unit = e- / ADU
                                        */
                          "-PIXEL_SCALE %f "
                          "-CHECKIMAGE_TYPE BACKGROUND "
                          "-CHECKIMAGE_NAME %s "
                          "-WEIGHT_TYPE MAP_RMS "
                          "-WEIGHT_IMAGE %s,%s "
                          "-CATALOG_TYPE FITS_1.0 "
                          "-CATALOG_NAME %s",
                          sex_exe,
                          filename_data, filename_data,
                          sex_config,
                          1.0/setting->average_gain,
                          setting->pixel_scale,
                          filename_bkg,
                          filename_sigma, filename_sigma,
                          filename_cat);

    cpl_msg_info(cpl_func, "Running '%s'", command);

    assure( system(command) == 0, return stars, "'%s' failed", command);

    /* 
     * The background map is here used just to evaluate a QC parameter,
     * and is also returned to the caller.
     */
    {
        int plane = 0;
        int extension = 0;

        *background = cpl_image_load(filename_bkg,
                                     CPL_TYPE_FLOAT, plane, extension);

        assure( !cpl_error_get_code(), return NULL,
                "Could not load SExtractor background image %s",
                filename_bkg );

        if (vignetted) {
            work_back = cpl_image_new(fors_image_get_size_x(image), 
                                      fors_image_get_size_y(image),
                                      CPL_TYPE_FLOAT);
            cpl_image_copy(work_back, *background, croplx, croply);

            assure( !cpl_error_get_code(), return NULL,
                    "Could not insert background image %s",
                    filename_bkg );

            cpl_image_delete(*background);
            *background = work_back;
            work_back = NULL;
        }
        
        /* 
         * To avoid using the non-illuminated parts, use only the central 
         * 100x100 window.
         *
         * It does not make too much sense to trend these QC parameters
         * anyway because they mostly describe the individual science 
         * exposures.
         */

        int nx = cpl_image_get_size_x(*background);
        int ny = cpl_image_get_size_y(*background);
        int xlo = nx/2 - 50;
        int xhi = nx/2 + 50;
        int ylo = ny/2 - 50;
        int yhi = ny/2 + 50;
        
        if (xlo <   0) xlo = 0;       /* Just in case... */
        if (xhi >= nx) xhi = nx - 1;
        if (ylo <   0) ylo = 0;
        if (yhi >= ny) yhi = ny - 1;

        work_back = cpl_image_duplicate(*background);
        
        sky_stats->mean   = cpl_image_get_mean_window(work_back, 
                                                      xlo, ylo, xhi, yhi);
        sky_stats->median = cpl_image_get_median_window(work_back, 
                                                 xlo, ylo, xhi, yhi);
        cpl_image_subtract_scalar(work_back, sky_stats->median);
        cpl_image_abs(work_back);
        sky_stats->rms    = cpl_image_get_median_window(work_back, 
                                                        xlo, ylo, xhi, yhi)
                          * STDEV_PR_MAD;

        cpl_image_delete(work_back); work_back = NULL;

        assure( !cpl_error_get_code(), return NULL,
                "Could not calculate sky statistics" );

    }
    
    cpl_msg_info(cpl_func, "Background = %f +- %f ADU",
                 sky_stats->median, sky_stats->rms);

    /* 
     * SExtractors background estimation is not reliable near
     * non-illuminated areas. The underestimated background 
     * leads to false detections.
     * Therefore, reject sources which are too close to the
     * illumination edge. The edge is masked using a max filter 
     * on the smoothed variance image of the background map. 
     * The filter must have comparable size to the SExtractor 
     * BACK_SIZE parameter. The discrimination level is half-way
     * between the high-variance and low-variance regions.
     */

    float level;
    {
        int xradius = 64;
        int yradius = 64;
        int plane = 0;
        int extension = 0;

        if (vignetted) {
            bsigma = cpl_image_load(filename_sigma_full,
                                    CPL_TYPE_FLOAT, plane, extension);
        }
        else {
            bsigma = cpl_image_load(filename_sigma,
                                    CPL_TYPE_FLOAT, plane, extension);
        }

        assure( !cpl_error_get_code(), return NULL,
                "Could not load SExtractor background error image %s",
                filename_sigma );

        /*
         * Duplication is necessary for creating the fors_image
         * to be passed to fors_image_filter_max_create()
         */
        
        /* this wraps the fors_image around the cpl images,
         * so set them to NULL */

        fbsigma = fors_image_new(cpl_image_duplicate(bsigma), bsigma);
        bsigma = NULL;

        {
            bool use_variance = true;
            bmaxsigma = fors_image_filter_max_create(fbsigma, xradius, 
                                                     yradius, use_variance);

     /* cpl_image_save(bmaxsigma, "/tmp/test.fits", CPL_BPP_IEEE_FLOAT, NULL,
                              CPL_IO_DEFAULT); */
            
        }

        fors_image_delete(&fbsigma);

        /*
          This is not robust if there are no non-illuminated areas.

          maxima = cpl_image_get_max(bmaxsigma);
          minima = cpl_image_get_min(bmaxsigma);
          level = (maxima + minima) / 2; 
        */
        
        /* 5 sigma rejection */
        level = cpl_image_get_median(bmaxsigma) * 5;

        cpl_msg_debug(cpl_func, "Threshold level = %f",
                      level);

    }

    out = cpl_table_load(filename_cat, 1, 1);

    assure( !cpl_error_get_code(), return NULL,
            "Could not load SExtractor output table %s",
            filename_cat); 

    /* Validate sextractor output */
    assure( cpl_table_has_column(out, "FLAGS"), return NULL,
            "%s: Missing column: %s", filename_cat, "FLAGS");

    assure( cpl_table_has_column(out, "CLASS_STAR"), return NULL,
            "%s: Missing column: %s", filename_cat, "CLASS_STAR");

    assure( cpl_table_has_column(out, "BACKGROUND"), return NULL,
            "%s: Missing column: %s", filename_cat, "BACKGROUND");

    assure( cpl_table_has_column(out, "X_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "X_IMAGE");

    assure( cpl_table_has_column(out, "Y_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "Y_IMAGE");

    assure( cpl_table_has_column(out, "FWHM_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "FWHM_IMAGE");

    assure( cpl_table_has_column(out, "A_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "A_IMAGE");

    assure( cpl_table_has_column(out, "B_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "B_IMAGE");

    assure( cpl_table_has_column(out, "THETA_IMAGE"), return NULL,
            "%s: Missing column: %s", filename_cat, "THETA_IMAGE");

    assure( cpl_table_has_column(out, sex_mag), return NULL,
            "%s: Missing column: %s", filename_cat, sex_mag);

    assure( cpl_table_has_column(out, sex_magerr), return NULL,
            "%s: Missing column: %s", filename_cat, sex_magerr);


    /* cpl_table_dump_structure(out, stdout); */

    if (vignetted) {
        cpl_table_add_scalar(out, "X_IMAGE", croplx - 1);
        cpl_table_add_scalar(out, "Y_IMAGE", croply - 1);
    }

    stars = fors_star_list_new();

    {
        int i;
        int bkg_rejected = 0;
        int rejected = 0;
        int nx = cpl_image_get_size_x(bmaxsigma);
        int ny = cpl_image_get_size_y(bmaxsigma);

        for (i = 0; i < cpl_table_get_nrow(out); i++) {
            fors_star       *s = NULL;
            unsigned int    flags = 0x0;
            int             x, y, tmp;
            double          bg_err;
            
            s = fors_star_new_from_table(   out,
                                            i,
                                            "X_IMAGE",
                                            "Y_IMAGE",
                                            "FWHM_IMAGE",
                                            "A_IMAGE",
                                            "B_IMAGE",
                                            "THETA_IMAGE", 
                                            sex_mag,
                                            sex_magerr,
                                            "CLASS_STAR");
            (*s).orientation *= M_PI/180;
            (*s).semi_major *= setting->binx;
            (*s).semi_minor *= setting->binx;
            (*s).fwhm *= setting->binx;
            (*s).dmagnitude = sqrt((*s).dmagnitude * (*s).dmagnitude + magsyserr * magsyserr);
            
            flags = cpl_table_get_int(      out, "FLAGS", i, NULL);
            
            x = (int)(s->pixel->x + 0.5);
            y = (int)(s->pixel->y + 0.5);
            if (x >= 1 && x <= nx && y >= 1 && y <= ny)
                bg_err = cpl_image_get(bmaxsigma, x, y, &tmp);
            else
                bg_err = -1.0;
            
            if (fors_extract_check_sex_flag(flags)
                && fors_extract_check_sex_star(s, bmaxsigma)
                && bg_err < level)
            {
                cpl_msg_debug(              cpl_func,
                                            "Source at (%f, %f): fwhm = %f px",
                                            s->pixel->x, s->pixel->y, s->fwhm);
                assure(                     !cpl_error_get_code(), return NULL,
                                            "Could not read SExtractor "
                                            "output table %s",
                                            filename_cat);
                fors_star_list_insert(stars, s);
            }
            else
            {
                cpl_msg_debug(              cpl_func,
                                            "Rejecting source at (%f, %f): "
                                            "flags = 0x%x; fwhm = %f pix; "
                                            "background error = %f; "
                                            "level = %f; "
                                            "background = %f",
                                            s->pixel->x, s->pixel->y,
                                            flags, s->fwhm,
                                            bg_err, level,
                                            cpl_table_get_float(
                                                out, "BACKGROUND", i, NULL));
                fors_star_delete(&s);
                rejected++;
                if (bg_err >= level)
                    bkg_rejected++;
            }
        }
        
        cpl_msg_info(cpl_func, "%d sources sextracted, %d rejected", 
                     fors_star_list_size(stars) + rejected,
                     rejected);
    }
    
    if (extracted_sources != NULL) {
        *extracted_sources = cpl_table_duplicate(out);
    }
    
    cleanup;
    return stars;
}

/*----------------------------------------------------------------------------*/
#undef cleanup
#define cleanup
/**
 * @brief   Extract sources
 * @param   sky_stats      (output) sky statistics
 * @param   background     (output) background
 * @param   extracted_sources (output) if non-NULL, table of sources
 * @return  newly allocated list of extracted stars
 *
 * The star positions are hardcoded!
 *
 * This method is used only for testing without relying on external
 * packages
 */
/*----------------------------------------------------------------------------*/
static fors_star_list *
extract_test(                               fors_extract_sky_stats *sky_stats,
                                            cpl_image **background,
                                            cpl_table **extracted_sources)
{
    assure( sky_stats != NULL, return NULL, NULL );
    assure( background != NULL, return NULL, NULL );
    
    sky_stats->mean = 1;
    sky_stats->median = 2;
    sky_stats->rms = 3;

    *background = cpl_image_new(10, 20, CPL_TYPE_FLOAT); /* Zero, wrong size */

    fors_star_list *stars = fors_star_list_new();

    struct {
        double x, y, magnitude, dmagnitude;
    } 
    data[] = {
        {100 ,  200, -10, 0.01},
        {1100,  200, -11, 0.01},
        {1   ,  5  , -10, 0.01},
        {100 , 1200, -12, 0.01},
        {1100, 1200, -13, 0.01}
    };
       
    int N = sizeof(data) / sizeof(*data);
    int i;
    double a = 2;
    double b = 1;
    double fwhm = 1.5;
    double orientation = 1.0; /* radians */

    for (i = 0; i < N; i++) {
        fors_star_list_insert(stars,
                              fors_star_new(data[i].x,
                                            data[i].y,
                                            fwhm,
                                            a, b, orientation,
                                            data[i].magnitude,
                                            data[i].dmagnitude,
                                            1.0));
    }

    if (extracted_sources != NULL) {
        *extracted_sources = fors_create_sources_table(stars);
        
        assure (!cpl_error_get_code(), return NULL, 
                "Could not create extracted sources table");
    }
    
    return stars;
}

/**@}*/