CR2RE Pipeline Reference Manual 1.7.0
hdrl_maglim-test.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2016 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
21#ifdef HAVE_CONFIG_H
22#include <config.h>
23#endif
24
25/*-----------------------------------------------------------------------------
26 Includes
27 -----------------------------------------------------------------------------*/
28
29#include "hdrl_maglim.h"
30
31#include <cpl.h>
32#include <stdint.h>
33#include <math.h>
34
35
36
37/*----------------------------------------------------------------------------*/
41/*----------------------------------------------------------------------------*/
42
51static cpl_image*
52hdrl_matrix_to_image_create(const cpl_matrix * kernel,
53 const cpl_boolean normalise) {
54
55 cpl_ensure(kernel != NULL, CPL_ERROR_NULL_INPUT, NULL);
56 cpl_size sx = cpl_matrix_get_ncol(kernel);
57 cpl_size sy = cpl_matrix_get_nrow(kernel);
58
59 cpl_image* image = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
60 double* pimage = cpl_image_get_data(image);
61 const double* pkernel = cpl_matrix_get_data_const(kernel);
62 cpl_size pix, piy;
63 for(cpl_size j = 0; j < sy; j++) {
64 piy = sx*j;
65 for(cpl_size i = 0; i < sx; i++) {
66 pix = piy + i;
67 pimage[pix] = pkernel[pix];
68 }
69 }
70 if(normalise) {
71 double sum = cpl_matrix_get_mean(kernel);
72 sum *= sx * sy;
73 cpl_image_divide_scalar(image, sum);
74 }
75 return image;
76
77}
78
79
80
81
82cpl_error_code
83test_hdrl_extend_image(void)
84{
85 cpl_image* ima;
86 cpl_image* ima1;
87 double* pinp;
88 const cpl_size border_nx = 3;
89 const cpl_size border_ny = 5;
90 hdrl_image_extend_method extend_method = HDRL_IMAGE_EXTEND_NEAREST;
91
92 cpl_size nx = 100;
93 cpl_size ny = 100;
94 ima = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
95 ima1 = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
96
97 for(cpl_size i = 1; i <= nx; i++) {
98 cpl_image_fill_window(ima1, i, 1, i, ny, i-1);
99 }
100 cpl_image_add(ima,ima1);
101
102 //cpl_image_save(ima, "ima1.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
103 for(cpl_size j = 1; j <= ny; j++) {
104 cpl_image_fill_window(ima1, 1, j, nx, j, (j-1) * 100);
105 }
106 //cpl_image_save(ima1, "ima2.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
107 cpl_image_add(ima,ima1);
108 //cpl_image_save(ima, "ima3.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
109 cpl_image_delete(ima1);
110
111 cpl_image* extended;
112
113 extended = hdrl_extend_image(NULL, border_nx, border_ny, extend_method);
114 cpl_test_error(CPL_ERROR_NULL_INPUT);
115 cpl_test_null(extended);
116
117 extended = hdrl_extend_image(ima, -1, border_ny, extend_method);
118 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
119 cpl_test_null(extended);
120
121 extended = hdrl_extend_image(ima, border_nx, -1, extend_method);
122 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
123 cpl_test_null(extended);
124
125 extended = hdrl_extend_image(ima, border_nx, border_ny, 3);
126 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
127 cpl_test_null(extended);
128
129 double* pext = cpl_image_get_data_double(ima);
130 for(cpl_size j = border_ny; j <= ny-border_ny; j++) {
131 for(cpl_size i = border_nx; i <= nx-border_nx; i++) {
132 cpl_test_abs(pext[j*nx+i], j*nx+i, HDRL_EPS_DATA);
133 }
134 }
135 extended = hdrl_extend_image(ima, border_nx, border_ny, extend_method);
136 cpl_test_error(CPL_ERROR_NONE);
137 cpl_test_nonnull(extended);
138
139 pext = cpl_image_get_data_double(extended);
140 pinp = cpl_image_get_data_double(ima);
141 /* the first (and last) border_nx X pixels and border_ny Y pixels of extended
142 * image have same value as 1st X pixel and Y pixel of input image
143 */
144 cpl_size sx = cpl_image_get_size_x(extended);
145 cpl_size sy = cpl_image_get_size_y(extended);
146 /* check central Y range */
147 for(cpl_size j = border_ny; j < sy-border_ny; j++) {
148 //cpl_msg_info(cpl_func,"j=%lld",j);
149 /* check early X range (i=0)*/
150 for(cpl_size i = 0; i < border_nx; i++) {
151 cpl_test_abs(pext[j*sx+i], pinp[(j-border_ny)*nx], HDRL_EPS_DATA);
152 }
153 // check end X range (i=nx-1)
154 for(cpl_size i = sx-border_nx; i < sx; i++) {
155 //cpl_msg_info(cpl_func,"i=%lld",i);
156 cpl_test_abs(pext[j*sx+i], pinp[(j-border_ny)*nx+(nx-1)], HDRL_EPS_DATA);
157 }
158 }
159
160
161 /* check early Y range */
162 for(cpl_size j = 0; j < border_ny; j++) {
163 // check central X range
164 for(cpl_size i = border_nx; i < sx-border_nx; i++) {
165 cpl_test_abs(pext[j*sx+i], pinp[0*nx+i-border_nx], HDRL_EPS_DATA);
166 } /* 0 is first j index of input image */
167 }
168
169
170 /* check end Y range */
171 for(cpl_size j = sy - border_ny; j < sy; j++) {
172 //cpl_msg_info(cpl_func,"j=%lld",j);
173 // check central X range
174 for(cpl_size i = border_nx; i < sx-border_nx; i++) {
175 //cpl_msg_info(cpl_func,"i=%lld",i);
176 cpl_test_abs(pext[j*sx+i], pinp[(ny-1)*nx+i-border_nx], HDRL_EPS_DATA);
177 } /* ny-1 is last j index of input image */
178 }
179
180 cpl_image_delete(extended);
181 cpl_image_delete(ima);
182 return cpl_error_get_code();
183}
184
185cpl_error_code
186test_hdrl_maglim_kernel_create(void)
187{
188 const cpl_size kernel_sx = 9;
189 const cpl_size kernel_sy = 9;
190 /* To simplify this unit test we choose
191 * fwhm_seeing = sigma_to_fwhm=sqrt(4*log(4))
192 * so that exponential: exp{-[(x^2+y^2)/2.0*(fwhm_seeing/sigma_to_fwhm)^2]}
193 * simplifies to exp[-(x^2+y^2)/2.0)]
194 * sqrt(4*log(4))=2,35482004503
195 */
196 const double fwhm_seeing = sqrt(4*log(4));
197 cpl_matrix* kernel_mat;
198 kernel_mat = hdrl_maglim_kernel_create(-1, kernel_sy, fwhm_seeing);
199 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
200 cpl_test_null(kernel_mat);
201
202 kernel_mat = hdrl_maglim_kernel_create(kernel_sx, -1, fwhm_seeing);
203 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
204 cpl_test_null(kernel_mat);
205
206 kernel_mat = hdrl_maglim_kernel_create(kernel_sx, kernel_sy, -1);
207 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
208 cpl_test_null(kernel_mat);
209
210 kernel_mat = hdrl_maglim_kernel_create(kernel_sx, kernel_sy, fwhm_seeing);
211 cpl_test_error(CPL_ERROR_NONE);
212 cpl_test_nonnull(kernel_mat);
213 double arg;
214 double x, y;
215
216 for(cpl_size j = 0; j < kernel_sy; j++) {
217 y = j - 0.5 * (kernel_sy-1);
218
219 for(cpl_size i = 0; i < kernel_sx; i++) {
220
221 x = i - 0.5 * (kernel_sx-1);
222 arg = (x * x + y * y) / 2;
223 cpl_test_abs(cpl_matrix_get(kernel_mat,j,i), exp(-arg), HDRL_EPS_DATA);
224
225 }
226 }
227
228 cpl_image*
229 kernel_ima = hdrl_matrix_to_image_create(kernel_mat, CPL_TRUE);
230 //cpl_image_save(kernel_ima, "kernel_ima.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
231 cpl_image_delete(kernel_ima);
232
233 cpl_matrix_delete(kernel_mat);
234 return cpl_error_get_code();
235}
236
237static cpl_image*
238test_util_crea_9x9_annular_image(void) {
239 /* We test with an image 9x9 with
240 * center at 100, then
241 * 1st annular at 90,
242 * 2nd annular at 80
243 */
244 cpl_image* image = cpl_image_new(9, 9, CPL_TYPE_DOUBLE);
245 cpl_image_add_scalar(image, 60);
246 cpl_image_fill_window(image, 2, 2, 8, 8, 70);
247 cpl_image_fill_window(image, 3, 3, 7, 7, 80);
248 cpl_image_fill_window(image, 4, 4, 6, 6, 90);
249 cpl_image_fill_window(image, 5, 5, 5, 5, 100);
250
251 return image;
252}
253cpl_error_code
254test_hdrl_image_convolve(void)
255{
256 cpl_image * input_image ;
257 cpl_matrix * kernel;
258 cpl_image * kernel_image;
259 hdrl_image_extend_method method = HDRL_IMAGE_EXTEND_MIRROR;
260 cpl_image* convolved;
261 const cpl_size kernel_sx = 9;
262 const cpl_size kernel_sy = 9;
263 const double fwhm_seeing = sqrt(4*log(4));
264 const int sx = 100;
265 const int sy = 100;
266
267 input_image = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
268 cpl_image_add_scalar(input_image,1);
269
270 kernel = hdrl_maglim_kernel_create(kernel_sx, kernel_sy, fwhm_seeing);
271 kernel_image = hdrl_matrix_to_image_create(kernel, CPL_TRUE);
272 //cpl_image_save(kernel_image, "kernel_image.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
273 cpl_image_delete(kernel_image);
274
275 convolved = hdrl_image_convolve(NULL, kernel, method);
276 cpl_test_error(CPL_ERROR_NULL_INPUT);
277 cpl_test_null(convolved);
278
279 convolved = hdrl_image_convolve(input_image, NULL, method);
280 cpl_test_error(CPL_ERROR_NULL_INPUT);
281 cpl_test_null(convolved);
282
283 convolved = hdrl_image_convolve(input_image, kernel, 2);
284 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
285 cpl_test_null(convolved);
286
287 //cpl_image_save(input_image, "input_image.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
288 convolved = hdrl_image_convolve(input_image, kernel, method);
289 //cpl_image_save(convolved, "convolved.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
290 cpl_test_error(CPL_ERROR_NONE);
291 cpl_test_nonnull(convolved);
292 cpl_image_delete(input_image);
293 cpl_matrix_delete(kernel);
294 cpl_image_delete(convolved);
295
296 /*
297 * Kernel 3x3, central pixel 1.
298 */
299
300 input_image = test_util_crea_9x9_annular_image();
301 //cpl_image_save(input_image, "convolved.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
302
303 kernel = cpl_matrix_new(3, 3);
304 cpl_matrix_set(kernel, 1, 1, 1);
305 cpl_matrix_set(kernel, 1, 0, 0.7);
306 cpl_matrix_set(kernel, 0, 1, 0.7);
307 cpl_matrix_set(kernel, 2, 1, 0.7);
308 cpl_matrix_set(kernel, 1, 2, 0.7);
309 kernel_image = hdrl_matrix_to_image_create(kernel, CPL_FALSE);
310 //cpl_image_save(kernel_image, "kernel_image.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
311 cpl_image_delete(kernel_image);
312
313 convolved = hdrl_image_convolve(input_image, kernel, method);
314 //cpl_image_save(convolved, "convolved.fits", CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
315
316 int status = 0;
317 //cpl_msg_info(cpl_func,"value: %16.16g",cpl_image_get(convolved, 4, 4, &status));
318 /* test some values */
319 cpl_test_abs(cpl_image_get(convolved, 5, 5, &status),92.63157894736842, HDRL_EPS_DATA);
320 cpl_test_abs(cpl_image_get(convolved, 4, 5, &status),90.0, HDRL_EPS_DATA);
321 cpl_test_abs(cpl_image_get(convolved, 3, 5, &status),80.0, HDRL_EPS_DATA);
322 cpl_test_abs(cpl_image_get(convolved, 2, 5, &status),70.0, HDRL_EPS_DATA);
323 cpl_test_abs(cpl_image_get(convolved, 1, 5, &status),61.8421052631579, HDRL_EPS_DATA);
324 cpl_test_abs(cpl_image_get(convolved, 1, 1, &status),60.0, HDRL_EPS_DATA);
325 cpl_test_abs(cpl_image_get(convolved, 2, 2, &status),66.31578947368422, HDRL_EPS_DATA);
326 cpl_test_abs(cpl_image_get(convolved, 3, 3, &status),76.31578947368422, HDRL_EPS_DATA);
327 cpl_test_abs(cpl_image_get(convolved, 4, 4, &status),86.31578947368422, HDRL_EPS_DATA);
328 cpl_image_delete(convolved);
329 cpl_image_delete(input_image);
330 cpl_matrix_delete(kernel);
331 return cpl_error_get_code();
332}
333
334cpl_error_code
335test_hdrl_matrix_to_image_create(void)
336{
337 cpl_matrix* kernel;
338 cpl_image* image;
339 const cpl_size kernel_sx = 9;
340 const cpl_size kernel_sy = 9;
341 const double fwhm_seeing = sqrt(4*log(4));
342
343
344 kernel = hdrl_maglim_kernel_create(kernel_sx, kernel_sy, fwhm_seeing);
345
346 image = hdrl_matrix_to_image_create(NULL, CPL_TRUE);
347 cpl_test_error(CPL_ERROR_NULL_INPUT);
348 cpl_test_null(image);
349
350 image = hdrl_matrix_to_image_create(kernel, CPL_TRUE);
351 cpl_test_nonnull(image);
352
353 cpl_test_abs(cpl_image_get_size_x(image), kernel_sx, HDRL_EPS_DATA);
354 cpl_test_abs(cpl_image_get_size_y(image), kernel_sy, HDRL_EPS_DATA);
355
356 cpl_image_delete(image);
357 cpl_matrix_delete(kernel);
358 return cpl_error_get_code();
359}
360
361cpl_error_code
362test_hdrl_maglim_compute(void) {
363 cpl_image* image;
364 double zeropoint = 0;
365 double fwhm_seeing = sqrt(4*log(4));;
366 const cpl_size kernel_sx = 9;
367 const cpl_size kernel_sy = 9;
368 double limiting_magnitude = 0;
369
370 hdrl_image_extend_method method = HDRL_IMAGE_EXTEND_MIRROR;
371 double histo_min = 0.;
372 double histo_max = 0.;
373 double bin_size = 0.;
374 cpl_size error_niter = 0;
375 hdrl_mode_type mode_method = HDRL_MODE_MEDIAN;
376 hdrl_parameter * mode_parameter = hdrl_collapse_mode_parameter_create(histo_min,
377 histo_max,
378 bin_size,
379 mode_method,
380 error_niter);
381 /* testinvalid input */
382 hdrl_maglim_compute(NULL, zeropoint, fwhm_seeing, kernel_sx, kernel_sy,
383 method, mode_parameter, &limiting_magnitude);
384 cpl_test_error(CPL_ERROR_NULL_INPUT);
385
386 image = test_util_crea_9x9_annular_image();
387 hdrl_maglim_compute(image, zeropoint, -1, kernel_sx, kernel_sy,
388 method, mode_parameter, &limiting_magnitude);
389 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
390
391
392 hdrl_maglim_compute(image, zeropoint, fwhm_seeing, -1, kernel_sy,
393 method, mode_parameter, &limiting_magnitude);
394 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
395
396 hdrl_maglim_compute(image, zeropoint, fwhm_seeing, kernel_sx, -1,
397 method, mode_parameter, &limiting_magnitude);
398 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
399
400
401 hdrl_maglim_compute(image, zeropoint, fwhm_seeing, kernel_sx, kernel_sy,
402 -1, mode_parameter, &limiting_magnitude);
403 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
404
405 hdrl_maglim_compute(image, zeropoint, fwhm_seeing, kernel_sx, kernel_sy,
406 method, NULL, &limiting_magnitude);
407 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
408
409
410 hdrl_maglim_compute(image, zeropoint, fwhm_seeing, kernel_sx, kernel_sy,
411 method, mode_parameter, &limiting_magnitude);
412
413
414 cpl_test_abs(limiting_magnitude,-5.591854160255954,HDRL_EPS_DATA);
415
416
417 cpl_image_delete(image);
418 hdrl_parameter_delete(mode_parameter);
419 return cpl_error_get_code();
420}
421
422
423/*----------------------------------------------------------------------------*/
427/*----------------------------------------------------------------------------*/
428int main(void)
429{
430 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_INFO);
431
432 test_hdrl_extend_image();
433 test_hdrl_maglim_kernel_create();
434 test_hdrl_image_convolve();
435 test_hdrl_matrix_to_image_create();
436 test_hdrl_maglim_compute();
437
438 return cpl_test_end(0);
439}
hdrl_collapse_mode_parameter_create
hdrl_parameter * hdrl_collapse_mode_parameter_create(double histo_min, double histo_max, double bin_size, hdrl_mode_type mode_method, cpl_size error_niter)
create a parameter object for the mode
Definition hdrl_collapse.c:567
hdrl_maglim_compute
cpl_error_code hdrl_maglim_compute(const cpl_image *image, const double zeropoint, const double fwhm, const cpl_size kernel_size_x, const cpl_size kernel_size_y, const hdrl_image_extend_method image_extend_method, const hdrl_parameter *mode_parameter, double *limiting_magnitude)
Computes the limiting magnitude of an image.
Definition hdrl_maglim.c:114
hdrl_parameter_delete
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
Definition hdrl_parameter.c:148
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2