CR2RE Pipeline Reference Manual 1.7.0
hdrl_cat_overlp-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_apio.h"
27#include "../hdrl_cat_areals.h"
28#include "../hdrl_cat_conf.h"
29#include "../hdrl_cat_moments.h"
30#include "../hdrl_cat_overlp.h"
31
32
33#define NTEST 10
34#define NT 117
35#define SIZE_IMAGE 512
36#define CENTER_IMAGE 256
37
38
39void standard_test(void)
40{
41 /* Intitialize */
42 double xpos[] = { 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000.};
43 double ypos[] = { 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000.};
44 double norm[] = {1000., 100., 200., 500., 550., 600., 650., 700., 750., 800.};
45
46 /* Create input data */
47 cpl_image *im = cpl_image_new(1024, 1024, CPL_TYPE_DOUBLE);
48 cpl_image *bkg = cpl_image_new(1024, 1024, CPL_TYPE_DOUBLE);
49 cpl_image *cnf = cpl_image_new(1024, 1024, CPL_TYPE_DOUBLE);
50
51 cpl_image_fill_noise_uniform(bkg, -10., 10. );
52 cpl_image_fill_noise_uniform(cnf, 99.9, 100.1);
53
54 hdrl_casu_fits *infErr = hdrl_casu_fits_wrap( im);
55 hdrl_casu_fits *inconf = hdrl_casu_fits_wrap(cnf);
56 hdrl_casu_fits *inf = hdrl_casu_fits_wrap(bkg);
57
58 double sigma = 2.;
59 double norm2 = 2. * CPL_MATH_PI * sigma * sigma;
60
61 double sky = 500.;
62 cpl_image_add_scalar(bkg, sky);
63
64 for (cpl_size i = 0; i < NTEST; i++) {
65 cpl_image_fill_gaussian(im, xpos[i], ypos[i], norm[i] * norm2, sigma, sigma);
66 cpl_image_add(bkg, im);
67 }
68
69 /*** Tests ***/
70 hdrl_casu_result *res = cpl_calloc(1, sizeof(hdrl_casu_result));
71 cpl_size ipix;
72 cpl_size icrowd;
73
74 /* Test 1: Error input null */
75 ipix = 5;
76 icrowd = 0;
77 cpl_test_eq(hdrl_catalogue_conf(infErr, NULL, ipix, 2.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
78 cpl_test_nonnull( res->catalogue);
79 hdrl_casu_tfits_delete(res->catalogue);
80 cpl_image_delete( res->segmentation_map);
81 cpl_image_delete( res->background);
82
83 /* Test 2 */
84 ipix = 5;
85 icrowd = 0;
86 cpl_test_eq(hdrl_catalogue_conf(infErr, inconf, ipix, 2.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
87 cpl_test_nonnull( res->catalogue);
88 hdrl_casu_tfits_delete(res->catalogue);
89 cpl_image_delete( res->segmentation_map);
90 cpl_image_delete( res->background);
91
92 /* Test 3 */
93 ipix = 10;
94 icrowd = 10;
95 cpl_test_eq(hdrl_catalogue_conf(infErr, NULL, ipix, 2.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
96 cpl_test_nonnull( res->catalogue);
97 hdrl_casu_tfits_delete(res->catalogue);
98 cpl_image_delete( res->segmentation_map);
99 cpl_image_delete( res->background);
100
101 /* Test 4 */
102 ipix = 10;
103 icrowd = 10;
104 cpl_test_eq(hdrl_catalogue_conf(infErr, inconf, ipix, 2.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
105 cpl_test_nonnull( res->catalogue);
106 hdrl_casu_tfits_delete(res->catalogue);
107 cpl_image_delete( res->segmentation_map);
108 cpl_image_delete( res->background);
109
110 /* Test 5 */
111 ipix = 0;
112 icrowd = 5;
113 cpl_test_eq(hdrl_catalogue_conf(inf, NULL, ipix, 1.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
114 cpl_test_nonnull( res->catalogue);
115 hdrl_casu_tfits_delete(res->catalogue);
116 cpl_image_delete( res->segmentation_map);
117 cpl_image_delete( res->background);
118
119 /* Test 6 */
120 ipix = 0;
121 icrowd = 5;
122 cpl_test_eq(hdrl_catalogue_conf(inf, inconf, ipix, 1.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
123 cpl_test_nonnull( res->catalogue);
124 hdrl_casu_tfits_delete(res->catalogue);
125 cpl_image_delete( res->segmentation_map);
126 cpl_image_delete( res->background);
127
128 /* Test 7 */
129 ipix = 10;
130 icrowd = 10;
131 cpl_test_eq(hdrl_catalogue_conf(inf, NULL, ipix, 1.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
132 cpl_test_nonnull( res->catalogue);
133 hdrl_casu_tfits_delete(res->catalogue);
134 cpl_image_delete( res->segmentation_map);
135 cpl_image_delete( res->background);
136
137 /* Test 8 */
138 ipix = 10;
139 icrowd = 10;
140 cpl_test_eq(hdrl_catalogue_conf(inf, inconf, ipix, 1.5, icrowd, 5., 1, 64, 6, 3., 1., INFINITY, res), CPL_ERROR_NONE);
141 cpl_test_nonnull( res->catalogue);
142 hdrl_casu_tfits_delete(res->catalogue);
143 cpl_image_delete( res->segmentation_map);
144 cpl_image_delete( res->background);
145
146 /* Clean up */
147 hdrl_casu_fits_delete(infErr);
148 hdrl_casu_fits_delete(inconf);
150 cpl_free(res);
151}
152
153void advanced_test(void)
154{
155 cpl_size x[] = {398,399,400,397,398,399,400,401,402,403,396,397,398,399,400,
156 401,402,403,404,395,396,397,398,399,400,401,402,403,404,405,
157 395,396,397,398,399,400,401,402,403,404,405,395,396,397,398,
158 399,400,401,402,403,404,405,395,396,397,398,399,400,401,402,
159 403,404,405,406,395,396,397,398,399,400,401,402,403,404,405,
160 395,396,397,398,399,400,401,402,403,404,405,395,396,397,398,
161 399,400,401,402,403,404,405,396,397,398,399,400,401,402,403,
162 404,397,398,399,400,401,402,403,398,399,400,401};
163
164 cpl_size y[] = {394,394,394,395,395,395,395,395,395,395,396,396,396,396,396,
165 396,396,396,396,397,397,397,397,397,397,397,397,397,397,397,
166 398,398,398,398,398,398,398,398,398,398,398,399,399,399,399,
167 399,399,399,399,399,399,399,400,400,400,400,400,400,400,400,
168 400,400,400,400,401,401,401,401,401,401,401,401,401,401,401,
169 402,402,402,402,402,402,402,402,402,402,402,403,403,403,403,
170 403,403,403,403,403,403,403,404,404,404,404,404,404,404,404,
171 404,405,405,405,405,405,405,405,406,406,406,406};
172
173 double z[] = {8.87152,12.515,7.69699,10.8527,22.2509,21.7368,13.0388,
174 12.1853,17.1976,7.43948,15.2245,29.1964,37.9117,57.9371,
175 71.5542,57.1288,34.7726,15.5934,11.5374,15.995,21.3606,
176 60.4006,103.46,147.55,168.274,147.476,98.9157,51.7186,20.188,
177 3.04248,5.77832,49.3103,98.2057,187.557,268.353,310.638,
178 274.295,183.969,94.6933,47.9889,20.245,26.3758,59.1781,
179 152.389,275.916,395.107,450.251,397.53,272.322,147.053,54.767,
180 11.8971,13.3888,73.3689,165.899,298.455,449.707,493.25,441.585,
181 299.31,157.474,70.1224,15.5313,8.76074,20.7188,54.5798,141.249,
182 264.87,382.736,435.452,393.871,268.175,138.485,65.9307,28.7812,
183 19.379,36.6449,93.5458,186.823,270.95,305.093,260.879,183.683,
184 100.676,32.6281,16.6497,5.94965,17.8105,57.256,106.32,145.264,
185 164.271,137.093,88.9384,60.7841,31.8582,10.0435,4.69162,
186 15.2187,32.5385,61.0381,74.5399,67.3727,43.3964,25.0956,
187 16.7595,-0.37323,21.3832,19.2497,18.5883,9.37448,19.6048,
188 11.5006,13.0159,14.5852,13.66,-1.04889};
189
190
191 /* Set up apm structure */
192 ap_t ap;
193 ap.lsiz = SIZE_IMAGE;
194 ap.csiz = SIZE_IMAGE;
195 ap.thresh = 11.0936;
196 ap.inframe = cpl_image_new(SIZE_IMAGE, SIZE_IMAGE, CPL_TYPE_DOUBLE);
197 ap.conframe = cpl_image_new(SIZE_IMAGE, SIZE_IMAGE, CPL_TYPE_DOUBLE);
198
199
200 /* Initialize */
201 hdrl_apinit(&ap);
202
203 ap.npl_pix = NT;
204 ap.plarray = cpl_realloc(ap.plarray, NT * sizeof(*ap.plarray));
205 for (cpl_size i = 0; i < NT; i++) {
206 ap.plarray[i].x = x[i];
207 ap.plarray[i].y = y[i];
208 ap.plarray[i].z = z[i];
209 ap.plarray[i].zsm = z[i];
210 }
211
212 ap.xintmin = 0.;
213 ap.areal_offset = 3.47165;
214 ap.thresh = 11.0936;
215 ap.fconst = 1.4427;
216
217 /* Run the test */
218 cpl_size iareal[NAREAL];
219 hdrl_areals(&ap, iareal);
220
221 ap.indata = cpl_image_get_data_double(ap.inframe);
222 ap.confdata = cpl_image_get_data_double(ap.conframe);
223 ap.mflag = cpl_calloc(2048 * 2048, sizeof(*ap.mflag));
224
225 /* Create a source */
226 double tmax = 1000.;
227 cpl_image_fill_gaussian(ap.inframe, CENTER_IMAGE, CENTER_IMAGE, tmax, 10., 10.);
228 cpl_image_fill_gaussian(ap.conframe, CENTER_IMAGE, CENTER_IMAGE, tmax, 100., 100.);
229
230 /* Do a basic moments analysis and work out the areal profiles*/
231 double momresults[8];
232 hdrl_moments(&ap, momresults);
233
234 double parmall[IMNUM][NPAR];
235 for (cpl_size i = 0; i < IMNUM; i++) {
236 for (cpl_size j = 0; j < NPAR; j++) {
237 parmall[i][j] = 0.;
238 }
239 }
240
241
242 /*** TESTS ***/
243 cpl_size nbit;
244 cpl_error_code e;
245
246 /* NOTE: Only if (iareal[0] >= ap->mulpix && ap->icrowd) */
247
248 /* Test 1 */
249 nbit = 1;
250 ap.ipnop = 1;
251 e = hdrl_overlp( &ap, parmall, &nbit,
252 momresults[1], momresults[2], momresults[3],
253 iareal[0], momresults[7]);
254 cpl_test_eq(e, CPL_ERROR_NONE);
255
256
257 /* Test 2 */
258 nbit = 1;
259 ap.ipnop = 1024;
260 e = hdrl_overlp( &ap, parmall, &nbit,
261 CENTER_IMAGE, CENTER_IMAGE, momresults[3],
262 SIZE_IMAGE * SIZE_IMAGE, tmax);
263 cpl_test_eq(e, CPL_ERROR_NONE);
264
265
266 nbit = 2;
267 ap.ipnop = 2;
268 ap.areal_offset = 1.5;
269 ap.thresh = 15.;
270 ap.fconst = 0.5;
271 hdrl_areals(&ap, iareal);
272 for (cpl_size i = 0; i < IMNUM; i++) {
273 for (cpl_size j = 0; j < NPAR; j++) {
274 parmall[i][j] = tmax;
275 }
276 }
277
278 /* Test 3 */
279 e = hdrl_overlp( &ap, parmall, &nbit,
280 CENTER_IMAGE, CENTER_IMAGE, momresults[3],
281 SIZE_IMAGE * SIZE_IMAGE, tmax);
282 cpl_test_eq(e, CPL_ERROR_NONE);
283
284 /* Test 4 */
285 e = hdrl_overlp_2orMore( &ap, parmall, &nbit,
286 CENTER_IMAGE, CENTER_IMAGE, momresults[3],
287 SIZE_IMAGE * SIZE_IMAGE, 1, 1, 1);
288 cpl_test_eq(e, CPL_ERROR_NONE);
289
290
291 /* Clean up */
292 hdrl_apclose(&ap);
293 cpl_free(ap.mflag);
294 cpl_image_delete(ap.inframe);
295 cpl_image_delete(ap.conframe);
296}
297
298
299int main(void)
300{
301 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
302
303 standard_test();
304 advanced_test();
305
306 return cpl_test_end(0);
307}
cpl_error_code hdrl_catalogue_conf(hdrl_casu_fits *infile, hdrl_casu_fits *conf, cpl_size ipix, double threshold, cpl_size icrowd, double rcore, cpl_size bkg_subtr, cpl_size nbsize, cpl_size cattype, double filtfwhm, double gain, double saturation, hdrl_casu_result *res)
Do source extraction.
void hdrl_apinit(ap_t *ap)
Initialize the ap structure.
void hdrl_apclose(ap_t *ap)
Close ap structure.
void hdrl_areals(ap_t *ap, cpl_size iareal[NAREAL])
Work out the areal profiles for an object.
hdrl_casu_fits * hdrl_casu_fits_wrap(cpl_image *im)
Wrap an image in a hdrl_casu_fits wrapper.
void hdrl_casu_tfits_delete(hdrl_casu_tfits *p)
Free all the workspace associated with a hdrl_casu_fits object.
void hdrl_casu_fits_delete(hdrl_casu_fits *p)
Free all the workspace associated with a hdrl_casu_fits object.
void hdrl_moments(ap_t *ap, double results[])
Do moments analysis on an object in a Plessey array.
cpl_error_code hdrl_overlp(ap_t *ap, double parm[IMNUM][NPAR], cpl_size *nbit, double xbar, double ybar, double total, cpl_size npix, double tmax)
Deblend overlapping images.
cpl_error_code hdrl_overlp_2orMore(ap_t *ap, double parm[IMNUM][NPAR], cpl_size *nbit, double xbar, double ybar, double total, cpl_size npix, double curthr_prev, double nexthr_prev, double lasthr_prev)
Deblend overlapping images 2 or more. If you know that exist more than 1 (come to the previous versio...