hdrl_persistence-test.c

/*
 * This file is part of the HDRL
 * Copyright (C) 2017 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
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
 
/*-----------------------------------------------------------------------------
                                    Includes
 -----------------------------------------------------------------------------*/
 
#include "hdrl.h"
#include <math.h>
#include <cpl.h>
#include "hdrl_persistence.c"
 
 
#if 0
static cpl_error_code get_error_and_reset(void){
    cpl_error_code r = cpl_error_get_code();
    cpl_error_reset();
    return r;
}
#endif
 
/*-----------------------------------------------------------------------------
                                    Define
 -----------------------------------------------------------------------------*/
//#define cpl_ensure_no_error cpl_test_eq(get_error_and_reset(), CPL_ERROR_NONE);
#define ARRAY_LEN(a) (sizeof((a))/sizeof((a)[0]))
 
#if 0
static void
test_example(void) {
 
    /* Below some example to verify results */
    cpl_frame* frm = cpl_frame_new();
    cpl_ensure_no_error; //one way to test no error
    cpl_test_error(CPL_ERROR_NONE); //another way to test no error
    cpl_test_nonnull(frm); // test created frame is not null.
 
    cpl_image_add_scalar(NULL, 1);
    cpl_test_error(CPL_ERROR_NULL_INPUT); //test expected error
 
    cpl_test_eq(1, 1);    //test results are as expected
    cpl_test_abs(1, 1, HDRL_EPS_DATA);//test results are as expected within accuracy
    cpl_test_rel(1, 1, HDRL_EPS_DATA);
 
 
    /* free memory */
    cpl_frame_delete(frm);
    /* verify no error is set by this unit test */
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
#endif
 
static hdrl_imagelist *
test_hdrl_persistence_mkhilist(
        const int                n,
        const cpl_size           nx,
        const cpl_size           ny,
        const bool               zero)
{
    hdrl_imagelist * list = hdrl_imagelist_new();
    for (int i=0; i<n; i++) {
        hdrl_image * hi = hdrl_image_new(nx, ny);
        cpl_image * im = hdrl_image_get_image(hi);
        if (!zero) cpl_image_fill_noise_uniform(im, -1, 1);
        hdrl_imagelist_set(list, hi, i);
    }
 
    return list;
}
 
static const char *
test_hdrl_persistence_make_trapMEF(
        const char * fn,
        int maxx,
        int maxy,
        cpl_type maxt,
        int denx,
        int deny,
        cpl_type dent,
        int fullwellx,
        int fullwelly,
        cpl_type fullwellt,
        int tabrows,
        const char * tabc1n,
        cpl_type tabc1t,
        const char * tabc2n,
        cpl_type tabc2t)
{
    cpl_image * max = cpl_image_new(maxx, maxy, maxt);
    cpl_image * den = cpl_image_new(denx, deny, dent);
    cpl_image * fullwell = cpl_image_new(fullwellx, fullwelly, fullwellt);
    cpl_table * frac = cpl_table_new(tabrows);
    cpl_table_new_column(frac, tabc1n, tabc1t);
    cpl_table_new_column(frac, tabc2n, tabc2t);
 
    cpl_propertylist * plist = cpl_propertylist_new();
 
    cpl_image_save(0, fn, CPL_TYPE_DOUBLE, 0, CPL_IO_CREATE);
 
    cpl_propertylist_update_string(plist, "EXTNAME", "PERSIST_TRAP_MAX");
    cpl_image_save(max, fn, cpl_image_get_type(max), plist, CPL_IO_EXTEND);
 
    cpl_propertylist_update_string(plist, "EXTNAME", "PERSIST_TRAP_DENSITY");
    cpl_image_save(den, fn, cpl_image_get_type(den), plist, CPL_IO_EXTEND);
 
    cpl_propertylist_update_string(plist, "EXTNAME", "PERSIST_FULL_WELL");
    cpl_image_save(fullwell, fn, cpl_image_get_type(fullwell), plist, CPL_IO_EXTEND);
 
    cpl_propertylist_update_string(plist, "EXTNAME", "PERSIST_TRAP_FRACTION");
    cpl_table_save(frac, 0, plist, fn, CPL_IO_EXTEND);
    cpl_propertylist_delete(plist);
    cpl_table_delete(frac);
    cpl_image_delete(den);
    cpl_image_delete(fullwell);
    cpl_image_delete(max);
    return fn;
}
 
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code test_hdrl_persistence_read_info(
        const char             * traps_,
        cpl_image             ** maximum,
        cpl_image             ** density,
        cpl_image             ** fullwell,
        cpl_table             ** frac)
{
    // maximum, density, fullwell, and frac checks: they should be non-NULL but pt to NULL
    cpl_error_ensure(
            maximum && !*maximum && density && !*density && fullwell && !*fullwell
            && frac && !*frac, CPL_ERROR_NULL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
                    " Input should be non-NULL but pointers to NULL");
 
 
    // traps_ checks
    cpl_size extmax = cpl_fits_find_extension(traps_, "PERSIST_TRAP_MAX");
    cpl_size extdens = cpl_fits_find_extension(traps_,"PERSIST_TRAP_DENSITY");
    cpl_size extfullwell = cpl_fits_find_extension(traps_,"PERSIST_FULL_WELL");
    cpl_size extfraction = cpl_fits_find_extension(traps_, "PERSIST_TRAP_FRACTION");
 
    cpl_error_ensure(!cpl_error_get_code(), CPL_ERROR_ILLEGAL_INPUT,
            return CPL_ERROR_DATA_NOT_FOUND,
                    "Can not find at least one of the required extensions");
 
    *maximum = cpl_image_load(traps_, CPL_TYPE_DOUBLE, 0, extmax);
    if(cpl_error_get_code() != CPL_ERROR_NONE) {
        return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
                "Couldn't load maximum trap map");
    }
 
    *density = cpl_image_load(traps_, CPL_TYPE_DOUBLE, 0, extdens);
    if(cpl_error_get_code() != CPL_ERROR_NONE) {
        cpl_image_delete(*maximum);
        *maximum = NULL;
        return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
                "Couldn't load trap density map");
    }
 
    *fullwell = cpl_image_load(traps_, CPL_TYPE_DOUBLE, 0, extfullwell);
    if(cpl_error_get_code() != CPL_ERROR_NONE) {
        cpl_image_delete(*maximum);
        cpl_image_delete(*density);
        *maximum = NULL;
        *density = NULL;
        return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
                "Couldn't load full well map");
    }
 
    *frac = cpl_table_load(traps_, extfraction, 0);
    if(cpl_error_get_code() != CPL_ERROR_NONE) {
        cpl_image_delete(*maximum);
        cpl_image_delete(*density);
        cpl_image_delete(*fullwell);
        *maximum = NULL;
        *density = NULL;
        *fullwell = NULL;
 
        return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
                "Couldn't load trap fraction table");
    }
    return CPL_ERROR_NONE;
}
 
 
 
 
 
static void
test_hdrl_persistence_verify_inputs(void)
{
    /* Test function:
     * static cpl_error_code hdrl_persistence_verify_inputs(
        const hdrl_parameter   * param_,
        const cpl_array        * dates_,
        const cpl_array        * exps_,
        const hdrl_imagelist   * illums_,
        const char             * traps_,
        const cpl_mask         * bpm_,
        const cpl_imagelist    * srcmasks_,
        cpl_image             ** maximum,
        cpl_image             ** density,
        cpl_table             ** frac);
 
        1) Verify that any case of invalid inputs (test each input) is properly
        handled (gracefully exit, set proper error code) and that error code is
        as expected
        2) Define some (dummy) inputs and verify the result error code is
        CPL_ERROR_NONE
     */
 
    /* LEGEND:
     *    bad1_*: NULL, or non-NULL for pointers to pointers to be populated
     *    bad2_*: too few
     *    bad3_*: too many
     *    bad4_*: wrong type
     *    bad5_*: wrong size (for images)
     */
 
    /* dates input variations */
    cpl_array * good_dates = cpl_array_new(2, CPL_TYPE_DOUBLE);
    cpl_array * bad1_dates = NULL;                               // NULL
    cpl_array * bad2_dates = cpl_array_new(0, CPL_TYPE_DOUBLE);  // too few
    cpl_array * bad3_dates = cpl_array_new(3, CPL_TYPE_DOUBLE);  // too many
    cpl_array * bad4_dates = cpl_array_new(2, CPL_TYPE_STRING);  // bad type
 
    /* exps input variations */
    cpl_array * good_exps = cpl_array_new(1, CPL_TYPE_DOUBLE);
    cpl_array * bad1_exps = NULL;                               // NULL
    cpl_array * bad2_exps = cpl_array_new(0, CPL_TYPE_DOUBLE);  // too few
    cpl_array * bad3_exps = cpl_array_new(2, CPL_TYPE_DOUBLE);  // too many
    cpl_array * bad4_exps = cpl_array_new(1, CPL_TYPE_STRING);  // bad type
 
    /* illums input variations */
    hdrl_imagelist * good_illums = test_hdrl_persistence_mkhilist(1, 3, 3, 0);
    hdrl_imagelist * bad1_illums = NULL;
    hdrl_imagelist * bad2_illums = hdrl_imagelist_new();
    hdrl_imagelist * bad3_illums = test_hdrl_persistence_mkhilist(2, 3, 3, 0);
    hdrl_imagelist * bad5_illums = test_hdrl_persistence_mkhilist(1, 4, 4, 0);
 
    /* traps input variations */
    const char * good_trap = test_hdrl_persistence_make_trapMEF("good_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
/*
    char * bad1_trap = NULL;
    char * bad2_trap = test_hdrl_persistence_make_trapMEF("bad2_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 5,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
    char * bad3_trap = test_hdrl_persistence_make_trapMEF("bad3_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 7,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
    char * bad4a_trap = test_hdrl_persistence_make_trapMEF("bad4a_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_DOUBLE, "NU", CPL_TYPE_FLOAT);
    char * bad4b_trap = test_hdrl_persistence_make_trapMEF("bad4b_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_DOUBLE);
    char * bad5a_trap = test_hdrl_persistence_make_trapMEF("bad5a_vi_trap.fits",
            3, 4, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
    char * bad5b_trap = test_hdrl_persistence_make_trapMEF("bad5b_vi_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 4, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
*/
 
    /* bpm input variations */
    cpl_mask * good_bpm = cpl_mask_new(3, 3);
    //cpl_mask * bad1_bpm = NULL;  // it's ok for bpm to be NULL
 
    /* srcmasks input variations */
    cpl_imagelist * good_srcmasks = cpl_imagelist_new();
    cpl_imagelist_set(good_srcmasks, cpl_image_new(3, 3, CPL_TYPE_INT), 0);
    //cpl_imagelist * bad1_srcmasks = NULL;  // it's ok for bpm to be NULL
    cpl_imagelist * bad2_srcmasks = cpl_imagelist_new();  // too few
    cpl_imagelist * bad3_srcmasks = cpl_imagelist_new();  // too many
    cpl_imagelist_set(bad3_srcmasks, cpl_image_new(3, 3, CPL_TYPE_INT), 0);
    cpl_imagelist_set(bad3_srcmasks, cpl_image_new(3, 3, CPL_TYPE_INT), 1);
    cpl_imagelist * bad4_srcmasks = cpl_imagelist_new();  // wrong type
    cpl_imagelist_set(bad4_srcmasks, cpl_image_new(3, 3, CPL_TYPE_FLOAT), 0);
 
    /* maximum trap map input variations */
    cpl_image * good_max = NULL;
    cpl_image * bad1_max = cpl_image_new(3, 3, CPL_TYPE_DOUBLE); // pre-exist
 
    /* density map input variations */
    cpl_image * good_den = NULL;
    cpl_image * bad1_den = cpl_image_new(3, 3, CPL_TYPE_DOUBLE); // pre-exist
 
    /* density map input variations */
    cpl_image * good_fullwell = NULL;
    cpl_image * bad1_fullwell = cpl_image_new(3, 3, CPL_TYPE_DOUBLE); // pre-exist
 
    /* fraction table input variations */
    cpl_table * good_frac = NULL;
    cpl_table * bad1_frac = cpl_table_new(3);  // pre-existing
 
    test_hdrl_persistence_read_info(good_trap, &good_max, &good_den, &good_fullwell, &good_frac);
 
     /* basic NULL tests (bad1_*) */
 
    hdrl_persistence_verify_inputs(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        10., 1000.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, bad1_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, bad1_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, bad1_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, NULL, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, NULL, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, NULL);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    /* non-NULL tests (bad1_*) */
/*
ToDo Now they should preexist - change this
    hdrl_persistence_verify_inputs(
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, bad1_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, bad1_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_verify_inputs(
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, bad1_frac);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
*/
 
    /* count requirement checks (bad2_*  bad3_*) */
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, bad2_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, bad3_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, bad2_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, bad3_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, bad2_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, bad3_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        bad2_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        bad3_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    /* type requirement checks (bad4_*) */
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, bad4_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INVALID_TYPE);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, bad4_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INVALID_TYPE);
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        bad4_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INVALID_TYPE);
 
    /* image size requirement checks (bad5_*) */
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, bad5_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
 
    /* everything is good check (this must come last) */
 
    hdrl_persistence_verify_inputs(1.0,
        1.0, 1.0, CPL_FALSE, good_dates, good_exps, good_illums,
        good_srcmasks, good_max, good_den, good_fullwell, good_frac);
    cpl_test_error(CPL_ERROR_NONE);
 
    /* cleanup */
 
    cpl_array_delete(good_dates);
    cpl_array_delete(bad2_dates);
    cpl_array_delete(bad3_dates);
    cpl_array_delete(bad4_dates);
    cpl_array_delete(bad2_exps);
    cpl_array_delete(bad3_exps);
    cpl_array_delete(bad4_exps);
    cpl_array_delete(good_exps);
    hdrl_imagelist_delete(good_illums);
    hdrl_imagelist_delete(bad2_illums);
    hdrl_imagelist_delete(bad3_illums);
    hdrl_imagelist_delete(bad5_illums);
    cpl_mask_delete(good_bpm);
    cpl_imagelist_delete(good_srcmasks);
    cpl_imagelist_delete(bad2_srcmasks);
    cpl_imagelist_delete(bad3_srcmasks);
    cpl_imagelist_delete(bad4_srcmasks);
    cpl_image_delete(good_max);
    cpl_image_delete(bad1_max);
    cpl_image_delete(good_den);
    cpl_image_delete(bad1_den);
    cpl_image_delete(good_fullwell);
    cpl_image_delete(bad1_fullwell);
    cpl_table_delete(good_frac);
    cpl_table_delete(bad1_frac);
 
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
 
static void
test_hdrl_persistence_threshhold_qi(void)
{
    /* Function to test:
     * static cpl_error_code hdrl_persistence_threshhold_img(
        hdrl_image            * hdrl_img,
        const double            lower,
        const cpl_image       * upper);
 
        1) Test invalid input and verify are handled properly and error code
           is proper
        2) Generate a hdrl_image with an image 10x10 pix of intensity
           increasing at each pixel from 0 to 100. Set as lower threshold a
           10x10 image set to 10, and as upper threshold a 10x10 image set to
           90.  Verify that the expected pixel position have proper value 10 and
           the other expected pixel positions have value 90 (it may be useful
           to save the image and check it wit ds9)
     */
    const cpl_size ny = 10;
    const cpl_size nx = 10;
 
    const hdrl_image * zero = hdrl_image_new(nx, ny);
    cpl_image * tmp = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
    cpl_image_fill_window(tmp, 1, 1, nx, ny, 90.0);
    const hdrl_image * upper = hdrl_image_create(tmp, NULL);
    cpl_image_fill_window(tmp, 1, 1, nx, ny, 10.0);
    const hdrl_image * lower = hdrl_image_create(tmp, NULL);
 
    hdrl_image * hi = hdrl_image_new(nx, ny);
    for (int i=0, xpos=1; xpos <= nx; xpos++) {
        for (int ypos=1 ; ypos <= ny; ypos++ ) {
            hdrl_value val = { i++, 0.0 };
            hdrl_image_set_pixel(hi, xpos, ypos, val);
        }
    }
 
    hdrl_persistence_threshhold_img(NULL, zero, upper);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_threshhold_img(hi, NULL, upper);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_threshhold_img(NULL, NULL, NULL);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_threshhold_img(hi, zero, upper);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 1, 1, NULL).data, 0.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 1, 2, NULL).data, 1.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 10, 9, NULL).data, 90.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 10, 10, NULL).data, 90.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    hdrl_persistence_threshhold_img(hi, lower, hi);  // uses itself as upper lim
 
    cpl_test_abs(  // 0.0 -> 10.0
        hdrl_image_get_pixel(hi, 1, 1, NULL).data, 10.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(  // 1.0 -> 10.0
        hdrl_image_get_pixel(hi, 1, 2, NULL).data, 10.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(  // last pixel val should remain unchanged
        hdrl_image_get_pixel(hi, 10, 10, NULL).data, 90.0, HDRL_EPS_DATA);
    cpl_test_error(CPL_ERROR_NONE);
 
    hdrl_image_delete(hi);
    hdrl_image_delete((hdrl_image *)zero);
    hdrl_image_delete((hdrl_image *)upper);
    hdrl_image_delete((hdrl_image *)lower);
    cpl_image_delete(tmp);
 
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
 
static void
test_hdrl_persistence_zero_flagged_pixels(void)
{
    /*Function to test:
     * static cpl_error_code hdrl_persistence_zero_flagged_pixels(
        hdrl_image            * hdrl_img,
        const cpl_mask        * mask,
        const cpl_binary        value);
 
        1) Test invalid input and verify are handled properly and error code
           is proper
        2) Simulate hdrl image and mask 3x3 pixels. Mask certain known
           position.  Verify that after operation hdrl image values have proper
           value either if masked or not.
     */
 
    const cpl_size ny = 3;
    const cpl_size nx = 3;
    cpl_image * tmp = cpl_image_new(ny, nx, CPL_TYPE_DOUBLE);
    cpl_image_add_scalar(tmp, 7.0);
    hdrl_image * hi = hdrl_image_create(tmp, NULL);
 
    cpl_mask * srcmask = cpl_mask_new(nx, ny);
    cpl_mask_set(srcmask, 2, 2, CPL_BINARY_1);
 
    hdrl_persistence_zero_flagged_pixels(NULL, srcmask, CPL_BINARY_1);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_zero_flagged_pixels(hi, NULL, CPL_BINARY_1);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_zero_flagged_pixels(hi, srcmask, 's');
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
 
    hdrl_persistence_zero_flagged_pixels(hi, srcmask, CPL_BINARY_1);
    cpl_test_error(CPL_ERROR_NONE);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 2, 2, NULL).data, 0.0, HDRL_EPS_DATA);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 1, 2, NULL).data, 7.0, HDRL_EPS_DATA);
 
    cpl_test_abs(
        hdrl_image_get_pixel(hi, 2, 3, NULL).data, 7.0, HDRL_EPS_DATA);
 
    hdrl_image_delete(hi);
    cpl_mask_delete(srcmask);
    cpl_image_delete(tmp);
 
    return;
}
 
static void
test_hdrl_persistence_compute_qi(void)
{
    /* Function to test:
     * static cpl_error_code hdrl_persistence_compute_qi(
        hdrl_imagelist         * Q,
        hdrl_imagelist         * Qacc,
        const double             delt,
        const hdrl_imagelist   * hirhos,
        const float            * taud,
        const double             exptime,
        const hdrl_image       * hiillum,
        const int                n_tau,
        const hdrl_imagelist   * himaxs)
 
        1) Verify invalid input
        2) Generate valid input and verify results are proper (inspect a few
           pixel values of the result)
     */
 
    const float taud[] = {1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0};
    const int n_tau = 6;
 
    const cpl_size ny = 20;
    const cpl_size nx = 20;
    hdrl_imagelist * hirhos = test_hdrl_persistence_mkhilist(n_tau, nx, ny, 0);
    hdrl_imagelist * himaxs = test_hdrl_persistence_mkhilist(n_tau, nx, ny, 0);
    hdrl_imagelist * Q = test_hdrl_persistence_mkhilist(n_tau, nx, ny, 1);
    hdrl_imagelist * Qacc = test_hdrl_persistence_mkhilist(n_tau, nx, ny, 1);
 
    /* Create an image with all pixels set to 1.0 and a single sharp Gaussian
     * peak in the image at coordinate (100, 100). */
    cpl_image * peak_img = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
    cpl_test_error(CPL_ERROR_NONE);
    cpl_image_fill_gaussian(peak_img, 10, 10, 1000.0, 1.0, 0.5);
    cpl_test_error(CPL_ERROR_NONE);
    cpl_image_add_scalar(peak_img, 1.0);
    cpl_test_error(CPL_ERROR_NONE);
    const hdrl_image * hiillum = hdrl_image_create(peak_img, NULL);
 
    double exptime = 300.0;
    double diff = 58967.353850 - 58967.315219;
    diff *= 24.0 * 3600.0;
 
    hdrl_persistence_compute_qi(
        NULL, Qacc, diff, hirhos, taud, exptime, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, NULL, diff, hirhos, taud, exptime, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, NULL, taud, exptime, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, NULL, exptime, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, taud, -1, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, taud, exptime, NULL, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, taud, exptime, hiillum, -1, himaxs);
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, taud, exptime, hiillum, n_tau, NULL);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
 
    hdrl_persistence_compute_qi(
        Q, Qacc, diff, hirhos, taud, exptime, hiillum, n_tau, himaxs);
    cpl_test_error(CPL_ERROR_NONE);
 
    // now we compare the hdrl_persistence_compute_qi() results to
    // manually-computed values for some randomly-chosen pixels
    typedef struct { int x; int y; } pixel_t;
    pixel_t pixs[] = { {4,6}, {9,2} };  // some random pixels to test
    for (size_t j = 0; j < ARRAY_LEN(pixs); ++j) {
        const pixel_t p = pixs[j];
 
        // first declare & initialise required inputs
        hdrl_data_t illum = hdrl_image_get_pixel(hiillum, p.x, p.y, NULL).data;
        hdrl_data_t q[] = { 0, 0, 0, 0, 0, 0 };     // starts as all zeros
        hdrl_data_t qacc[] = { 0, 0, 0, 0, 0, 0 };  // starts as all zeros
        hdrl_data_t rhos[] = { 0, 0, 0, 0, 0, 0 };  // next loop populates
        hdrl_data_t maxs[] = { 0, 0, 0, 0, 0, 0 };  // next loop populates
        for (size_t i = 0; i < ARRAY_LEN(rhos); ++i) {
            hdrl_image * hi = hdrl_imagelist_get(hirhos, i);
            rhos[i] = hdrl_image_get_pixel(hi, p.x, p.y, NULL).data;
            hi = hdrl_imagelist_get(himaxs, i);
            maxs[i] = hdrl_image_get_pixel(hi, p.x, p.y, NULL).data;
        }
 
        // now perform the calculations & comparisons
        for (int i=0; i<n_tau; i++) {
            // compute Qi
            hdrl_data_t tmp = qacc[i] * exp(-diff/taud[i]);
            tmp -= q[i];
            tmp *= diff/taud[i];
            q[i] += tmp;
            // compare Qi
            hdrl_image * hi = hdrl_imagelist_get(Q, i);
            hdrl_data_t val = hdrl_image_get_pixel(hi, p.x, p.y, NULL).data;
            //cpl_test_abs(q[i], val, HDRL_EPS_DATA);  // not useful: always 0
                                                       // vs 0 after 1st frame
                                                       // as this test simulates
 
            // compute Qacci
            tmp = illum * rhos[i] * (1.0 - exp(-exptime/taud[i]));
            if (tmp > maxs[i]) tmp = maxs[i];  // threshhold
            qacc[i] = tmp;
            // compare Qacci
            hi = hdrl_imagelist_get(Qacc, i);
            val = hdrl_image_get_pixel(hi, p.x, p.y, NULL).data;
            cpl_test_abs(qacc[i], val, HDRL_EPS_DATA);
        }
    }
 
    hdrl_image_delete((hdrl_image *)hiillum);
    cpl_image_delete(peak_img);
    hdrl_imagelist_delete(Qacc);
    hdrl_imagelist_delete(Q);
    hdrl_imagelist_delete(himaxs);
    hdrl_imagelist_delete(hirhos);
 
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
 
static void
test_hdrl_persistence_calc_stats(void)
{
 
    /* Function to test:
     * static cpl_propertylist * hdrl_persistence_calc_stats(
        const hdrl_image * Qtot,
        const double trim_perc)
 
        1) Test invalid input
        2) Generate valid input and verify statistics are as expected.
     */
 
    //cpl_image * im = cpl_image_create();
    //cpl_propertylist * pl = hdrl_persistence_calc_stats(0, 0.1);
    //cpl_propertylist_delete(pl);
 
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
 
static void
test_hdrl_persistence_compute(void)
{
    /* Function to test:
       hdrl_persistence_result * hdrl_persistence_compute(
        const hdrl_parameter   * param_,
        const cpl_array        * dates_,
        const cpl_array        * exps_,
        const hdrl_imagelist   * illums_,
        const char             * traps_,
        const cpl_mask         * bpm_,
        const cpl_imagelist    * srcmasks_)
 
        1) Verify invalid input
        2) generate valid inputs and verify results are as expected
        (for example that parameters have the expected value)
     */
 
    double                  det_turnover = 1.0;
    const double            mean_trim = 1.0;
    const cpl_boolean       cleanQ = CPL_FALSE;
    const cpl_array * dates = cpl_array_new(2, CPL_TYPE_FLOAT);
    const cpl_array * exps = cpl_array_new(1, CPL_TYPE_FLOAT);
    hdrl_imagelist * illums = test_hdrl_persistence_mkhilist(1, 3, 3, 0);
    const char * traps = test_hdrl_persistence_make_trapMEF("compute_trap.fits",
            3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 3, 3, CPL_TYPE_DOUBLE, 6,
            "TAU", CPL_TYPE_FLOAT, "NU", CPL_TYPE_FLOAT);
    const cpl_mask * bpm = cpl_mask_new(3, 3);
    cpl_imagelist * srcmasks = cpl_imagelist_new();
    cpl_imagelist_set(srcmasks, cpl_image_new(3, 3, CPL_TYPE_INT), 0);
 
    // no point testing NULL or other invalid inputs here as that is already
    // performed by test_hdrl_persistence_verify_inputs (hdrl_persistence_compute
    // calls hdrl_persistence_verify_inputs as its 1st step)
 
    det_turnover = -45727230.0;
 
 
    cpl_table * frac = NULL;
    cpl_image * maximum = NULL, * density = NULL, * fullwell = NULL;
 
    test_hdrl_persistence_read_info(traps, &maximum, &density, &fullwell, &frac);
 
    hdrl_image * persistence = NULL;
    cpl_propertylist * persistence_qc = NULL;
    double gain = 1.0;
    hdrl_persistence_compute(gain, det_turnover,
            mean_trim, cleanQ, dates,
            exps, illums, srcmasks, maximum, density, fullwell, frac, &persistence, &persistence_qc);
    /*ToDo Fix later*/
    cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
 
    hdrl_image_delete(persistence);
    cpl_propertylist_delete(persistence_qc);
    cpl_image_delete(maximum);
    cpl_image_delete(density);
    cpl_image_delete(fullwell);
    cpl_table_delete(frac);
 
    cpl_imagelist_delete(srcmasks);
    cpl_mask_delete((cpl_mask *)bpm);
    hdrl_imagelist_delete(illums);
    cpl_array_delete((cpl_array *)exps);
    cpl_array_delete((cpl_array *)dates);
 
    cpl_test_error(CPL_ERROR_NONE);
    return;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int main(void)
{
    cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
    //test_example();
 
    test_hdrl_persistence_verify_inputs();
    //test_hdrl_persistence_parameter_verify();
    test_hdrl_persistence_threshhold_qi();
    test_hdrl_persistence_zero_flagged_pixels();
    test_hdrl_persistence_compute_qi();
    test_hdrl_persistence_calc_stats();
    //test_hdrl_persistence_parameter_create();
    //test_hdrl_persistence_parameter_check();
    test_hdrl_persistence_compute();
    //cpl_test_error(CPL_ERROR_NONE);
 
    return cpl_test_end(0);
}