CR2RE Pipeline Reference Manual 1.7.0
hdrl_cat_sim_montecarlo-test.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2017 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include "config.h"
22#endif
23
24#include <cpl_test.h>
25
26#include "../hdrl_cat_casu.h"
27
28#include "hdrl_random.h"
29
30
31#define COMP_TOL_REL 1. / 3.
32#define COMP_TOL_ABS 1e-2
33#define IMG_XSIZE 120
34#define IMG_YSIZE 180
35
36
37void test_basic(void)
38{
39 /* Generate a field with some stars and a confidence map */
40 cpl_image *bkg = cpl_image_new(IMG_XSIZE, IMG_YSIZE, CPL_TYPE_DOUBLE);
41 cpl_image *im = cpl_image_new(IMG_XSIZE, IMG_YSIZE, CPL_TYPE_DOUBLE);
42 cpl_image *cnf = cpl_image_new(IMG_XSIZE, IMG_YSIZE, CPL_TYPE_DOUBLE);
43
44 /* Generate confidence uniform confidence map */
45 cpl_image_add_scalar(cnf, 100.);
46
47 /* Generate sky and a gaussian source without errors */
48 double sigma = 2.;
49 double norm2 = 2. * CPL_MATH_PI * sigma * sigma;
50
51 double sky = 500.;
52 cpl_image_add_scalar(bkg, sky);
53
54 double xpos = 80.;
55 double ypos = 100.;
56 double norm = 3000.;
57 cpl_image_fill_gaussian(im, xpos, ypos, norm * norm2, sigma, sigma);
58 cpl_image_add(bkg, im);
59
60 hdrl_casu_fits *inf = hdrl_casu_fits_wrap( im);
61 hdrl_casu_fits *inconf = hdrl_casu_fits_wrap(cnf);
62
63 /* Give it a WCS */
64 cpl_propertylist *pl = hdrl_casu_fits_get_ehu(inf);
65 cpl_propertylist_update_string(pl, "CTYPE1", "RA---TAN" );
66 cpl_propertylist_update_string(pl, "CTYPE2", "DEC--TAN" );
67 cpl_propertylist_update_double(pl, "CRVAL1", 30. );
68 cpl_propertylist_update_double(pl, "CRVAL2", 12. );
69 cpl_propertylist_update_double(pl, "CRPIX1", 512. );
70 cpl_propertylist_update_double(pl, "CRPIX2", 512. );
71 cpl_propertylist_update_double(pl, "CD1_1", -1. / 3600.);
72 cpl_propertylist_update_double(pl, "CD1_2", 0. );
73 cpl_propertylist_update_double(pl, "CD2_1", 0. );
74 cpl_propertylist_update_double(pl, "CD2_2", 1. / 3600.);
75 cpl_propertylist_update_int( pl, "NAXIS1", IMG_XSIZE );
76 cpl_propertylist_update_int( pl, "NAXIS2", IMG_YSIZE );
77 cpl_wcs *wcs = cpl_wcs_new_from_propertylist(pl);
78
79 /* Perform the montecarlo simulation */
80 hdrl_casu_result *res = cpl_calloc(1, sizeof(hdrl_casu_result));
81 cpl_table *tabfinal = NULL;
82 for (cpl_size iterate = 0; iterate < 100; iterate++) {
83
84 /* Add realistic poisonnian noise */
85 hdrl_random_state * rng = hdrl_random_state_new(1, NULL);
86 cpl_size size = cpl_image_get_size_x(bkg) * cpl_image_get_size_y(bkg);
87 double *pim = cpl_image_get_data_double(im);
88 double *pbkg= cpl_image_get_data_double(bkg);
89
90 for (cpl_size i = 0; i < size; i++) {
91 pim[i] = hdrl_random_poisson(rng, pbkg[i]);
92 }
93 hdrl_random_state_delete(rng);
94
95 /* Run casu catalogue */
96 hdrl_casu_catalogue(inf, inconf, wcs, 5, 2.5, 0, 3., 1, 32,
97 HDRL_CATALOGUE_ALL, 3.0, 1.0, HDRL_SATURATION_INIT, res);
98
99 cpl_table *tab = hdrl_casu_tfits_get_table(res->catalogue);
100
101 if (iterate == 0) {
102 tabfinal = cpl_table_duplicate(tab);
103 } else {
104 cpl_table_insert(tabfinal, tab, iterate);
105 }
106
107 cpl_image_delete( res->segmentation_map);
108 cpl_image_delete( res->background);
109 hdrl_casu_tfits_delete(res->catalogue);
110 }
111
112 /* Do the checks */
113 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "X_coordinate_err"),
114 cpl_table_get_column_stdev(tabfinal, "X_coordinate" ), COMP_TOL_REL);
115 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Y_coordinate_err"),
116 cpl_table_get_column_stdev(tabfinal, "Y_coordinate" ), COMP_TOL_REL);
117 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Peak_height_err" ),
118 cpl_table_get_column_stdev(tabfinal, "Peak_height" ), COMP_TOL_REL);
119 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_1_err" ),
120 cpl_table_get_column_stdev(tabfinal, "Aper_flux_1" ), COMP_TOL_REL);
121 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_2_err" ),
122 cpl_table_get_column_stdev(tabfinal, "Aper_flux_2" ), COMP_TOL_REL);
123 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_3_err" ),
124 cpl_table_get_column_stdev(tabfinal, "Aper_flux_3" ), COMP_TOL_REL);
125 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_4_err" ),
126 cpl_table_get_column_stdev(tabfinal, "Aper_flux_4" ), COMP_TOL_REL);
127 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_5_err" ),
128 cpl_table_get_column_stdev(tabfinal, "Aper_flux_5" ), COMP_TOL_REL);
129 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_6_err" ),
130 cpl_table_get_column_stdev(tabfinal, "Aper_flux_6" ), COMP_TOL_REL);
131 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_7_err" ),
132 cpl_table_get_column_stdev(tabfinal, "Aper_flux_7" ), COMP_TOL_REL);
133 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_8_err" ),
134 cpl_table_get_column_stdev(tabfinal, "Aper_flux_8" ), COMP_TOL_REL);
135 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_9_err" ),
136 cpl_table_get_column_stdev(tabfinal, "Aper_flux_9" ), COMP_TOL_REL);
137 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_10_err"),
138 cpl_table_get_column_stdev(tabfinal, "Aper_flux_10" ), COMP_TOL_REL);
139 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_11_err"),
140 cpl_table_get_column_stdev(tabfinal, "Aper_flux_11" ), COMP_TOL_REL);
141 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_12_err"),
142 cpl_table_get_column_stdev(tabfinal, "Aper_flux_12" ), COMP_TOL_REL);
143 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Aper_flux_13_err"),
144 cpl_table_get_column_stdev(tabfinal, "Aper_flux_13" ), COMP_TOL_REL);
145 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Petr_flux_err" ),
146 cpl_table_get_column_stdev(tabfinal, "Petr_flux" ), COMP_TOL_REL);
147 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Kron_flux_err" ),
148 cpl_table_get_column_stdev(tabfinal, "Kron_flux" ), COMP_TOL_REL);
149 cpl_test_rel(cpl_table_get_column_mean( tabfinal, "Half_flux_err" ),
150 cpl_table_get_column_stdev(tabfinal, "Half_flux" ), COMP_TOL_REL);
151
152 /* Clean up */
153 hdrl_casu_fits_delete(inf);
154 hdrl_casu_fits_delete(inconf);
155 cpl_wcs_delete(wcs);
156 cpl_table_delete(tabfinal);
157 cpl_image_delete(bkg);
158 cpl_free(res);
159}
160
161/*----------------------------------------------------------------------------*/
165/*----------------------------------------------------------------------------*/
166int main(void)
167{
168 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
169
170 test_basic();
171
172 return cpl_test_end(0);
173}
hdrl_casu_fits_wrap
hdrl_casu_fits * hdrl_casu_fits_wrap(cpl_image *im)
Wrap an image in a hdrl_casu_fits wrapper.
Definition hdrl_cat_casu.c:165
hdrl_casu_tfits_delete
void hdrl_casu_tfits_delete(hdrl_casu_tfits *p)
Free all the workspace associated with a hdrl_casu_fits object.
Definition hdrl_cat_casu.c:139
hdrl_casu_fits_delete
void hdrl_casu_fits_delete(hdrl_casu_fits *p)
Free all the workspace associated with a hdrl_casu_fits object.
Definition hdrl_cat_casu.c:267
hdrl_casu_tfits_get_table
cpl_table * hdrl_casu_tfits_get_table(hdrl_casu_tfits *p)
Get the CPL table from the hdrl_casu_tfits object.
Definition hdrl_cat_casu.c:121
hdrl_casu_fits_get_ehu
cpl_propertylist * hdrl_casu_fits_get_ehu(hdrl_casu_fits *p)
Get the propertylist for the extension header for a given hdrl_casu_fits image.
Definition hdrl_cat_casu.c:226
hdrl_casu_catalogue
cpl_error_code hdrl_casu_catalogue(hdrl_casu_fits *infile, hdrl_casu_fits *conf, const cpl_wcs *wcs, cpl_size ipix, double threshold, cpl_size icrowd, double rcore, cpl_size bkg_subtr, cpl_size nbsize, hdrl_catalogue_options cattype, double filtfwhm, double gainloc, double saturation, hdrl_casu_result *res)
Generate object catalogues from input images.
Definition hdrl_cat_casu.c:321
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2